Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_fleur_vdW
       8             :   IMPLICIT NONE
       9             :   PRIVATE
      10             :   PUBLIC fleur_vdW_mCallsen
      11             : CONTAINS
      12           0 :   SUBROUTINE fleur_vdW_mCallsen(fmpi,atoms,sphhar,stars,input,      &
      13             :    cell,sym ,juphon,vacuum,results,    &
      14           0 :    den,vpw_total,vr_total)
      15             :     !Interface to Juelich vdW-code
      16             :     USE m_types
      17             :     USE m_constants
      18             :     USE m_psqpw
      19             :     USE m_fft3d
      20             :     USE m_qpwtonmt
      21             :     USE m_convol
      22             :     USE m_cdn_io
      23             : 
      24             :     IMPLICIT NONE
      25             : 
      26             :     TYPE(t_mpi),INTENT(IN)       :: fmpi
      27             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      28             :     TYPE(t_sphhar),INTENT(IN)    :: sphhar
      29             :     TYPE(t_stars),INTENT(IN)     :: stars
      30             :     TYPE(t_vacuum),INTENT(IN)    :: vacuum
      31             :     
      32             :     TYPE(t_cell),INTENT(IN)      :: cell
      33             :     TYPE(t_input),INTENT(IN)     :: input
      34             :     TYPE(t_sym),INTENT(IN)       :: sym
      35             :     TYPE(t_juphon),INTENT(IN)    :: juphon
      36             :     TYPE(t_potden),INTENT(INOUT) :: den
      37             :      
      38             :     TYPE(t_results),INTENT(INOUT) :: results
      39             :     !COMPLEX,INTENT(in)     :: qpw(:)
      40             :     !REAL,INTENT(inout)     :: rho(:,:,:)
      41             :     COMPLEX,INTENT(inout)  :: vpw_total(:)
      42             :     REAL,INTENT(inout)     :: vr_total(:,:,:,:)
      43             : 
      44             : 
      45             :     !locals
      46           0 :     TYPE(t_atoms)      :: atoms_tmp
      47           0 :     REAL,ALLOCATABLE   :: n_grid(:),v_grid(:),rhc(:,:,:)
      48           0 :     COMPLEX,ALLOCATABLE:: vpw(:),psq(:)
      49             :     INTEGER            :: n,ncmsh,j,i
      50             :     LOGICAL            :: l_core
      51           0 :     REAL tec(atoms%ntype,input%jspins),qintc(atoms%ntype,input%jspins)
      52             :     complex            :: sigma_loc(2)
      53             : 
      54           0 :     if (.not.btest(input%vdw,1)) return ! No vdw contribution
      55             : 
      56           0 :     l_core=.FALSE. !try to subtract core charge?
      57           0 :     ALLOCATE(n_grid(27*stars%mx1*stars%mx2*stars%mx3),v_grid(27*stars%mx1*stars%mx2*stars%mx3))
      58           0 :     ALLOCATE(vpw(stars%ng3))
      59           0 :     ALLOCATE(psq(stars%ng3),rhc(atoms%jmtd,atoms%ntype,input%jspins))
      60             : 
      61             :     IF (l_core) l_core = isCoreDensityPresent()
      62             : 
      63             :     IF (l_core.and.btest(input%vdw,3)) THEN
      64             :        WRITE(oUnit,*) "VdW contribution without core charge"       
      65             :        ! read the core charge
      66             :        CALL readCoreDensity(input,atoms,rhc,tec,qintc)
      67             :        DO j=1,input%jspins
      68             :           DO n=1,atoms%ntype
      69             :              ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/atoms%rmsh(1,n) ) / atoms%dx(n) + 1 )
      70             :              ncmsh = MIN( ncmsh, atoms%msh )
      71             :              den%mt(:,1,n,1) = den%mt(:,1,n,1) - rhc(:SIZE(den%mt,1),n,j)/(4. * SQRT( ATAN (1.) ))
      72             :           ENDDO
      73             :        ENDDO
      74             :     ENDIF
      75             : 
      76             :     ! Construct the pseudo charge
      77           0 :     atoms_tmp=atoms
      78           0 :     atoms_tmp%zatom=0.0
      79             :     CALL psqpw(fmpi,&
      80             :          atoms_tmp,sphhar,stars,vacuum,&
      81             :          cell,input,sym , juphon, &
      82           0 :          den,1,.TRUE.,2,psq,sigma_loc)
      83             : 
      84             :     !put pseudo charge on real-space grid
      85             :     !use v_pot for imaginary part
      86             :     CALL fft3d(                        &
      87             :          n_grid,v_grid,psq,         &
      88           0 :          stars,1)
      89             : 
      90             : 
      91             :     CALL priv_fleur_vdW(cell,stars, &
      92           0 :          n_grid,results%e_vdW,v_grid,.TRUE.)
      93             : 
      94           0 :     WRITE(oUnit,*) "------  vdW-Potential code by M. Callsen included-------"
      95           0 :     WRITE(oUnit,*) "vdW-Energy contribution:",results%e_vdW
      96             : 
      97           0 :     if (.not.btest(input%vdw,2)) return ! No vdw contribution to potential
      98             :     
      99             : 
     100             :     !Put potential on rez. grid
     101           0 :     n_grid=0.0
     102             :     CALL fft3d(                           &
     103             :          v_grid,n_grid,vpw,            &
     104           0 :          stars,-1)
     105             : 
     106             :     !Calculate MT-contribution to the potential
     107             : 
     108             :     CALL qpw_to_nmt(                                                     &
     109             :          sphhar,atoms,stars,sym,cell ,fmpi,  &
     110           0 :          1,4,vpw,vr_total)
     111             : 
     112           0 :     WRITE(oUnit,*) "vdW average Potential  :",vpw(1)
     113             : 
     114             : 
     115           0 :     CALL convol(stars,psq,vpw)
     116             : 
     117             :     ! Add to total potential
     118           0 :     vpw_total(:)=vpw_total(:)+psq
     119             : 
     120           0 :   END SUBROUTINE fleur_vdW_mCallsen
     121             : 
     122             : 
     123             : 
     124           0 :   SUBROUTINE priv_fleur_vdW(cell,stars,n_pseudo,e_vdw,v_vdw,l_vdW_v1)
     125             : 
     126             :     USE m_types
     127             :     USE m_constants
     128             :     USE m_juDFT
     129             :     USE param, ONLY:  Zab_v1,Zab_v2
     130             : 
     131             :     USE nonlocal_data, ONLY: nx,ny,nz,                &
     132             :          n_grid,                  &
     133             :          a1,a2,a3,                &
     134             :          b1,b2,b3,                &
     135             :          G_cut,Zab,               &
     136             :          lambda,m_c,omega,tpibya
     137             :     USE nonlocal_funct,ONLY: soler
     138             :     IMPLICIT NONE
     139             :     TYPE(t_cell),INTENT(IN)  :: cell
     140             :     TYPE(t_stars),INTENT(IN) :: stars
     141             :     REAL,INTENT(inout)       :: n_pseudo(:)
     142             :     LOGICAL,INTENT(in)       :: l_vdW_v1
     143             :     REAL,INTENT(out)         :: e_vdw
     144             :     REAL,INTENT(out)         :: v_vdw(:)
     145             : 
     146           0 :     IF (SIZE(n_pseudo).NE.SIZE(v_vdw)) CALL juDFT_error("BUG in fleur_to_vdW")
     147             : 
     148           0 :     a1=cell%amat(:,1)
     149           0 :     a2=cell%amat(:,2)
     150           0 :     a3=cell%amat(:,3)
     151             : 
     152           0 :     b1=cell%bmat(:,1)
     153           0 :     b2=cell%bmat(:,2)
     154           0 :     b3=cell%bmat(:,3)
     155             : 
     156           0 :     nx=3*stars%mx1
     157           0 :     ny=3*stars%mx2
     158           0 :     nz=3*stars%mx3
     159           0 :     n_grid=nx*ny*nz
     160           0 :     IF (SIZE(n_pseudo).NE.n_grid) CALL juDFT_error("BUG2 in fleur_to_vdW")
     161           0 :     g_cut=9*stars%gmax**2  !????
     162             : 
     163           0 :     IF (l_vdW_v1) THEN
     164           0 :        Zab=Zab_v1
     165             :     ELSE
     166           0 :        Zab=Zab_v2
     167             :     ENDIF
     168             : 
     169           0 :     lambda=1.2
     170           0 :     m_c=12
     171           0 :     omega=cell%vol
     172           0 :     tpibya = 2.0*pi_const/omega
     173             : 
     174           0 :     WRITE(oUnit,*)
     175           0 :     WRITE(oUnit,'(A)')      'lattice vectors in Bohr:'
     176           0 :     WRITE(oUnit,'(3F16.8)') a1(:)
     177           0 :     WRITE(oUnit,'(3F16.8)') a2(:)
     178           0 :     WRITE(oUnit,'(3F16.8)') a3(:)
     179           0 :     WRITE(oUnit,*)
     180             : 
     181           0 :     WRITE(oUnit,'(A)')     'reciprocal lattice vectors in Bohr:'
     182           0 :     WRITE(oUnit,'(3F16.8)') b1(:)
     183           0 :     WRITE(oUnit,'(3F16.8)') b2(:)
     184           0 :     WRITE(oUnit,'(3F16.8)') b3(:)
     185           0 :     WRITE(oUnit,'(3F16.8)')
     186             : 
     187           0 :     WRITE(oUnit,'(A,F18.12)') '(2 pi/a) in 1/Bohr^3:',tpibya
     188           0 :     WRITE(oUnit,'(A,F18.12)') 'omega in 1/Bohr^3:   ',omega
     189           0 :     WRITE(oUnit,*)
     190           0 :     WRITE(oUnit,'(A,F18.12)') 'G_cut: ',G_cut
     191             : 
     192           0 :     CALL soler(n_pseudo,e_vdW,v_vdW)
     193             : 
     194           0 :   END SUBROUTINE priv_fleur_vdW
     195             : 
     196             :   SUBROUTINE test_charge(qpw,stars)
     197             :     USE m_types
     198             :     IMPLICIT NONE
     199             :     TYPE(t_stars),INTENT(IN) :: stars
     200             :     COMPLEX,INTENT(out)      :: qpw(:)
     201             :     REAL,SAVE:: pos=0.1
     202             :     COMPLEX,ALLOCATABLE,SAVE::testcharge(:)
     203             :     COMPLEX:: ph
     204             :     INTEGER:: n
     205             :     IF (pos==0.1) THEN
     206             :        ALLOCATE(testcharge(SIZE(qpw)))
     207             :        DO n=1,SIZE(qpw)
     208             :           testcharge(n)=EXP(-.1*DOT_PRODUCT(stars%kv3(:,n),stars%kv3(:,n)))
     209             :        ENDDO
     210             :     ENDIF
     211             :     DO n=1,SIZE(qpw)
     212             :        ph=EXP(CMPLX(0,stars%kv3(3,n)*pos))
     213             :        qpw(n)=(ph+CONJG(ph))*testcharge(n)
     214             :     ENDDO
     215             :     pos=pos+0.01
     216             :     IF (pos>0.4) STOP "testcharge"
     217             : 
     218             :   END SUBROUTINE test_charge
     219             : END MODULE m_fleur_vdW
     220             : 
     221             :