Line data    Source code

       1             :       MODULE DFT_D2
       2             : !
       3             :       IMPLICIT NONE
       4             :       PRIVATE
       5             : !
       6             :       REAL,DIMENSION(94,2) :: c6r
       7             : !
       8             :       PUBLIC :: driver_DFT_D2
       9             : !
      10             :       CONTAINS
      11             : !
      12           0 :       SUBROUTINE initialize_c6r
      13             : !
      14             : ! the first index is "corresponds" to atomic number Z, i.e., 1  to H, 
      15             : !   2 to He etc.
      16             : ! the second index corresponds to 
      17             : !   1 c6 coefficient
      18             : !   2 r_vdW of that element
      19             : !
      20           0 :       c6r(1:94,1)=0.0
      21           0 :       c6r(1:94,1)=1.0
      22             : !
      23           0 :       c6r(1,1) =  1.451020284   ! H atom
      24           0 :       c6r(1,2) =  1.001
      25           0 :       c6r(2,1) =  0.829154448   ! He atom
      26           0 :       c6r(2,2) =  1.012
      27           0 :       c6r(5,1) = 32.440667778   ! B atom
      28           0 :       c6r(5,2) =  1.485
      29           0 :       c6r(6,1) = 18.137753550   ! C atom
      30           0 :       c6r(6,2) =  1.452
      31           0 :       c6r(7,1) = 12.748249638   ! N atom
      32           0 :       c6r(7,2) =  1.397
      33           0 :       c6r(8,1) =  7.25510142            ! O atom
      34           0 :       c6r(8,2) =  1.342
      35           0 :       c6r(9,1) =  7.77332295            ! F atom
      36           0 :       c6r(9,2) =  1.287
      37           0 :       c6r(11,1)= 59.180898726   ! Na atom
      38           0 :       c6r(11,2)=  1.144
      39           0 :       c6r(16,1)= 57.729878442   ! S atom 
      40           0 :       c6r(16,2)=  1.683
      41           0 :       c6r(17,1)= 52.547663142   ! Cl atom
      42           0 :       c6r(17,2)=  1.639
      43           0 :       c6r(19,1)=111.93585048            ! K atom
      44           0 :       c6r(19,2)=  1.485
      45           0 :       c6r(26,1)=111.93585048            ! Fe atom
      46           0 :       c6r(26,2)=  1.562
      47           0 :       c6r(27,1)=111.93585048            ! Co atom
      48           0 :       c6r(27,2)=  1.562
      49           0 :       c6r(29,1)=111.93585048            ! Cu atom
      50           0 :       c6r(29,2)=  1.562
      51           0 :       c6r(35,1)=129.244449582   ! Br atom
      52           0 :       c6r(35,2)=  1.749
      53           0 :       c6r(47,1)=255.690502902   ! Ag atom
      54           0 :       c6r(47,2)=  1.639
      55           0 :       c6r(63,1)=  0.0           ! Eu atom
      56           0 :       c6r(63,2)=  1.0
      57           0 :       c6r(77,1)=130.0           ! Ir atom
      58           0 :       c6r(77,2)=  1.848
      59             : !
      60           0 :       END SUBROUTINE initialize_c6r
      61             : !
      62           0 :       SUBROUTINE calc_ene_DFT_D2(cart_init,z_init,cart_large,z_large)
      63             :       INTEGER, INTENT(IN),DIMENSION(:)   :: z_init,z_large
      64             :       REAL,INTENT(IN),DIMENSION(:,:)     :: cart_init,cart_large
      65             : !
      66             :       INTEGER  :: i1,i2
      67             :       REAL     :: dx,dy,dz
      68             :       REAL     :: c6,c6i,c6j,ri,rj
      69             :       REAL     :: dist,damp,vdW_energy
      70           0 :       REAL,ALLOCATABLE,DIMENSION(:,:) :: force_vdW
      71             : !
      72           0 :       vdW_energy=0.0
      73             : !
      74             : !      print*,'1: ',size(cart_init(3,:)),size(cart_large(3,:))
      75             : !      print*,z_init
      76             : !      print*,'2: ',size(z_init),size(z_large)
      77             : !
      78           0 :       if(size(cart_init(3,:)) <= 0) then
      79           0 :         STOP 'size(cart_init(3,:)) <= 0'
      80             :       else
      81           0 :         allocate(force_vdW(3,size(cart_init(3,:))))
      82           0 :         force_vdW(:,:)=0.0
      83             :       endif
      84             : !
      85           0 :       print'(20x,A)','Forces'
      86             : !
      87           0 :       do i1=1,size(cart_init(3,:))
      88           0 :         c6i=c6r(z_init(i1),1)
      89           0 :         ri =c6r(z_init(i1),2)
      90           0 :         do i2=1,size(cart_large(3,:))
      91           0 :           c6j=c6r(z_large(i2),1)
      92           0 :           rj =c6r(z_large(i2),2)
      93             : !
      94           0 :           c6=0.75*sqrt(c6i*c6j)
      95             : !
      96             :           dist=sqrt( (cart_init(1,i1)-cart_large(1,i2))**2 + 
      97             :      +               (cart_init(2,i1)-cart_large(2,i2))**2 + 
      98           0 :      +               (cart_init(3,i1)-cart_large(3,i2))**2 )
      99             : !
     100             : ! avoid self-interaction
     101             : !
     102           0 :           if(dist > 0.0001) then
     103             : !
     104           0 :             damp = 1.0/(1.0+exp(-20.0*(dist/(ri+rj)-1.0)))
     105             : !
     106           0 :             vdW_energy = vdW_energy - 0.5*damp*c6/dist**6
     107             : !
     108           0 :             dx=cart_init(1,i1)-cart_large(1,i2)
     109           0 :             dy=cart_init(2,i1)-cart_large(2,i2)
     110           0 :             dz=cart_init(3,i1)-cart_large(3,i2)
     111             : !
     112             :             force_vdW(1,i1)=force_vdW(1,i1)+ c6*dx*                    
     113             :      +        ((damp**2)*20.0*exp(-20.0*(dist/(ri+rj)-1.0))/ 
     114           0 :      +        ((ri+rj)*dist**7) - 6.0*damp/dist**8)
     115             :             force_vdW(2,i1)=force_vdW(2,i1)+ c6*dy*
     116             :      +        ((damp**2)*20.0*exp(-20.0*(dist/(ri+rj)-1.0))/ 
     117           0 :      +        ((ri+rj)*dist**7) - 6.0*damp/dist**8)
     118             :             force_vdW(3,i1)=force_vdW(3,i1)+ c6*dz*                    
     119             :      +        ((damp**2)*20.0*exp(-20.0*(dist/(ri+rj)-1.0))/ 
     120           0 :      +        ((ri+rj)*dist**7) - 6.0*damp/dist**8)
     121             : !
     122             :           endif
     123             : !
     124             :         enddo
     125             : !
     126           0 :         print'(A5,I3,3F14.8)','Atom ',i1,force_vdW(1:3,i1)
     127             : !
     128             :       enddo
     129             : !
     130             :       print'(A47,4x,F16.12)',
     131           0 :      &  'vdW energy calculated by DFT-D2 method (in eV):',vdW_energy
     132             : !
     133             :       print'(A47,4x,F16.12)',
     134           0 :      &  'vdW energy calculated by DFT-D2 method (in au):'
     135           0 :      &  ,vdW_energy/27.21138386
     136             : !
     137           0 :       if(allocated(force_vdW)) deallocate(force_vdW)
     138             : !
     139           0 :       END SUBROUTINE calc_ene_DFT_D2
     140             : !
     141           0 :       SUBROUTINE driver_DFT_D2(cart_init,z_init,cart_large,z_large)
     142             :       INTEGER, DIMENSION(:)   :: z_init,z_large
     143             :       REAL    ,DIMENSION(:,:) :: cart_init,cart_large
     144             : !
     145           0 :       call initialize_c6r
     146           0 :       call calc_ene_DFT_D2(cart_init,z_init,cart_large,z_large)
     147             : !
     148           0 :       END SUBROUTINE driver_DFT_D2
     149             : !
     150             :       END MODULE DFT_D2
     151             : !