LCOV - code coverage report
Current view: top level - force - force_sf.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 152 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 4 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_force_sf
       8             : !     *****************************************************************
       9             : !     This routine calculates a contribution to the forces stemming
      10             : !     from the discontinuity of density and potential at the muffin-tin
      11             : !     boundary oint n [ rho V (IS) - rho V (MT) ] dS
      12             : !     Klueppelberg May 13
      13             : !     *****************************************************************
      14             : 
      15             : !     To enable debug code that compares potential and density on the
      16             : !     muffin-tin boundary, uncomment the following line:
      17             : !#define debug
      18             : 
      19             :       IMPLICIT NONE
      20             : 
      21             :       COMPLEX, PRIVATE, SAVE, ALLOCATABLE :: force_mt(:,:)
      22             :       COMPLEX, PRIVATE, SAVE, ALLOCATABLE :: force_is(:,:)
      23             :       LOGICAL, PRIVATE, SAVE :: isdone=.false.,mtdone=.false.
      24             : 
      25             :       CONTAINS
      26             : 
      27           0 :       SUBROUTINE force_sf_is(atoms_in,stars,sym,jsp,cell,qpw,vpw,excpw,vxcpw )
      28             : !     *****************************************************************
      29             : !     This subroutine calculates the contribution evaluated with
      30             : !     quantities from the interstital oint n rho V dS
      31             : !     The Fourier transform of density and potential is first calculated
      32             : !     as an expansion in spherical harmonics. Then the normal vector,
      33             : !     which is proportional to Y_1t connects the l component of rho
      34             : !     with the l+-1 components of V. This is done up to a cutoff lmax.
      35             : !     It is called in a spin loop at the end of vgen.F
      36             : !     *****************************************************************
      37             : 
      38             :       USE m_constants, ONLY : tpi_const
      39             :       USE m_sphbes
      40             :       USE m_phasy1
      41             :       USE m_gaunt
      42             : #ifdef debug
      43             :       USE m_ylm
      44             : #endif
      45             :     USE m_types
      46             :       IMPLICIT NONE
      47             :       TYPE(t_sym),INTENT(IN)     :: sym
      48             :       TYPE(t_stars),INTENT(IN)   :: stars
      49             :       TYPE(t_cell),INTENT(IN)    :: cell
      50             :       TYPE(t_atoms),INTENT(IN)   :: atoms_in
      51           0 :       TYPE(t_atoms)              :: atoms !copy with modified data
      52             : !     .. Scalar Arguments ..
      53             :       INTEGER, INTENT (IN) :: jsp
      54             : 
      55             : !     .. Array Arguments ..
      56             :       COMPLEX, INTENT (IN) :: qpw(:,:) !(stars%ng3,input%jspins)
      57             :       COMPLEX, INTENT (IN) :: vpw(:,:)!(n3d,jspd)
      58             :       COMPLEX, INTENT (IN) :: excpw(stars%ng3)
      59             :       COMPLEX, INTENT (IN) :: vxcpw(:,:)!(stars%ng3,input%jspins)
      60             : 
      61             : !     .. Local Scalars ..
      62             :       INTEGER :: n,j,itype,s,l ,lm,t,lp,mp,lmp,jp,natom,m
      63             :       REAL    :: r,r2
      64             :       COMPLEX :: img,rhoprep,Vprep
      65             :       LOGICAL :: isthere
      66             : 
      67             : !     .. Local Arrays ..
      68           0 :       INTEGER :: lmaxb(atoms_in%ntype)
      69           0 :       COMPLEX :: coeff(3,-1:1),qpw2(stars%ng3,size(qpw,2)),qpwcalc(stars%ng3,size(qpw,2))
      70           0 :       REAL   , ALLOCATABLE :: bsl(:,:,:)
      71           0 :       COMPLEX, ALLOCATABLE :: pylm(:,:,:),rho(:),V(:),pylm2(:,:)
      72             : !       COMPLEX, ALLOCATABLE :: qpw2(:,:),qpwcalc(:,:)
      73             : #ifdef debug
      74             :       REAL    :: vec(3)
      75             :       COMPLEX :: factorrho,factorv
      76             :       COMPLEX, ALLOCATABLE :: ylm(:),testrho(:,:),testV(:,:)
      77             : #endif
      78           0 :       atoms=atoms_in
      79           0 :       atoms%lmax = 2*atoms_in%lmaxd!60!
      80           0 :       lmaxb = atoms%lmax
      81           0 :       img = cmplx(0.0,1.0)
      82             : 
      83           0 :       CALL init_sf(sym,cell,atoms)
      84           0 :       isdone = .true.
      85             : 
      86           0 :       ALLOCATE ( bsl(stars%ng3,0:atoms%lmaxd,atoms%ntype) )
      87             : 
      88           0 :       ALLOCATE ( pylm2((atoms%lmaxd+1)**2,atoms%ntype ))
      89           0 :       ALLOCATE ( rho((atoms%lmaxd+1)**2),V((atoms%lmaxd+1)**2) )
      90             : 
      91             : #ifdef debug
      92             :       ALLOCATE ( ylm((atoms%lmaxd+1)**2),testrho((atoms%lmaxd+1)**2,atoms%ntype) )
      93             :       ALLOCATE ( testV((atoms%lmaxd+1)**2,atoms%ntype) )
      94             : #endif
      95             : 
      96           0 :       coeff(:, :) =   cmplx(0.0,0.0)
      97           0 :       coeff(1,-1) =     sqrt(tpi_const/3.)
      98           0 :       coeff(1, 1) =    -sqrt(tpi_const/3.)
      99           0 :       coeff(2,-1) = img*sqrt(tpi_const/3.)
     100           0 :       coeff(2, 1) = img*sqrt(tpi_const/3.)
     101           0 :       coeff(3, 0) =  sqrt(2.*tpi_const/3.)
     102           0 :       WRITE (1704,*) 'excpw:',excpw
     103           0 :       WRITE (1704,*) 'vxcpw:',vxcpw
     104             : !     load in density without coretails
     105           0 :       INQUIRE (FILE='qpw',EXIST=isthere)
     106             :       IF (isthere.and..false.) THEN
     107             :         qpw2 = 0.0
     108             :         OPEN (15,file='qpw',form='formatted',status='unknown')
     109             :         DO jp = 1,size(qpw,2)
     110             :         DO s = 1,stars%ng3
     111             :           READ (15,'(i2,i10,2f20.14)') n,j,qpw2(s,jp)
     112             :         END DO ! s
     113             :         END DO ! jp
     114             :         CLOSE (15)
     115             :         IF (any(abs(qpw2).gt.10**(-6))) THEN
     116             :           qpwcalc = qpw2
     117             :         ELSE
     118             :           qpwcalc = qpw
     119             :         END IF
     120             :       ELSE
     121           0 :         qpwcalc = qpw
     122             :       END IF
     123             : 
     124             : !     prepare star quantities
     125             : !       DO s = 1,ng3
     126             : !         DO itype = 1,ntype
     127             : !           r = sk3(s)*rmt(itype)
     128             : !           CALL sphbes(lmax,r,bsl(s,:,itype))
     129             : !         END DO ! itype
     130             : !       END DO ! s
     131           0 :       DO itype = 1,atoms%ntype
     132           0 :         CALL sphbes(atoms%lmax(itype),0.0,bsl(1,:,itype))
     133           0 :         DO s = 2,stars%ng3
     134             : !         Only call sphbes if the length of the star changed
     135           0 :           IF (abs(stars%sk3(s)-stars%sk3(s-1)).gt.1.0e-14) THEN
     136           0 :             r = stars%sk3(s)*atoms%rmt(itype)
     137           0 :             CALL sphbes(atoms%lmax(itype),r,bsl(s,:,itype))
     138             :           ELSE
     139           0 :             bsl(s,:,itype) = bsl(s-1,:,itype)
     140             :           END IF
     141             :         END DO ! s
     142             :       END DO ! itype
     143             : 
     144             : 
     145           0 :       force_is = 0.0
     146             : 
     147           0 :       natom = 1
     148           0 :       DO itype = 1,atoms%ntype
     149           0 :         r2  = atoms%rmt(itype)**2
     150           0 :         rho = 0.0
     151           0 :         V   = 0.0
     152             : !         DO l = 0,lmax-1
     153             : !           DO s = 1,ng3
     154             : ! !           calculate phase factors for the current atom type to prevent
     155             : ! !           total overhead. it is still calculated lmax times even 
     156             : ! !           though phasy1 generates pylm for all values of l at once.
     157             : ! !           but this loop sequence is more convenient.
     158             : ! !           allocating pylm as pylm(lm,s,itype) and precalculating it
     159             : ! !           leads to exhaustive use of memory in larger systems
     160             : ! !             CALL phasy1(
     161             : ! ! !      >              ntypd,n3d,natd,nop,lmaxd,ntype,neq,lmax,
     162             : ! !      >              1,n3d,1,nop,lmax,1,neq(itype),lmaxb(itype),0,
     163             : ! !      >              2.*tpi,taual(1:3,natom),bmat,kv3,!1+sum(neq(1:itype-1))),bmat,kv3,
     164             : ! !      >              tau,mrot,symor,s,invtab,
     165             : ! !      <              pylm2(:))
     166             : !             pylm2 = 1.
     167             : ! !             IF ((s.eq.1).and.(l.eq.0)) WRITE (851,*) pylm2(:)
     168             : !             rhoprep = nstr(s) * bsl(s,l,itype) * qpwcalc(s,jsp)
     169             : !             IF (l.eq.0) THEN
     170             : !               Vprep = nstr(s) * bsl(s,l,itype)
     171             : !      *              * (1*vpw(s,jsp)-1*vxcpw(s,jsp)+1*excpw(s))
     172             : ! !      *              * vpw(s,jsp)
     173             : !             END IF
     174             : ! !           for l = 0 we calculate rho_00 and V_00...
     175             : !             DO m = -l,l
     176             : !               lm = l*(l+1) + m + 1
     177             : !               rho(lm) = rho(lm) + rhoprep * pylm2(lm)!pylm(lm,s,itype)!
     178             : !               IF (l.gt.0) CYCLE
     179             : !                 V(lm) =   V(lm) +   Vprep * pylm2(lm)!pylm(lm,s,itype)!
     180             : !             END DO ! m
     181             : ! !           ... and V_1mp, for l > 0, we calculate rho_lm, V_l+1,mp ...
     182             : !               Vprep = nstr(s) * bsl(s,l+1,itype)
     183             : !      *              * (1*vpw(s,jsp)-1*vxcpw(s,jsp)+1*excpw(s))
     184             : ! !      *              * vpw(s,jsp)
     185             : !             DO m = -l-1,l+1
     186             : !               lm = (l+1)*(l+2) + m + 1
     187             : !                 V(lm) =   V(lm) +   Vprep * pylm2(lm)!pylm(lm,s,itype)!
     188             : !             END DO ! m
     189             : !           END DO ! s
     190             : 
     191           0 :         DO s = 1,stars%ng3 !l = 0,atoms%lmax-1
     192             : !         calculate phase factors for the current atom type to prevent
     193             : !         total overhead. it is still calculated lmax times even 
     194             : !         though phasy1 generates pylm for all values of l at once.
     195             : !         but this loop sequence is more convenient.
     196             : !         allocating pylm as pylm(lm,s,itype) and precalculating it
     197             : !         leads to exhaustive use of memory in larger systems
     198           0 :           CALL phasy1(atoms,stars,sym,cell,s,pylm2(:,:))
     199             : 
     200           0 :           DO l = 0,atoms%lmax(itype)-1 !s = 1,stars%ng3
     201             : !           calculate phase factors for the current atom type to prevent
     202             : !           total overhead. it is still calculated lmax times even 
     203             : !           though phasy1 generates pylm for all values of l at once.
     204             : !           but this loop sequence is more convenient.
     205             : !           allocating pylm as pylm(lm,s,itype) and precalculating it
     206             : !           leads to exhaustive use of memory in larger systems
     207             : !             CALL phasy1(
     208             : ! !      >              ntypd,n3d,natd,nop,lmaxd,ntype,neq,lmax,
     209             : !      >              1,n3d,1,nop,lmax,1,neq(itype),lmaxb(itype),0,
     210             : !      >              2.*tpi,taual(1:3,natom),bmat,kv3,!1+sum(neq(1:itype-1))),bmat,kv3,
     211             : !      >              tau,mrot,symor,s,invtab,
     212             : !      <              pylm2(:))
     213             : !             IF ((s.eq.1).and.(l.eq.0)) WRITE (851,*) pylm2(:)
     214           0 :             rhoprep = stars%nstr(s) * bsl(s,l,itype) * qpwcalc(s,jsp)
     215           0 :             IF (l.eq.0) THEN
     216           0 :               Vprep = stars%nstr(s) * bsl(s,l,itype) * (1*vpw(s,jsp)-1*vxcpw(s,jsp)+1*excpw(s)) ! Switching between Veff and VCoul + exc
     217             : !      *              * vpw(s,jsp)
     218             :             END IF
     219             : !           for l = 0 we calculate rho_00 and V_00...
     220           0 :             DO m = -l,l
     221           0 :               lm = l*(l+1) + m + 1
     222           0 :               rho(lm) = rho(lm) + rhoprep * pylm2(lm,itype)!pylm(lm,s,itype)!
     223           0 :               IF (l.gt.0) CYCLE
     224           0 :                 V(lm) =   V(lm) +   Vprep * pylm2(lm,itype)!pylm(lm,s,itype)!
     225             :             END DO ! m
     226             : !           ... and V_1mp, for l > 0, we calculate rho_lm, V_l+1,mp ...
     227           0 :               Vprep = stars%nstr(s) * bsl(s,l+1,itype) * (1*vpw(s,jsp)-1*vxcpw(s,jsp)+1*excpw(s))
     228             : !      *              * vpw(s,jsp)
     229           0 :             DO m = -l-1,l+1
     230           0 :               lm = (l+1)*(l+2) + m + 1
     231           0 :                 V(lm) =   V(lm) +   Vprep * pylm2(lm,itype)!pylm(lm,s,itype)!
     232             :             END DO ! m
     233             :           END DO ! l
     234             :           END DO ! s
     235             : 
     236             : !           V = 0.0
     237             : !           V(1) = sqrt(2.*tpi)
     238             : !           rho = 0.0
     239             : !           rho(1) = sqrt(2.*tpi)
     240             : 
     241           0 :         DO l = 0,atoms%lmax(itype)-1 ! new: altered s and l loop above
     242             : 
     243           0 :           DO m = -l,l
     244           0 :             lm = l*(l+1) + m + 1
     245           0 :             WRITE (1705,*) itype,l,m,lm-1,rho(lm),V(lm)
     246             : !           ... because rho_lm occurs with V_l-1,mp and V_l+1,mp
     247           0 :             DO lp = abs(l-1),l+1,2
     248           0 :               DO t = -1,1
     249           0 :                 mp = t-m
     250           0 :                 IF (lp.lt.abs(mp)) CYCLE
     251           0 :                 lmp = lp*(lp+1) + mp + 1
     252             :                 force_is(:,itype) = force_is(:,itype) + r2&
     253           0 :                      * rho(lm) * V(lmp) * conjg(coeff(:,t)) * gaunt1(1,l,lp,t,m,mp,atoms%lmax(itype))
     254             :               END DO ! t
     255             :             END DO ! lp
     256             :           END DO ! m
     257             :         END DO ! l
     258             : #ifdef debug
     259             :         testrho(:,itype) = rho(:)
     260             :           testV(:,itype) =   V(:)
     261             : #endif
     262           0 :         WRITE (849,'(3(2(f20.14),1x))') (force_is(s,itype),s=1,3)
     263           0 :         natom = natom + atoms%neq(itype)
     264             :       END DO ! itype
     265             : 
     266             : #ifdef debug
     267             : !     test output
     268             : !     Evaluate plane wave density in spherical harmonic representation
     269             : !     on MT surface along a vector (1,1,1)
     270             :       DO itype = 1,atoms%ntype
     271             :         WRITE (180,'(2i3)') jsp,itype
     272             :         WRITE (180,'(3f20.14)') real(force_is(1,itype)), real(force_is(2,itype)), real(force_is(3,itype))
     273             :         vec(1) =-0.730! 1.0
     274             :         vec(2) = 0.938! 1.0
     275             :         vec(3) =-1.625! 1.0
     276             :         CALL ylm4(atoms%lmax,vec,ylm)
     277             :         factorrho = 0.0
     278             :         factorv   = 0.0
     279             :         DO l = 0,atoms%lmax
     280             :         DO m = -l,l
     281             :           lm = l*(l+1) + m + 1
     282             :           factorrho = factorrho + testrho(lm,itype) * ylm(lm)
     283             :           factorv   = factorv   +   testV(lm,itype) * ylm(lm)
     284             :           WRITE (180,'(3i4,2f20.14)') l,m,lm,testV(lm,itype)
     285             :         END DO ! jp
     286             :         END DO ! lapw%kp
     287             :         WRITE (180,'(a5,2f20.14)') 'den: ',factorrho
     288             :         WRITE (180,'(a5,2f20.14)') 'pot: ',factorv
     289             :       END DO ! itype
     290             :       DEALLOCATE ( ylm,testrho,testV )
     291             : #endif
     292             : 
     293           0 :       DEALLOCATE ( bsl,rho,V )
     294             : !       DEALLOCATE ( pylm )
     295           0 :       DEALLOCATE ( pylm2 )
     296             : !       DEALLOCATE ( qpw2 )
     297             : !       DEALLOCATE ( qpwcalc )
     298             : 
     299           0 :       END SUBROUTINE force_sf_is
     300             : 
     301             : 
     302             : 
     303           0 :       SUBROUTINE force_sf_mt(&
     304             :                             atoms,sphhar,jspin,&
     305             :                             ispin,mpi,&
     306           0 :                             vr,excr,&
     307           0 :                             vxcr,rho,&
     308             :                             sym,cell )
     309             : !     *****************************************************************
     310             : !     This subroutine calculates the contribution evaluated with
     311             : !     quantities from the muffin tin
     312             : !     n rho V = sum(nu,nup) rho(nu)V(nup)
     313             : !             * sum(m,mu,mup) c_1m* Y_1m* c_lnu Y_numu c_lnup Y_nupmup
     314             : !     It is called in a spin loop at the end of cdnval.F
     315             : !     *****************************************************************
     316             : 
     317             :       USE m_constants, ONLY : tpi_const,sfp_const
     318             :       USE m_gaunt
     319             :       USE m_ylm
     320             :       USE m_types
     321             :       IMPLICIT NONE
     322             :       TYPE(t_mpi),INTENT(IN)   :: mpi
     323             :       TYPE(t_sym),INTENT(IN)   :: sym
     324             :       TYPE(t_cell),INTENT(IN)  :: cell
     325             :       TYPE(t_sphhar),INTENT(IN):: sphhar
     326             :       TYPE(t_atoms),INTENT(IN) :: atoms
     327             : 
     328             : !     .. Scalar Arguments ..
     329             :       INTEGER, INTENT (IN) :: jspin
     330             :       INTEGER, INTENT (IN) :: ispin 
     331             : 
     332             : !     .. Array Arguments ..
     333             :       REAL   , INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) ! 
     334             :       REAL   , INTENT (IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
     335             :       REAL   , INTENT (IN) :: excr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
     336             :       REAL   , INTENT (IN) :: vxcr(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
     337             : 
     338             : !     .. Local Scalars ..
     339             :       INTEGER :: natom,itype,nd,lh,l,lhp,lp,mem,m,memp,mp,t,i,lmp
     340             :       REAL    :: pot,den
     341             :       COMPLEX :: img,factor
     342             : 
     343             : !     .. Local Arrays ..
     344             :       COMPLEX :: coeff(3,-1:1)
     345           0 :       COMPLEX :: d1((atoms%lmaxd+1)**2,atoms%ntype ),d2((atoms%lmaxd+1)**2,atoms%ntype )
     346             : #ifdef debug
     347             :       COMPLEX :: testrho((atoms%lmaxd+1)**2 ,atoms%ntype),testv((atoms%lmaxd+1)**2,atoms%ntype )
     348             :       COMPLEX :: ylm((atoms%lmaxd+1)**2),factorrho,factorv
     349             :       REAL    :: r,Gext(3)
     350             :       INTEGER :: lm
     351             :       testrho = 0.0
     352             :       testv   = 0.0
     353             : 
     354             :       WRITE (181,'(a2)') 'MT'
     355             : #endif
     356           0 :       WRITE (1704,*) 'vxcr:',vxcr
     357           0 :       WRITE (1704,*) 'excr:',excr
     358             : 
     359           0 :       CALL init_sf(sym,cell,atoms)
     360           0 :       mtdone = .true.
     361             : 
     362           0 :       img = cmplx(0.0,1.0)
     363           0 :       force_mt = 0.0
     364             : 
     365           0 :       coeff(:, :) =   cmplx(0.0,0.0)
     366           0 :       coeff(1,-1) =     sqrt(tpi_const/3.)
     367           0 :       coeff(1, 1) =    -sqrt(tpi_const/3.)
     368           0 :       coeff(2,-1) = img*sqrt(tpi_const/3.)
     369           0 :       coeff(2, 1) = img*sqrt(tpi_const/3.)
     370           0 :       coeff(3, 0) =  sqrt(2.*tpi_const/3.)
     371             : 
     372           0 :       d1 = 0
     373           0 :       d2 = 0
     374             : 
     375             : !     Calculate forces: For each atom, loop over all lattice harmonics.
     376           0 :       natom = 1
     377           0 :       DO itype = 1,atoms%ntype
     378             : 
     379           0 :         nd = atoms%ntypsy(natom)
     380             : 
     381           0 :         DO lh = 0,sphhar%nlh(nd)
     382           0 :           l = sphhar%llh(lh,nd)
     383             : 
     384             : !         The l=0 component of the potential array is saved with a
     385             : !         factor r/sfp in front of it. For this calculation, we need
     386             : !         the pure potential
     387           0 :           pot = vr(atoms%jri(itype),lh,itype)
     388           0 :           IF (lh.eq.0) THEN
     389           0 :             pot = pot*sfp_const/atoms%rmt(itype)
     390             :           END IF
     391           0 :           pot = +1*excr(atoms%jri(itype),lh,itype) -1*vxcr(atoms%jri(itype),lh,itype,ispin) +1*pot
     392             : !           pot = 0
     393             : !           IF (l.eq.0) pot = sfp
     394             : 
     395           0 :           WRITE (400,'(3(i4,1x),f20.14)') itype,lh,l,vr(atoms%jri(itype),lh,itype)
     396           0 :           WRITE (401,'(3(i4,1x),f20.14)') itype,lh,l,rho(atoms%jri(itype),lh,itype,ispin)
     397           0 :           DO mem = 1,sphhar%nmem(lh,nd)
     398           0 :             m = sphhar%mlh(mem,lh,nd)
     399           0 :             lmp = l*(l+1) + m + 1
     400             :             d1(lmp,itype) = d1(lmp,itype) + sphhar%clnu(mem,lh,nd)&
     401           0 :                  *rho(atoms%jri(itype),lh,itype,ispin)/atoms%rmt(itype)**2
     402           0 :             d2(lmp,itype) = d2(lmp,itype) + sphhar%clnu(mem,lh,nd) * pot
     403             :           END DO ! mem
     404             : 
     405           0 :           DO lhp = 0,sphhar%nlh(nd)
     406           0 :             lp = sphhar%llh(lhp,nd)
     407           0 :             IF (abs(l-lp).ne.1) CYCLE
     408             : 
     409           0 :             den = rho(atoms%jri(itype),lhp,itype,ispin) ! jspin ?
     410             : !             den = 0
     411             : !             IF (lp.eq.0) den = sfp*rmt(itype)**2
     412             : 
     413           0 :             DO mem = 1,sphhar%nmem(lh,nd)
     414           0 :               m = sphhar%mlh(mem,lh,nd)
     415           0 :             DO memp = 1,sphhar%nmem(lhp,nd)
     416           0 :               mp = sphhar%mlh(memp,lhp,nd)
     417           0 :               IF (abs(m+mp).gt.1) CYCLE
     418             : 
     419             : !             Due to the normal vector n, the lattice harmonics meet
     420             : !             with a Y_1m resulting in a gaunt coefficient.
     421             :               factor = pot * den * sphhar%clnu(mem,lh,nd) * sphhar%clnu(memp,lhp,nd)&
     422           0 :                     * gaunt1(1,l,lp,m+mp,m,mp,atoms%lmaxd)
     423           0 :               force_mt(:,itype) = force_mt(:,itype) + factor * conjg(coeff(:,m+mp))
     424             : 
     425             :             END DO ! memp
     426             :             END DO ! mem
     427             : 
     428             :           END DO ! lhp
     429             : 
     430             : 
     431             : #ifdef debug
     432             : !         testrho/v debug code
     433             : !         construct spherical harmonic expansion of lattice harmonic
     434             : !         density and potential on the muffin-tin boundary
     435             :           DO mem = 1,sphhar%nmem(lh,nd)
     436             :             m = sphhar%mlh(mem,lh,nd)
     437             :           DO lp = 0,atoms%lmaxd
     438             :             IF (l.ne.lp) CYCLE
     439             :           DO mp = -lp,lp
     440             :             IF (m.ne.mp) CYCLE
     441             :             lm = lp*(lp+1) + mp + 1
     442             :             testrho(lm,itype) = testrho(lm,itype)&
     443             :                  + rho(atoms%jri(itype),lh,itype,ispin) &
     444             :                  * sphhar%clnu(mem,lh,nd) / atoms%rmt(itype)**2
     445             :             testv(lm,itype)   = testv(lm,itype) + pot * sphhar%clnu(mem,lh,nd)
     446             :           END DO ! mp
     447             :           END DO ! lp
     448             :           END DO ! mem
     449             : #endif
     450             : 
     451             :         END DO ! lh
     452             : 
     453           0 :         DO l = 0,atoms%lmaxd
     454           0 :         DO m = -l,l
     455           0 :           lmp = l*(l+1) + m + 1
     456           0 :           WRITE (1706,*) itype,l,m,lmp-1,d1(lmp,itype),d2(lmp,itype)
     457             :         END DO ! m
     458             :         END DO ! l
     459             : 
     460           0 :         IF (mpi%irank.eq.0) THEN
     461           0 :           WRITE (850,'(3(2(f20.14),1x))') (force_mt(t,itype),t=1,3)
     462             :         END IF
     463           0 :         natom = natom + atoms%neq(itype)
     464             :       END DO ! itype
     465             : 
     466             : !     debug
     467           0 :       WRITE (*,*) 'sanity check'
     468           0 :       WRITE (*,*) 'look is',force_is
     469           0 :       WRITE (*,*) 'look mt',force_mt
     470             : 
     471             : #ifdef debug
     472             : !     test output
     473             : !     Evaluate lattice harmonic density in spherical harmonic
     474             : !     representation in along a reciprocal vector (1,1,1)
     475             :       DO itype = 1,atoms%ntype
     476             :         WRITE (181,'(2i3)') jspin,itype
     477             :         WRITE (181,'(3f20.14)') real(force_mt(1,itype)), real(force_mt(2,itype)), real(force_mt(3,itype))
     478             :         Gext(1) =-0.730! 1.0
     479             :         Gext(2) = 0.938! 1.0
     480             :         Gext(3) =-1.625! 1.0
     481             :         CALL ylm4(atoms%lmaxd,Gext,ylm)
     482             :         factorrho = 0.0
     483             :         factorv   = 0.0
     484             :         DO lp = 0,atoms%lmaxd
     485             :         DO mp = -lp,lp
     486             :           lm = lp*(lp+1) + mp + 1
     487             :           factorrho = factorrho + testrho(lm,itype) * ylm(lm)
     488             :           factorv   = factorv   + testv(lm,itype)   * ylm(lm)
     489             :           WRITE (181,'(3i4,2f20.14)') lp,mp,lm,testv(lm,itype)
     490             :         END DO ! mp
     491             :         END DO ! lp
     492             :         WRITE (181,'(a5,2f20.14)') 'den: ',factorrho
     493             :         WRITE (181,'(a5,2f20.14)') 'pot: ',factorv
     494             :       END DO ! itype
     495             : #endif
     496             : 
     497           0 :       END SUBROUTINE force_sf_mt
     498             : 
     499           0 :       SUBROUTINE init_sf(sym,cell,atoms)
     500             : !     Initialize results arrays if neither force_sf_is nor force_sf_mt
     501             : !     were executed up until now.
     502             : !     Called at the beginning of cdnval.F to once fill the wigner array
     503             : !     Also called in force_sf_is/mt to guarantee that
     504             : !     all arrays are allocated.
     505             : 
     506             :     USE m_types
     507             :       IMPLICIT NONE
     508             :       TYPE(t_sym),INTENT(IN)   :: sym
     509             :       TYPE(t_cell),INTENT(IN)  :: cell
     510             :       TYPE(t_atoms),INTENT(IN) :: atoms
     511             : 
     512             : !     debug
     513           0 :       IF (ALLOCATED(force_is)) THEN
     514           0 :         WRITE (*,*) 'init is:',force_is
     515           0 :         WRITE (*,*) 'init mt:',force_mt
     516             :       ELSE
     517           0 :         WRITE (*,*) 'init: not initialized'
     518             :       END IF
     519             : 
     520           0 :       IF (isdone.OR.mtdone) RETURN
     521           0 :       IF (.not.ALLOCATED(force_is)) THEN
     522           0 :         ALLOCATE ( force_is(3,atoms%ntype),force_mt(3,atoms%ntype) )
     523             :       END IF
     524           0 :       force_is = 0.0
     525           0 :       force_mt = 0.0
     526             : 
     527           0 :         WRITE (*,*) 'init: end'
     528             :       END SUBROUTINE init_sf
     529             : 
     530           0 :       SUBROUTINE exit_sf(isp,atoms,force)
     531             : !     Write out force contribution from surface and deallocate arrays if
     532             : !     all force_sf_is and force_sf_mt were executed
     533             : !     Called at the end of cdnval.F and totale.f for writing purposes
     534             : !     and to deallocate the arrays.
     535             : 
     536             :     USE m_types
     537             :       IMPLICIT NONE
     538             :       INTEGER,INTENT(IN)         :: isp
     539             :       TYPE(t_atoms),INTENT(IN)   :: atoms
     540             :       REAL,INTENT(INOUT)         :: force(:,:,:)
     541             : 
     542             :       INTEGER :: itype,dir
     543           0 :       COMPLEX :: force_sf(3,atoms%ntype)
     544             : 
     545             : !     debug
     546           0 :       IF (ALLOCATED(force_is)) THEN
     547           0 :         WRITE (*,*) 'exit is:',force_is
     548           0 :         WRITE (*,*) 'exit mt:',force_mt
     549             :       ELSE
     550           0 :         WRITE (*,*) 'exit: not initialized'
     551             :       END IF
     552             : 
     553           0 :       IF (isdone.AND.mtdone) THEN
     554           0 :         force_sf(:,:) = force_is(:,:) - force_mt(:,:)
     555           0 :         force(:,:,isp) = force(:,:,isp) + real(force_sf(:,:))
     556           0 :         WRITE (6,*)
     557           0 :         DO itype = 1,atoms%ntype
     558           0 :           WRITE (6,FMT=8010) itype
     559           0 :           WRITE (6,FMT=8020) (force_sf(dir,itype),dir=1,3)
     560             :         END DO ! itype
     561           0 :         isdone = .false.
     562           0 :         mtdone = .false.
     563           0 :         DEALLOCATE ( force_is,force_mt )
     564             :       END IF
     565             : 
     566             : 8010  FORMAT (' FORCES: SURFACE CORRECTION FOR ATOM TYPE',i4)
     567             : 8020  FORMAT (' FX_SF=',2f10.6,' FY_SF=',2f10.6,' FZ_SF=',2f10.6)
     568             : 
     569           0 :       END SUBROUTINE exit_sf
     570             : 
     571             :       END MODULE m_force_sf

Generated by: LCOV version 1.13