LCOV - code coverage report
Current view: top level - force - force_sf.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 125 126 99.2 %
Date: 2024-03-29 04:21:46 Functions: 4 4 100.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_force_sf
       8             :    USE m_constants
       9             :    USE m_types
      10             :    !-----------------------------------------------------------------------------
      11             :    ! This routine calculates a contribution to the forces stemming
      12             :    ! from the discontinuity of density and potential at the muffin-tin
      13             :    ! boundary oint n [ rho V (IS) - rho V (MT) ] dS
      14             :    ! Klueppelberg May 13
      15             :    !
      16             :    ! The surface integral is found in the first two lines of equation (14) in 
      17             :    ! 
      18             :    ! Klüppelberg et al., PRB 91 035105 (2015).
      19             :    ! 
      20             :    ! Klueppelberg (force level 3) 
      21             :    !-----------------------------------------------------------------------------
      22             : 
      23             :    IMPLICIT NONE
      24             : 
      25             :    COMPLEX, PRIVATE, SAVE, ALLOCATABLE :: force_mt(:,:)
      26             :    COMPLEX, PRIVATE, SAVE, ALLOCATABLE :: force_is(:,:)
      27             :    LOGICAL, PRIVATE, SAVE :: isdone=.false.,mtdone=.false.
      28             : 
      29             : CONTAINS
      30           1 :    SUBROUTINE force_sf_is(atoms_in,stars,sym,jsp,cell,qpw,vpw,excpw,vxcpw)
      31             :       !--------------------------------------------------------------------------
      32             :       ! This subroutine calculates oint n rho V dS evaluated with quantities 
      33             :       ! from the interstital. The Fourier transform of density and potential is 
      34             :       ! first calculated as an expansion in spherical harmonics. Then the normal 
      35             :       ! vector, which is proportional to Y_1t, connects the l component of rho
      36             :       ! with the l+-1 components of V. This is done up to a cutoff lmax.
      37             :       ! It is called in a spin loop at the end of vgen.F90-
      38             :       !--------------------------------------------------------------------------
      39             : 
      40             :       USE m_sphbes
      41             :       USE m_phasy1
      42             :       USE m_gaunt
      43             : 
      44             :       TYPE(t_sym),INTENT(IN)     :: sym
      45             :       TYPE(t_stars),INTENT(IN)   :: stars
      46             :       TYPE(t_cell),INTENT(IN)    :: cell
      47             :       TYPE(t_atoms),INTENT(IN)   :: atoms_in
      48             : 
      49             :       INTEGER, INTENT (IN) :: jsp
      50             : 
      51             :       COMPLEX, INTENT (IN) :: qpw(:,:) !(stars%ng3,input%jspins)
      52             :       COMPLEX, INTENT (IN) :: vpw(:,:) !(n3d,jspd)
      53             :       COMPLEX, INTENT (IN) :: excpw(stars%ng3)
      54             :       COMPLEX, INTENT (IN) :: vxcpw(:,:) !(stars%ng3,input%jspins)
      55             : 
      56           1 :       TYPE(t_atoms)              :: atoms
      57             : 
      58             :       INTEGER :: n,j,itype,s,l ,lm,t,lp,mp,lmp,jp,natom,m
      59             :       REAL    :: r,r2
      60             :       COMPLEX :: img,rhoprep,Vprep
      61             :       LOGICAL :: isthere
      62             : 
      63           1 :       INTEGER :: lmaxb(atoms_in%ntype)
      64           1 :       COMPLEX :: coeff(3,-1:1),qpw2(stars%ng3,size(qpw,2)),qpwcalc(stars%ng3,size(qpw,2))
      65           1 :       REAL   , ALLOCATABLE :: bsl(:,:,:)
      66           1 :       COMPLEX, ALLOCATABLE :: pylm(:,:,:),rho(:),V(:),pylm2(:,:)
      67             : 
      68           1 :       CALL timestart("Force level 3 (IS)")
      69             : 
      70           4 :       atoms=atoms_in
      71           5 :       atoms%lmax = 2*atoms_in%lmax
      72           1 :       atoms%lmaxd = 2*atoms_in%lmaxd
      73           4 :       lmaxb = atoms%lmax
      74           1 :       img = cmplx(0.0,1.0)
      75             : 
      76           1 :       CALL init_sf(sym,cell,atoms)
      77           1 :       isdone = .true.
      78             : 
      79           5 :       ALLOCATE ( bsl(0:atoms%lmaxd,stars%ng3,atoms%ntype) )
      80             : 
      81           4 :       ALLOCATE ( pylm2((atoms%lmaxd+1)**2,atoms%ntype ))
      82           4 :       ALLOCATE ( rho((atoms%lmaxd+1)**2),V((atoms%lmaxd+1)**2) )
      83             : 
      84           1 :       coeff(:, :) =   cmplx(0.0,0.0)
      85           1 :       coeff(1,-1) =     sqrt(tpi_const/3.)
      86           1 :       coeff(1, 1) =    -sqrt(tpi_const/3.)
      87           1 :       coeff(2,-1) = img*sqrt(tpi_const/3.)
      88           1 :       coeff(2, 1) = img*sqrt(tpi_const/3.)
      89           1 :       coeff(3, 0) =  sqrt(2.*tpi_const/3.)
      90             : 
      91        3556 :       qpwcalc = qpw
      92             : 
      93           4 :       DO itype = 1,atoms%ntype
      94           3 :          CALL sphbes(atoms%lmax(itype),0.0,bsl(:,1,itype))
      95       10663 :          DO s = 2,stars%ng3
      96             :             ! Only call sphbes if the length of the star changed.
      97       10662 :             IF (abs(stars%sk3(s)-stars%sk3(s-1)).gt.1.0e-14) THEN
      98       10659 :                r = stars%sk3(s)*atoms%rmt(itype)
      99       10659 :                CALL sphbes(atoms%lmax(itype),r,bsl(:,s,itype))
     100             :             ELSE
     101           0 :                bsl(:,s,itype) = bsl(:,s-1,itype)
     102             :             END IF
     103             :          END DO ! s
     104             :       END DO ! itype
     105             : 
     106             : 
     107          13 :       force_is = 0.0
     108             : 
     109           4 :       DO itype = 1,atoms%ntype
     110           3 :          natom = atoms%firstAtom(itype)
     111           3 :          r2  = atoms%rmt(itype)**2
     112         870 :          rho = 0.0
     113         870 :          V   = 0.0
     114             : 
     115       10665 :          DO s = 1,stars%ng3
     116       10662 :             CALL phasy1(atoms,stars,sym,cell,s,pylm2(:,:))
     117             : 
     118      181257 :             DO l = 0,atoms%lmax(itype)-1
     119      170592 :                rhoprep = stars%nstr(s) * bsl(l,s,itype) * qpwcalc(s,jsp)
     120      170592 :                IF (l.eq.0) THEN
     121             :                   Vprep = stars%nstr(s) * bsl(l,s,itype) * &
     122       10662 :                                           (vpw(s,jsp)-vxcpw(s,jsp)+excpw(s)) 
     123             :                                           ! Switching between Veff and VCoul + exc
     124             :                END IF
     125             :                ! For l = 0 we calculate rho_00 and V_00:
     126     2900064 :                DO m = -l,l
     127     2729472 :                   lm = l*(l+1) + m + 1
     128     2729472 :                   rho(lm) = rho(lm) + rhoprep * pylm2(lm,itype)!pylm(lm,s,itype)!
     129     2729472 :                   IF (l.gt.0) CYCLE
     130     2900064 :                   V(lm) =   V(lm) +   Vprep * pylm2(lm,itype)!pylm(lm,s,itype)!
     131             :                END DO ! m
     132             :                ! And V_1mp, for l > 0, we calculate rho_lm, V_l+1,mp:
     133             :                Vprep = stars%nstr(s) * bsl(l+1,s,itype) * &
     134      170592 :                                        (vpw(s,jsp)-vxcpw(s,jsp)+excpw(s))
     135     3251910 :                DO m = -l-1,l+1
     136     3070656 :                   lm = (l+1)*(l+2) + m + 1
     137     3241248 :                   V(lm) =   V(lm) +   Vprep * pylm2(lm,itype)!pylm(lm,s,itype)!
     138             :                END DO ! m
     139             :             END DO ! l
     140             :          END DO ! s
     141             : 
     142          52 :          DO l = 0,atoms%lmax(itype)-1 ! new: altered s and l loop above
     143         819 :             DO m = -l,l
     144         768 :                lm = l*(l+1) + m + 1
     145             :                ! Because rho_lm occurs with V_l-1,mp and V_l+1,mp:
     146        2349 :                DO lp = abs(l-1),l+1,2
     147        6900 :                   DO t = -1,1
     148        4599 :                      mp = t-m
     149        4599 :                      IF (lp.lt.abs(mp)) CYCLE
     150        4329 :                      lmp = lp*(lp+1) + mp + 1
     151             :                      force_is(:,itype) = force_is(:,itype) + r2 * rho(lm) * &
     152             :                                          V(lmp) * conjg(coeff(:,t)) * &
     153       18849 :                                          gaunt2(1,l,lp,t,m,mp,atoms%lmax(itype))
     154             :                   END DO ! t
     155             :                END DO ! lp
     156             :             END DO ! m
     157             :          END DO ! l
     158             :       END DO ! itype
     159             : 
     160           1 :       DEALLOCATE ( bsl,rho,V )
     161           1 :       DEALLOCATE ( pylm2 )
     162             : 
     163           1 :       CALL timestop("Force level 3 (IS)")
     164             : 
     165           4 :    END SUBROUTINE force_sf_is
     166             : 
     167           1 :    SUBROUTINE force_sf_mt(atoms,sphhar,jspin,ispin,fmpi,vr,excr,vxcr,rho,sym,cell )
     168             :       !--------------------------------------------------------------------------
     169             :       ! This subroutine calculates the contribution evaluated with quantities
     170             :       ! from the Muffin Tins.
     171             :       ! 
     172             :       ! n rho V = sum(nu,nup) rho(nu)V(nup)
     173             :       !           * sum(m,mu,mup) c_1m* Y_1m* c_lnu Y_numu c_lnup Y_nupmup
     174             :       !
     175             :       ! It is called in a spin loop at the end of cdngen.F90
     176             :       !--------------------------------------------------------------------------
     177             : 
     178             :       USE m_gaunt
     179             :       USE m_ylm
     180             : 
     181             :       TYPE(t_mpi),INTENT(IN)   :: fmpi
     182             :       TYPE(t_sym),INTENT(IN)   :: sym
     183             :       TYPE(t_cell),INTENT(IN)  :: cell
     184             :       TYPE(t_sphhar),INTENT(IN):: sphhar
     185             :       TYPE(t_atoms),INTENT(IN) :: atoms
     186             : 
     187             :       INTEGER, INTENT (IN) :: jspin
     188             :       INTEGER, INTENT (IN) :: ispin 
     189             : 
     190             :       REAL   , INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) ! 
     191             :       REAL   , INTENT (IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
     192             :       REAL   , INTENT (IN) :: excr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
     193             :       REAL   , INTENT (IN) :: vxcr(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
     194             : 
     195             :       INTEGER :: natom,itype,nd,lh,l,lhp,lp,mem,m,memp,mp,t,i,lmp
     196             :       REAL    :: pot,den
     197             :       COMPLEX :: img,factor
     198             : 
     199             :       COMPLEX :: coeff(3,-1:1)
     200           1 :       COMPLEX :: d1((atoms%lmaxd+1)**2,atoms%ntype ),d2((atoms%lmaxd+1)**2,atoms%ntype )
     201             : 
     202           1 :       CALL timestart("Force level 3 (MT)")
     203             : 
     204           1 :       CALL init_sf(sym,cell,atoms)
     205           1 :       mtdone = .true.
     206             : 
     207           1 :       img = cmplx(0.0,1.0)
     208          13 :       force_mt = 0.0
     209             : 
     210           1 :       coeff(:, :) =   cmplx(0.0,0.0)
     211           1 :       coeff(1,-1) =     sqrt(tpi_const/3.)
     212           1 :       coeff(1, 1) =    -sqrt(tpi_const/3.)
     213           1 :       coeff(2,-1) = img*sqrt(tpi_const/3.)
     214           1 :       coeff(2, 1) = img*sqrt(tpi_const/3.)
     215           1 :       coeff(3, 0) =  sqrt(2.*tpi_const/3.)
     216             : 
     217         247 :       d1 = 0
     218         247 :       d2 = 0
     219             : 
     220             :       ! Calculate forces: For each atom, loop over all lattice harmonics.
     221           4 :       DO itype = 1,atoms%ntype
     222           3 :          natom = atoms%firstAtom(itype)
     223           3 :          nd = sym%ntypsy(natom)
     224             : 
     225         247 :          DO lh = 0,sphhar%nlh(nd)
     226         243 :             l = sphhar%llh(lh,nd)
     227             : 
     228             :             ! The l=0 component of the potential array is saved with an additional 
     229             :             ! factor r/sfp. For this calculation, we need the pure potential.
     230         243 :             pot = vr(atoms%jri(itype),lh,itype)
     231         243 :             IF (lh.eq.0) THEN
     232           3 :                pot = pot*sfp_const/atoms%rmt(itype)
     233             :             END IF
     234             : 
     235         243 :             pot = excr(atoms%jri(itype),lh,itype)-vxcr(atoms%jri(itype),lh,itype,ispin)+pot
     236             : 
     237         702 :             DO mem = 1,sphhar%nmem(lh,nd)
     238         459 :                m = sphhar%mlh(mem,lh,nd)
     239         459 :                lmp = l*(l+1) + m + 1
     240             :                d1(lmp,itype) = d1(lmp,itype) + sphhar%clnu(mem,lh,nd) * &
     241         459 :                                rho(atoms%jri(itype),lh,itype,ispin)/atoms%rmt(itype)**2
     242         702 :                d2(lmp,itype) = d2(lmp,itype) + sphhar%clnu(mem,lh,nd) * pot
     243             :             END DO ! mem
     244             : 
     245       19929 :             DO lhp = 0,sphhar%nlh(nd)
     246       19683 :                lp = sphhar%llh(lhp,nd)
     247       19683 :                IF (abs(l-lp).ne.1) CYCLE
     248             : 
     249        4848 :                den = rho(atoms%jri(itype),lhp,itype,ispin)
     250             : 
     251       14355 :                DO mem = 1,sphhar%nmem(lh,nd)
     252        9264 :                   m = sphhar%mlh(mem,lh,nd)
     253       46659 :                   DO memp = 1,sphhar%nmem(lhp,nd)
     254       17712 :                      mp = sphhar%mlh(memp,lhp,nd)
     255       17712 :                      IF (abs(m+mp).gt.1) CYCLE
     256             : 
     257             :                      ! Due to the normal vector n, the lattice harmonics form a
     258             :                      ! Gaunt coefficient with the Y_1m.
     259             :                      factor = pot * den * sphhar%clnu(mem,lh,nd) * sphhar%clnu(memp,lhp,nd)&
     260        4104 :                                   * gaunt1(1,l,lp,m+mp,m,mp,atoms%lmaxd)
     261             :                
     262       25680 :                      force_mt(:,itype) = force_mt(:,itype) + factor * conjg(coeff(:,m+mp))
     263             : 
     264             :                   END DO ! memp
     265             :                END DO ! mem
     266             : 
     267             :             END DO ! lhp
     268             : 
     269             :          END DO ! lh
     270             :       END DO ! itype
     271             : 
     272           1 :       CALL timestop("Force level 3 (MT)")
     273             : 
     274           1 :    END SUBROUTINE force_sf_mt
     275             : 
     276           4 :    SUBROUTINE init_sf(sym,cell,atoms)
     277             :       ! Initialize force arrays if neither force_sf_is nor force_sf_mt were
     278             :       ! executed yet.
     279             :       ! Called at the beginning of cdnval.F90 to once fill the wigner array
     280             :       ! Also called in force_sf_is/mt to guarantee that all arrays are allocated.
     281             : 
     282             :       TYPE(t_sym),INTENT(IN)   :: sym
     283             :       TYPE(t_cell),INTENT(IN)  :: cell
     284             :       TYPE(t_atoms),INTENT(IN) :: atoms
     285             : 
     286           4 :       IF (isdone.OR.mtdone) RETURN
     287           2 :       IF (.NOT.ALLOCATED(force_is)) THEN
     288           8 :          ALLOCATE ( force_is(3,atoms%ntype),force_mt(3,atoms%ntype) )
     289             :       END IF
     290          26 :       force_is = 0.0
     291          26 :       force_mt = 0.0
     292             : 
     293             :    END SUBROUTINE init_sf
     294             : 
     295           2 :    SUBROUTINE exit_sf(isp,atoms,force)
     296             :       ! Write out the force contribution from the surface terms and deallocate
     297             :       ! arrays if all force_sf_is and force_sf_mt were executed.
     298             :       ! Called at the end of totale.f90.
     299             : 
     300             :       INTEGER,INTENT(IN)         :: isp
     301             :       TYPE(t_atoms),INTENT(IN)   :: atoms
     302             :       REAL,INTENT(INOUT)         :: force(:,:,:)
     303             : 
     304             :       INTEGER :: itype,dir
     305           2 :       COMPLEX :: force_sf(3,atoms%ntype)
     306             : 
     307           2 :       IF (isdone.AND.mtdone) THEN
     308          13 :          force_sf(:,:) = force_is(:,:) - force_mt(:,:)
     309          13 :          force(:,:,isp) = force(:,:,isp) + real(force_sf(:,:))
     310           1 :          WRITE (oUnit,*)
     311           4 :          DO itype = 1,atoms%ntype
     312           3 :             WRITE (oUnit,FMT=8010) itype
     313           4 :             WRITE (oUnit,FMT=8020) (force_sf(dir,itype),dir=1,3)
     314             :          END DO ! itype
     315           1 :          isdone = .false.
     316           1 :          mtdone = .false.
     317           1 :          DEALLOCATE ( force_is,force_mt )
     318             :       END IF
     319             : 
     320             : 8010  FORMAT (' FORCES: SURFACE CORRECTION FOR ATOM TYPE',i4)
     321             : 8020  FORMAT (' FX_SF=',2f10.6,' FY_SF=',2f10.6,' FZ_SF=',2f10.6)
     322             : 
     323           2 :    END SUBROUTINE exit_sf
     324             : 
     325             : END MODULE m_force_sf

Generated by: LCOV version 1.14