mod_DFT_FMT_d.F90

Go to the documentation of this file.

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
!
!All equations are taken from the appendix of:
!Gross: 'A density functional theory for vapor-liquid interfaces using the PCP-SAFT eos'
!But: here we dont treat chain-molecules!! -> no multiplications with segment number
MODULE mod_dft_fmt_d
  IMPLICIT NONE
  PRIVATE 
  PUBLIC fmt_weighted_densities_d
  PUBLIC fmt_dfdrho_d

CONTAINS
!  Differentiation of fmt_weighted_densities in forward (tangent) mode:
!   variations   of useful results: phi_dn0 phi_dn1 phi_dn2 phi_dnv1
!                phi_dn3 phi_dnv2
!   with respect to varying inputs: rhop
  SUBROUTINE fmt_weighted_densities_d(rhop, rhopd, n0, n1, n2, n3, nv1, &
&   nv2, phi_dn0, phi_dn0d, phi_dn1, phi_dn1d, phi_dn2, phi_dn2d, &
&   phi_dn3, phi_dn3d, phi_dnv1, phi_dnv1d, phi_dnv2, phi_dnv2d, user)
    USE parameters, ONLY : pi
    USE basic_variables, ONLY : ncomp,parame
    USE eos_variables, Only: dhs
    USE mod_dft, ONLY : zp, dzp, fa


!PETSc module
Use petscmanagement
IMPLICIT NONE

#include <finclude/petscsys.h>

!passed
type(userctx)                                         :: user
petscscalar                                            :: rhop(ncomp,user%gxs:user%gxe)
petscscalar                                            :: rhopd(ncomp,user%gxs:user%gxe)
REAL,dimension(user%gxs:user%gxe),Intent(OUT)          :: n0,n1,n2,n3,nv1,nv2  !ngp muss groesser als fa+fa/2 sein!!
REAL,dimension(user%gxs:user%gxe),Intent(OUT)          :: phi_dn0,phi_dn1,phi_dn2,phi_dn3,phi_dnv1,phi_dnv2
REAL,dimension(user%gxs:user%gxe),Intent(OUT)          :: phi_dn0d,phi_dn1d,phi_dn2d,phi_dn3d,phi_dnv1d,phi_dnv2d



    
!local
    REAL, DIMENSION(user%gxs:user%gxe) :: n0d,n1d,n2d,n3d,nv1d,nv2d
    INTEGER :: k, i, j
    INTEGER :: fa2
    REAL :: dz, d2, zz
    INTEGER :: n
    INTEGER, parameter :: nmax = 800
    REAL, DIMENSION(NMAX) :: x_int, n2_int, n3_int, nv2_int
    REAL, DIMENSION(NMAX) :: n2_intd, n3_intd, nv2_intd
    REAL, DIMENSION(NMAX) :: y2_n2, y2_n3, y2_nv2
    REAL, DIMENSION(NMAX) :: y2_n2d, y2_n3d, y2_nv2d
    REAL :: int_n2, int_n3, int_nv2, xhi, xlo
    REAL :: int_n2d, int_n3d, int_nv2d
    REAL :: zms, zms2, zms3, logzms
    REAL :: zmsd, zms2d, zms3d, logzmsd
    REAL :: nn0, nn1, nn2, nn3, nnv1, nnv2
    REAL :: nn0d, nn1d, nn2d, nn3d, nnv1d, nnv2d
    REAL :: rhopjk, rhopjp1k
    REAL :: rhopjkd, rhopjp1kd
    INTRINSIC maxval
    INTRINSIC log
    REAL :: result1
    


    n0 = 0.0
    n1 = 0.0
    n2 = 0.0
    n3 = 0.0
    nv1 = 0.0
    nv2 = 0.0
    result1 = maxval(fa(1:ncomp) + 5)
    fa2 = result1/2
    nv1d = 0.0
    nv2d = 0.0
    n0d = 0.0
    n1d = 0.0
    n2d = 0.0
    n3d = 0.0
    y2_n2d = 0.0
    y2_n3d = 0.0
    y2_nv2d = 0.0
    DO k=1,ncomp
!fa2 = (fa(k) + 5) / 2   !grid points in sig/2
!half of dhs [A] 
      d2 = dhs(k)/2.0




       Do i=user%xs-fa2,user%xe+fa2 !to evaluate dF/drho at any point, the derivatives dphi/dn have to be evaluated at +-d/2 around this point
                 
        n = 1
        x_int = 0.0
        n2_int = 0.0
        n3_int = 0.0
        nv2_int = 0.0
        n2_intd = 0.0
        nv2_intd = 0.0
        n3_intd = 0.0
!to evaluate dphi/dn at a given point, the weighted densities are needed at +-d/2 around this point
        DO j=i-fa2,i+fa2
          rhopjkd = parame(k,1)*rhopd(k, j)
          rhopjk = rhop(k, j)*parame(k,1)
          rhopjp1kd = parame(k,1)*rhopd(k, j+1)
          rhopjp1k = rhop(k, j+1)*parame(k,1)
          IF (zp(i) - zp(j+1) .LT. d2 .AND. zp(i) - zp(j) .GE. d2) THEN
!the position of j+1 is already within i-d/2 while j is still outside this range in this case, the integration steplength (dz) is
! just the distance, which j+1 overlaps with i-d/2 and what is integrated is the interpolated value of the integrand
!here always n=1!
            IF (n .NE. 1) THEN
              GOTO 100
            ELSE
!distance between grid points j and i
              zz = zp(j) - zp(i)
!the part of the intervall between zp(j) and zp(j+1) which is already within i-d/2
              dz = zp(j+1) - (zp(i)-d2)
!if(dz < epsilon(dz)) dz = epsilon(dz) !bei unguenstiger Kombination von sig und ngrid kann dz unter Machinengenauigkeit epsilon 
!liegen, dann ist x(2) = x(1) + dz = x(1) -> das fuehrt zu Abbruch in Spline Interpolation           
!array containing x-values for spline integration 
              x_int(n) = 0.0
!integrand für n2 (=rhop) linear interpoliert für den Punk zp(i)-d/2             
              n2_intd(n) = rhopjkd + (dzp-dz)*(rhopjp1kd-rhopjkd)/dzp
              n2_int(n) = rhopjk + (rhopjp1k-rhopjk)/dzp*(dzp-dz)
!integrand für n3: rhop*(d2**2 - z'**2), da hier gerade z' = d2 -> integrand wird hier = 0!! 
              n3_intd(n) = 0.0
              n3_int(n) = 0.0
!analog
              nv2_intd(n) = zz*rhopjkd + (dzp-dz)*((zp(j+1)-zp(i))*&
&               rhopjp1kd-zz*rhopjkd)/dzp
              nv2_int(n) = rhopjk*zz + (rhopjp1k*(zp(j+1)-zp(i))-rhopjk*&
&               zz)/dzp*(dzp-dz)
            END IF
          ELSE IF (zp(j) .GT. zp(i) - d2 .AND. zp(j) .LE. zp(i) + d2) &
&         THEN
!grid point j  within i+-d2 
            n = n + 1
!first time in this If condition, dz is stil the old value from above!
            x_int(n) = x_int(n-1) + dz
            zz = zp(j) - zp(i)
            dz = dzp
            n2_intd(n) = rhopjkd
            n2_int(n) = rhopjk
            n3_intd(n) = (d2**2-zz**2)*rhopjkd
            n3_int(n) = rhopjk*(d2**2-zz**2)
            nv2_intd(n) = zz*rhopjkd
            nv2_int(n) = rhopjk*zz
            IF (zp(j) .LT. zp(i) + d2 .AND. zp(j+1) .GE. zp(i) + d2) &
&           THEN
!zp(j) is still within zp(i)+d2 but zp(j+1) is already outside zp(i)+d2
              dz = zp(i) + d2 - zp(j)
!If(dz <= epsilon(dz)) exit   !wie oben, kann auch hier bei ungluecklicher Wahl von sig und ngrid dz < eps werden und somit x(n) 
!= x(n-1) -> Abbruch in Spline interpolation. Dann einfach ngrid aendern!
              zz = zp(j) - zp(i)
              n = n + 1
              x_int(n) = x_int(n-1) + dz
!            If(x_int(n) == x_int(n-1)) Then
!              n = n - 1
!              exit
!            End If  
              n2_intd(n) = rhopjkd + dz*(rhopjp1kd-rhopjkd)/dzp
              n2_int(n) = rhopjk + (rhopjp1k-rhopjk)/dzp*dz
!Begründung wie oben
              n3_intd(n) = 0.0
              n3_int(n) = 0.0
              nv2_intd(n) = zz*rhopjkd + dz*((zp(j+1)-zp(i))*rhopjp1kd-&
&               zz*rhopjkd)/dzp
              nv2_int(n) = rhopjk*zz + (rhopjp1k*(zp(j+1)-zp(i))-rhopjk*&
&               zz)/dzp*dz
            END IF
          END IF
        END DO
!spline integration
        xlo = x_int(1)
        xhi = x_int(n)
        CALL spline_d(x_int, n2_int, n2_intd, n, 1.e30, 1.e30, y2_n2, &
&               y2_n2d)
        CALL spline_d(x_int, n3_int, n3_intd, n, 1.e30, 1.e30, y2_n3, &
&               y2_n3d)
        CALL spline_d(x_int, nv2_int, nv2_intd, n, 1.e30, 1.e30, y2_nv2&
&               , y2_nv2d)
        CALL splint_integral_d(x_int, n2_int, n2_intd, y2_n2, y2_n2d, n&
&                        , xlo, xhi, int_n2, int_n2d)
        CALL splint_integral_d(x_int, n3_int, n3_intd, y2_n3, y2_n3d, n&
&                        , xlo, xhi, int_n3, int_n3d)
        CALL splint_integral_d(x_int, nv2_int, nv2_intd, y2_nv2, y2_nv2d&
&                        , n, xlo, xhi, int_nv2, int_nv2d)
!weighted densities
        n2d(i) = n2d(i) + pi*dhs(k)*int_n2d
        n2(i) = n2(i) + pi*dhs(k)*int_n2
        n1d(i) = n1d(i) + 0.5*int_n2d
        n1(i) = n1(i) + 0.5*int_n2
        n0d(i) = n0d(i) + int_n2d/dhs(k)
        n0(i) = n0(i) + int_n2/dhs(k)
        n3d(i) = n3d(i) + pi*int_n3d
        n3(i) = n3(i) + pi*int_n3
        nv2d(i) = nv2d(i) - 2.0*pi*int_nv2d
        nv2(i) = nv2(i) - 2.0*pi*int_nv2
        nv1d(i) = nv1d(i) - int_nv2d/dhs(k)
        nv1(i) = nv1(i) - int_nv2/dhs(k)
      END DO
    END DO
    phi_dn0d = 0.0
    phi_dn1d = 0.0
    phi_dn2d = 0.0
    phi_dnv1d = 0.0
    phi_dn3d = 0.0
    phi_dnv2d = 0.0
!     If(i < 5) Then 
!     write(*,*)'i dens',i,nv2(i),n1(i),n0(i)
!     write(*,*)'i dens',i,n3(i),nv1(i),nv2(i)
!     end if
!    
!    If(i > 95) then 
!       write(*,*)'i dens',i,nv2(i),n1(i),n0(i)
!       write(*,*)'i dens',i,n3(i),nv1(i),nv2(i)
!    End If 
!    
! pause


!derivatives of FMT helmholtz energy density w.r.t. weighted densities 
Do i=user%xs-maxval((fa(1:ncomp)+1)/2),user%xe+maxval((fa(1:ncomp)+1)/2)     
!weils kürzer ist
      nn0d = n0d(i)
      nn0 = n0(i)
      nn1d = n1d(i)
      nn1 = n1(i)
      nn2d = n2d(i)
      nn2 = n2(i)
      nn3d = n3d(i)
      nn3 = n3(i)
      nnv1d = nv1d(i)
      nnv1 = nv1(i)
      nnv2d = nv2d(i)
      nnv2 = nv2(i)
      zmsd = -nn3d
      zms = 1.0 - nn3
      zms2d = zmsd*zms + zms*zmsd
      zms2 = zms*zms
      zms3d = zms2d*zms + zms2*zmsd
      zms3 = zms2*zms
      logzmsd = zmsd/zms
      logzms = log(zms)
!if(isnan(logzms)) stop 'zms < 0, log(zms) undefined FMT_Weighted_Densities'    
      phi_dn0d(i) = -logzmsd
      phi_dn0(i) = -logzms
      phi_dn1d(i) = (nn2d*zms-nn2*zmsd)/zms**2
      phi_dn1(i) = nn2/zms
      phi_dn2d(i) = (nn1d*zms-nn1*zmsd)/zms**2 + (3.0*((nn2d*nn2+nn2*&
&       nn2d-nnv2d*nnv2-nnv2*nnv2d)*(nn3+zms2*logzms)+(nn2*nn2-nnv2*nnv2&
&       )*(nn3d+zms2d*logzms+zms2*logzmsd))*36.0*pi*nn3**2*zms2-3.0*(nn2&
&       *nn2-nnv2*nnv2)*(nn3+zms2*logzms)*36.0*pi*((nn3d*nn3+nn3*nn3d)*&
&       zms2+nn3**2*zms2d))/(36.0*pi*nn3*nn3*zms2)**2
      phi_dn2(i) = nn1/zms + 3.0*(nn2*nn2-nnv2*nnv2)*(nn3+zms2*logzms)/(&
&       36.0*pi*nn3*nn3*zms2)
      phi_dn3d(i) = (nn0d*zms-nn0*zmsd)/zms**2 + ((nn1d*nn2+nn1*nn2d-&
&       nnv1d*nnv2-nnv1*nnv2d)*zms2-(nn1*nn2-nnv1*nnv2)*zms2d)/zms2**2 -&
&       (((3*nn2**2*nn2d-3.0*((nn2d*nnv2+nn2*nnv2d)*nnv2+nn2*nnv2*nnv2d)&
&       )*(nn3*(nn3**2-5.0*nn3+2.0)+2.0*zms3*logzms)+(nn2**3-3.0*nn2*&
&       nnv2*nnv2)*(nn3d*(nn3**2-5.0*nn3+2.0)+nn3*(2*nn3*nn3d-5.0*nn3d)+&
&       2.0*(zms3d*logzms+zms3*logzmsd)))*36.0*pi*nn3**3*zms3-(nn2**3-&
&       3.0*nn2*nnv2*nnv2)*(nn3*(nn3**2-5.0*nn3+2.0)+2.0*zms3*logzms)*&
&       36.0*pi*(3*nn3**2*nn3d*zms3+nn3**3*zms3d))/(36.0*pi*nn3**3*zms3)&
&       **2
      phi_dn3(i) = nn0/zms + (nn1*nn2-nnv1*nnv2)/zms2 - (nn2**3-3.0*nn2*&
&       nnv2*nnv2)*(nn3*(nn3**2-5.0*nn3+2.0)+2.0*zms3*logzms)/(36.0*pi*&
&       nn3**3*zms3)
      phi_dnv1d(i) = -((nnv2d*zms-nnv2*zmsd)/zms**2)
      phi_dnv1(i) = -(nnv2/zms)
      phi_dnv2d(i) = -((nnv1d*zms-nnv1*zmsd)/zms**2) - (6.0*((nn2d*nnv2+&
&       nn2*nnv2d)*(nn3+zms2*logzms)+nn2*nnv2*(nn3d+zms2d*logzms+zms2*&
&       logzmsd))*36.0*pi*nn3**2*zms2-6.0*nn2*nnv2*(nn3+zms2*logzms)*&
&       36.0*pi*(2*nn3*nn3d*zms2+nn3**2*zms2d))/(36.0*pi*nn3**2*zms2)**2
      phi_dnv2(i) = -(nnv1/zms) - 6.0*nn2*nnv2*(nn3+zms2*logzms)/(36.0*&
&       pi*nn3**2*zms2)
    END DO
    GOTO 110
 100 stop
 110 CONTINUE
  END SUBROUTINE fmt_weighted_densities_d



!  Differentiation of fmt_dfdrho in forward (tangent) mode:
!   variations   of useful results: df_drho_fmt
!   with respect to varying inputs: df_drho_fmt phi_dn0 phi_dn1
!                phi_dn2 phi_dnv1 phi_dn3 phi_dnv2
  SUBROUTINE fmt_dfdrho_d(i, fa, user, phi_dn0, phi_dn0d, phi_dn1, &
&   phi_dn1d, phi_dn2, phi_dn2d, phi_dn3, phi_dn3d, phi_dnv1, phi_dnv1d&
&   , phi_dnv2, phi_dnv2d, df_drho_fmt, df_drho_fmtd)
    USE basic_variables, ONLY : ncomp,parame
    USE eos_variables, Only: dhs
    USE mod_dft, ONLY : zp, dzp
    USE parameters, ONLY : pi
   

!PETSc module
Use petscmanagement

 IMPLICIT NONE

!passed
INTEGER, INTENT(IN) :: i         !the grid point at which to calculate dFdrho
INTEGER, INTENT(IN) :: fa(ncomp)
type(userctx)                                     :: user
REAL,dimension(user%gxs:user%gxe),INTENT(IN)       :: phi_dn0,phi_dn1,phi_dn2,phi_dn3,phi_dnv1,phi_dnv2
REAL,dimension(user%gxs:user%gxe),INTENT(IN)       :: phi_dn0d,phi_dn1d,phi_dn2d,phi_dn3d,phi_dnv1d,phi_dnv2d

REAL,dimension(user%xs:user%xe,ncomp), INTENT(OUT)       :: df_drho_fmt
REAL,dimension(user%xs:user%xe,ncomp), INTENT(OUT)       :: df_drho_fmtd


!local
    INTEGER :: j, k, n
    INTEGER :: fa2
    REAL :: dz, d2, zz, zz_jp1
    INTEGER, parameter :: nmax = 800
    REAL, DIMENSION(NMAX) :: x_int, y_int, y0_int, y1_int, y2_int, y3_int&
&   , yv1_int, yv2_int, y2
    REAL, DIMENSION(NMAX) :: y_intd, y2d
    REAL :: xhi, xlo, integral, int0, int1, int2, int3, intv1, intv2
    REAL :: integrald
    REAL :: at_j, at_jp1
    REAL :: at_jd, at_jp1d


    y2d = 0.0
    y_intd = 0.0
!Das Integral (Gleichung A1 in Gross DFT 2009) wird hier in einem Schlag berechnet!
!Falls die einzelnen Terme einzeln integriert werden sollen, einfach die auskommentierte Version verwenden
!das einzige, das hier von k abhaengt, sind fa und dhs!! die Ableitungen phi_dn... sind nicht Komponentenspez, da die 
    DO k=1,ncomp
!gewichteten Dichten ja uch nicht mehr Komponentenspez sind (n_i = sum n_i(k))
      n = 1
!number of grid points within dhs/2
      fa2 = (fa(k)+5)/2
      d2 = dhs(k)/2.0
      DO j=i-fa2,i+fa2
        IF (zp(i) - zp(j+1) .LT. d2 .AND. zp(i) - zp(j) .GE. d2) THEN
!the position of j+1 is already within i-d/2 while j is still outside this range in this case, the integration steplength (dz) is
! just the distance, which j+1 overlaps with i-d/2 and what is integrated is the interpolated value of the integrand
          x_int(n) = 0.0
          IF (n .NE. 1) THEN
            GOTO 100
          ELSE
            dz = zp(j+1) - (zp(i)-d2)
            zz = zp(j) - zp(i)
            zz_jp1 = zp(j+1) - zp(i)
            at_jd = phi_dn0d(j)/dhs(k) + 0.5*phi_dn1d(j) + pi*dhs(k)*&
&             phi_dn2d(j) + pi*(d2**2-zz**2)*phi_dn3d(j) + zz*phi_dnv1d(&
&             j)/dhs(k) + 2.0*pi*zz*phi_dnv2d(j)
            at_j = phi_dn0(j)/dhs(k) + 0.5*phi_dn1(j) + pi*dhs(k)*&
&             phi_dn2(j) + pi*phi_dn3(j)*(d2**2-zz**2) + phi_dnv1(j)*zz/&
&             dhs(k) + 2.0*pi*phi_dnv2(j)*zz
            at_jp1d = phi_dn0d(j+1)/dhs(k) + 0.5*phi_dn1d(j+1) + pi*dhs(&
&             k)*phi_dn2d(j+1) + pi*(d2**2-zz_jp1**2)*phi_dn3d(j+1) + &
&             zz_jp1*phi_dnv1d(j+1)/dhs(k) + 2.0*pi*zz_jp1*phi_dnv2d(j+1&
&             )
            at_jp1 = phi_dn0(j+1)/dhs(k) + 0.5*phi_dn1(j+1) + pi*dhs(k)*&
&             phi_dn2(j+1) + pi*phi_dn3(j+1)*(d2**2-zz_jp1**2) + &
&             phi_dnv1(j+1)*zz_jp1/dhs(k) + 2.0*pi*phi_dnv2(j+1)*zz_jp1
!lineare interpolation genau, wie in FMT_Weighted_Densities         
            y_intd(n) = at_jd + (dzp-dz)*(at_jp1d-at_jd)/dzp
            y_int(n) = at_j + (at_jp1-at_j)/dzp*(dzp-dz)
!           y0_int(n) = phi_dn0(j) + (phi_dn0(j+1) -phi_dn0(j))/dzp * (dzp-dz)
!           y1_int(n) = phi_dn1(j) + (phi_dn1(j+1) -phi_dn1(j))/dzp * (dzp-dz)          
!           y2_int(n) = phi_dn2(j) + (phi_dn2(j+1) -phi_dn2(j))/dzp * (dzp-dz)          
!           y3_int(n) = phi_dn3(j)*(d2**2-zz**2) + (phi_dn3(j+1)*(d2**2-zz_jp1**2) - phi_dn3(j)*(d2**2-zz**2) )/dzp * (dzp-dz)
!           yv1_int(n) = phi_dnv1(j)*zz + (phi_dnv1(j+1)*zz_jp1 - phi_dnv1(j)*zz)/dzp * (dzp-dz)          
!           yv2_int(n) = phi_dnv2(j)*zz + (phi_dnv2(j+1)*zz_jp1 - phi_dnv2(j)*zz)/dzp * (dzp-dz)
!           
          END IF
        ELSE IF (zp(j) .GT. zp(i) - d2 .AND. zp(j) .LE. zp(i) + d2) THEN
!grid points j and j+1 are completely within i+-d2 
          n = n + 1
          zz = zp(j) - zp(i)
          x_int(n) = x_int(n-1) + dz
          y_intd(n) = phi_dn0d(j)/dhs(k) + 0.5*phi_dn1d(j) + pi*dhs(k)*&
&           phi_dn2d(j) + pi*(d2**2-zz**2)*phi_dn3d(j) + zz*phi_dnv1d(j)&
&           /dhs(k) + 2.0*pi*zz*phi_dnv2d(j)
          y_int(n) = phi_dn0(j)/dhs(k) + 0.5*phi_dn1(j) + pi*dhs(k)*&
&           phi_dn2(j) + pi*phi_dn3(j)*(d2**2-zz**2) + phi_dnv1(j)*zz/&
&           dhs(k) + 2.0*pi*phi_dnv2(j)*zz
!           y0_int(n) = phi_dn0(j)
!           y1_int(n) = phi_dn1(j)
!           y2_int(n) = phi_dn2(j)
!           y3_int(n) = phi_dn3(j)*(d2**2-zz**2)
!           yv1_int(n) = phi_dnv1(j)*zz
!           yv2_int(n) = phi_dnv2(j)*zz
!           
          dz = dzp
          IF (zp(j) .LT. zp(i) + d2 .AND. zp(j+1) .GT. zp(i) + d2) THEN
!zp(j) is still within zp(i)+d2 but zp(j+1) is already out side zp(i)+d2
            n = n + 1
            zz = zp(j) - zp(i)
            zz_jp1 = zp(j+1) - zp(i)
            dz = zp(i) + d2 - zp(j)
            x_int(n) = x_int(n-1) + dz
            at_jd = phi_dn0d(j)/dhs(k) + 0.5*phi_dn1d(j) + pi*dhs(k)*&
&             phi_dn2d(j) + pi*(d2**2-zz**2)*phi_dn3d(j) + zz*phi_dnv1d(&
&             j)/dhs(k) + 2.0*pi*zz*phi_dnv2d(j)
            at_j = phi_dn0(j)/dhs(k) + 0.5*phi_dn1(j) + pi*dhs(k)*&
&             phi_dn2(j) + pi*phi_dn3(j)*(d2**2-zz**2) + phi_dnv1(j)*zz/&
&             dhs(k) + 2.0*pi*phi_dnv2(j)*zz
            at_jp1d = phi_dn0d(j+1)/dhs(k) + 0.5*phi_dn1d(j+1) + pi*dhs(&
&             k)*phi_dn2d(j+1) + pi*(d2**2-zz_jp1**2)*phi_dn3d(j+1) + &
&             zz_jp1*phi_dnv1d(j+1)/dhs(k) + 2.0*pi*zz_jp1*phi_dnv2d(j+1&
&             )
            at_jp1 = phi_dn0(j+1)/dhs(k) + 0.5*phi_dn1(j+1) + pi*dhs(k)*&
&             phi_dn2(j+1) + pi*phi_dn3(j+1)*(d2**2-zz_jp1**2) + &
&             phi_dnv1(j+1)*zz_jp1/dhs(k) + 2.0*pi*phi_dnv2(j+1)*zz_jp1
            y_intd(n) = at_jd + dz*(at_jp1d-at_jd)/dzp
            y_int(n) = at_j + (at_jp1-at_j)/dzp*dz
!           y0_int(n) = phi_dn0(j) + (phi_dn0(j+1) - phi_dn0(j))/dzp * dz
!           y1_int(n) = phi_dn1(j) + (phi_dn1(j+1) - phi_dn1(j))/dzp * dz
!           y2_int(n) = phi_dn2(j) + (phi_dn2(j+1) - phi_dn2(j))/dzp * dz
!           y3_int(n) = phi_dn3(j)*(d2**2 - zz**2) + (phi_dn3(j+1)*(d2**2 - zz_jp1**2) - phi_dn3(j)*(d2**2 - zz**2))/dzp * dz
!           yv1_int(n) = phi_dnv1(j)*zz + (phi_dnv1(j+1)*zz_jp1 - phi_dnv1(j)*zz)/dzp * dz  
!           yv2_int(n) = phi_dnv2(j)*zz + (phi_dnv2(j+1)*zz_jp1 - phi_dnv2(j)*zz)/dzp * dz  
!           
          END IF
        END IF
      END DO
!spline integration
      xlo = x_int(1)
      xhi = x_int(n)
      CALL spline_d(x_int, y_int, y_intd, n, 1.e30, 1.e30, y2, y2d)
      CALL splint_integral_d(x_int, y_int, y_intd, y2, y2d, n, xlo, xhi&
&                      , integral, integrald)
!   call spline ( x_int, y0_int, n, 1.E30, 1.E30, y2 )
!   call splint_integral ( x_int, y0_int, y2, n, xlo, xhi, int0 )
!   call spline ( x_int, y1_int, n, 1.E30, 1.E30, y2 )
!   call splint_integral ( x_int, y1_int, y2, n, xlo, xhi, int1 )
!   call spline ( x_int, y2_int, n, 1.E30, 1.E30, y2 )
!   call splint_integral ( x_int, y2_int, y2, n, xlo, xhi, int2 )
!   call spline ( x_int, y3_int, n, 1.E30, 1.E30, y2 )
!   call splint_integral ( x_int, y3_int, y2, n, xlo, xhi, int3 )
!   call spline ( x_int, yv1_int, n, 1.E30, 1.E30, y2 )
!   call splint_integral ( x_int, yv1_int, y2, n, xlo, xhi, intv1 )
!   call spline ( x_int, yv2_int, n, 1.E30, 1.E30, y2 )
!   call splint_integral ( x_int, yv2_int, y2, n, xlo, xhi, intv2 )
!   
!dF_drho_FMT(i,k) = int0/dhs(k) + 0.5*int1 + PI*dhs(k)*int2 + PI*int3 + intv1/dhs(k) + 2.0*PI*intv2
      df_drho_fmtd(i, k) = parame(k,1)*integrald
      df_drho_fmt(i, k) = integral*parame(k,1)
    END DO
    GOTO 110
 100 stop
 110 CONTINUE
  END SUBROUTINE fmt_dfdrho_d
 
END MODULE mod_dft_fmt_d