LCOV - code coverage report
Current view: top level - vgen - fleur_vdW.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 76 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 3 0.0 %

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

Generated by: LCOV version 1.13