LCOV - code coverage report
Current view: top level - force - force_a12.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 67 73 91.8 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_forcea12
       2             : ! ************************************************************
       3             : ! Pulay 1st term  force contribution a la Rici et al., eq. A12
       4             : !
       5             : ! ************************************************************
       6             : !
       7             : CONTAINS
       8          12 :   SUBROUTINE force_a12(atoms,nobd,sym, DIMENSION, cell,oneD,&
       9          12 :                        we,jsp,ne,usdus,eigVecCoeffs,acoflo,bcoflo,e1cof,e2cof,f_a12,results)
      10             :     USE m_types_setup
      11             :     USE m_types_misc
      12             :     USE m_types_usdus
      13             :     USE m_types_cdnval
      14             :     USE m_constants
      15             :     USE m_juDFT
      16             :     IMPLICIT NONE
      17             : 
      18             :     TYPE(t_results),INTENT(INOUT)   :: results
      19             :     TYPE(t_dimension),INTENT(IN)    :: DIMENSION
      20             :     TYPE(t_oneD),INTENT(IN)         :: oneD
      21             :     TYPE(t_sym),INTENT(IN)          :: sym
      22             :     TYPE(t_cell),INTENT(IN)         :: cell
      23             :     TYPE(t_atoms),INTENT(IN)        :: atoms
      24             :     TYPE(t_usdus),INTENT(IN)        :: usdus
      25             :     TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
      26             :     !     ..
      27             :     !     .. Scalar Arguments ..
      28             :     INTEGER, INTENT (IN) :: nobd    
      29             :     INTEGER, INTENT (IN) :: ne ,jsp 
      30             :     !     ..
      31             :     !     .. Array Arguments ..
      32             :     REAL,    INTENT(IN)    :: we(nobd)
      33             :     COMPLEX, INTENT(IN)    :: acoflo(-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
      34             :     COMPLEX, INTENT(IN)    :: bcoflo(-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
      35             :     COMPLEX, INTENT(IN)    :: e1cof(ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
      36             :     COMPLEX, INTENT(IN)    :: e2cof(ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
      37             :     COMPLEX, INTENT(INOUT) :: f_a12(3,atoms%ntype)
      38             :     !     ..
      39             :     !     .. Local Scalars ..
      40             :     COMPLEX a12,cil1,cil2
      41             :     REAL,PARAMETER:: zero=0.0
      42             :     COMPLEX,PARAMETER:: czero=CMPLX(0.0,0.0)
      43             :     INTEGER i,ie,irinv,is,isinv,it,j,l,l1,l2,lm1,lm2 ,m1,m2,n,natom,natrun,ilo,m
      44             :     !     ..
      45             :     !     .. Local Arrays ..
      46             :     COMPLEX forc_a12(3),gv(3)
      47          24 :     COMPLEX acof_flapw(nobd,0:DIMENSION%lmd),bcof_flapw(nobd,0:DIMENSION%lmd)
      48             :     REAL aaa(2),bbb(2),ccc(2),ddd(2),eee(2),fff(2),gvint(3),starsum(3),vec(3),vecsum(3)
      49             :     !     ..
      50             :     !     .. Statement Functions ..
      51             :     REAL   alpha,beta,delta,epslon,gamma,phi 
      52             :     INTEGER krondel
      53             : 
      54             :     ! Kronecker delta for arguments >=0 AND <0
      55             : 
      56             :     krondel(i,j) = MIN(ABS(i)+1,ABS(j)+1)/MAX(ABS(i)+1,ABS(j)+1)* (1+SIGN(1,i)*SIGN(1,j))/2
      57             :     alpha(l,m) = (l+1)*0.5e0*SQRT(REAL((l-m)* (l-m-1))/ REAL((2*l-1)* (2*l+1)))
      58             :     beta(l,m) = l*0.5e0*SQRT(REAL((l+m+2)* (l+m+1))/ REAL((2*l+1)* (2*l+3)))
      59             :     GAMMA(l,m) = (l+1)*0.5e0*SQRT(REAL((l+m)* (l+m-1))/ REAL((2*l-1)* (2*l+1)))
      60             :     delta(l,m) = l*0.5e0*SQRT(REAL((l-m+2)* (l-m+1))/ REAL((2*l+1)* (2*l+3)))
      61             :     epslon(l,m) = (l+1)*SQRT(REAL((l-m)* (l+m))/ REAL((2*l-1)* (2*l+1)))
      62             :     phi(l,m) = l*SQRT(REAL((l-m+1)* (l+m+1))/REAL((2*l+1)* (2*l+3)))
      63             : 
      64          12 :     CALL timestart("force_a12")
      65             : 
      66          12 :     natom = 1
      67          36 :     DO  n = 1,atoms%ntype
      68          24 :        IF (atoms%l_geo(n)) THEN
      69          24 :           forc_a12(:) = czero
      70             : 
      71             :           !
      72          48 :           DO natrun = natom,natom + atoms%neq(n) - 1
      73             : 
      74          24 :              gv(:) = czero
      75             : 
      76             :              !--->       the local orbitals do not contribute to
      77             :              !--->       the term a12, because they vanish at the
      78             :              !--->       mt-boundary. Therefore, the LO-contribution
      79             :              !--->       to the a and b coefficients has to be
      80             :              !--->       substracted before calculation a12.
      81             :              !
      82         240 :              DO l1 = 0,atoms%lmax(n)
      83        2184 :                 DO m1 = -l1,l1
      84        1944 :                    lm1 = l1* (l1+1) + m1
      85       29376 :                    DO ie = 1,ne
      86       27216 :                       acof_flapw(ie,lm1) = eigVecCoeffs%acof(ie,lm1,natrun,jsp)
      87       29160 :                       bcof_flapw(ie,lm1) = eigVecCoeffs%bcof(ie,lm1,natrun,jsp)
      88             :                    ENDDO
      89             :                 ENDDO
      90             :              ENDDO
      91          24 :              DO ilo = 1,atoms%nlo(n)
      92           0 :                 l1 = atoms%llo(ilo,n)
      93          24 :                 DO m1 = -l1,l1
      94           0 :                    lm1 = l1* (l1+1) + m1
      95           0 :                    DO ie = 1,ne
      96           0 :                       acof_flapw(ie,lm1) = acof_flapw(ie,lm1) - acoflo(m1,ie,ilo,natrun)
      97           0 :                       bcof_flapw(ie,lm1) = bcof_flapw(ie,lm1) - bcoflo(m1,ie,ilo,natrun)
      98             :                    ENDDO
      99             :                 ENDDO
     100             :              ENDDO
     101             :              !
     102         456 :              DO l1 = 0,atoms%lmax(n)
     103         216 :                 cil1 = ImagUnit**l1
     104        2184 :                 DO m1 = -l1,l1
     105        1944 :                    lm1 = l1* (l1+1) + m1
     106       19656 :                    DO l2 = 0,atoms%lmax(n)
     107       17496 :                       cil2 = ImagUnit**l2
     108      176904 :                       DO m2 = -l2,l2
     109      157464 :                          lm2 = l2* (l2+1) + m2
     110             :                          !
     111      157464 :                          a12 = czero
     112     2361960 :                          DO ie = 1,ne
     113             :                             !
     114             :                             a12 = a12 + CONJG(cil1*&
     115             :                                  (acof_flapw(ie,lm1)*usdus%us(l1,n,jsp) + bcof_flapw(ie,lm1)*usdus%uds(l1,n,jsp) ))*cil2*&
     116     2361960 :                                  (e1cof(ie,lm2,natrun)*usdus%us(l2,n,jsp)+ e2cof(ie,lm2,natrun)*usdus%uds(l2,n,jsp))*we(ie)
     117             : 
     118             :                          END DO
     119      157464 :                          aaa(1) = alpha(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1+1)
     120      157464 :                          aaa(2) = alpha(l2,m2)*krondel(l1,l2-1)* krondel(m1,m2+1)
     121      157464 :                          bbb(1) = beta(l1,m1)*krondel(l2,l1+1)* krondel(m2,m1+1)
     122      157464 :                          bbb(2) = beta(l2,m2)*krondel(l1,l2+1)* krondel(m1,m2+1)
     123      157464 :                          ccc(1) = GAMMA(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1-1)
     124      157464 :                          ccc(2) = GAMMA(l2,m2)*krondel(l1,l2-1)* krondel(m1,m2-1)
     125      157464 :                          ddd(1) = delta(l1,m1)*krondel(l2,l1+1)* krondel(m2,m1-1)
     126      157464 :                          ddd(2) = delta(l2,m2)*krondel(l1,l2+1)* krondel(m1,m2-1)
     127      157464 :                          eee(1) = epslon(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1)
     128      157464 :                          eee(2) = epslon(l2,m2)*krondel(l1,l2-1)* krondel(m1,m2)
     129      157464 :                          fff(1) = phi(l1,m1)*krondel(l2,l1+1)* krondel(m2,m1)
     130      157464 :                          fff(2) = phi(l2,m2)*krondel(l1,l2+1)* krondel(m1,m2)
     131             :                          !
     132             :                          gv(1) = gv(1) + (aaa(1)+bbb(1)-ccc(1)-ddd(1)+&
     133      157464 :                               aaa(2)+bbb(2)-ccc(2)-ddd(2))*0.5* atoms%rmt(n)**2*a12
     134             :                          !
     135             :                          gv(2) = gv(2) + ImagUnit* (aaa(1)+bbb(1)+ccc(1)+&
     136      157464 :                               ddd(1)-aaa(2)-bbb(2)-ccc(2)-ddd(2))*0.5* atoms%rmt(n)**2*a12
     137             :                          !
     138      174960 :                          gv(3) = gv(3) + (eee(1)+eee(2)-fff(1)-fff(2))* 0.5*atoms%rmt(n)**2*a12
     139             :                          !
     140             :                          !  m1,m2 loops end
     141             :                       END DO
     142             :                    END DO
     143             :                    !  l1,l2 loops end
     144             :                 END DO
     145             :              END DO
     146             :              !
     147             :              !  to complete summation over stars of k now sum
     148             :              !  over all operations which leave (k+G)*R(natrun)*taual(natrun)
     149             :              !  invariant. We sum over ALL these operations and not only
     150             :              !  the ones needed for the actual star of k. Should be
     151             :              !  ok if we divide properly by the number of operations
     152             :              !  First, we find operation S where RS=T. T -like R- leaves
     153             :              !  the above scalar product invariant (if S=1 then R=T).
     154             :              !  R is the operation which generates position of equivalent atom
     155             :              !  out of position of representative
     156             :              !  S=R^(-1) T
     157             :              !  number of ops which leave (k+G)*op*taual invariant: invarind
     158             :              !  index of inverse operation of R: irinv
     159             :              !  index of operation T: invarop
     160             :              !  now, we calculate index of operation S: is
     161             :              !
     162             :              !  transform vector gv into internal coordinates
     163             : 
     164          24 :              vec(:) = REAL(gv(:)) /atoms%neq(n)
     165             : 
     166          24 :              gvint=MATMUL(cell%bmat,vec)/tpi_const
     167             :              !
     168             : 
     169          24 :              vecsum(:) = zero
     170             : 
     171             :              !-gb2002
     172             :              !            irinv = invtab(ngopr(natrun))
     173             :              !            DO it = 1,invarind(natrun)
     174             :              !               is = multab(irinv,invarop(natrun,it))
     175             :              !c  note, actually we need the inverse of S but -in principle
     176             :              !c  because {S} is agroup and we sum over all S- S should also
     177             :              !c  work; to be lucid we take the inverse:
     178             :              !                isinv = invtab(is)
     179             :              !!               isinv = is
     180             :              ! Rotation is alreadt done in to_pulay, here we work only in the
     181             :              ! coordinate system of the representative atom (natom):
     182             :              !!        
     183         168 :              DO it = 1,sym%invarind(natom)
     184         144 :                 is =sym%invarop(natom,it)
     185         144 :                 isinv = sym%invtab(is)
     186         144 :                 IF (oneD%odi%d1) isinv = oneD%ods%ngopr(natom)
     187             :                 !-gb 2002
     188             :                 !  now we have the wanted index of operation with which we have
     189             :                 !  to rotate gv. Note gv is given in cart. coordinates but
     190             :                 !  mrot acts on internal ones
     191         576 :                 DO i = 1,3
     192         432 :                    vec(i) = zero
     193        1872 :                    DO j = 1,3
     194        1728 :                       IF (.NOT.oneD%odi%d1) THEN
     195        1296 :                          vec(i) = vec(i) + sym%mrot(i,j,isinv)*gvint(j)
     196             :                       ELSE
     197           0 :                          vec(i) = vec(i) + oneD%ods%mrot(i,j,isinv)*gvint(j)
     198             :                       END IF
     199             :                    END DO
     200             :                 END DO
     201        1032 :                 DO i = 1,3
     202         576 :                    vecsum(i) = vecsum(i) + vec(i)
     203             :                 END DO
     204             :                 !   end operator loop
     205             :              END DO
     206             :              !
     207             :              !   transform from internal to cart. coordinates
     208          24 :              starsum=MATMUL(cell%amat,vecsum)
     209         192 :              DO i = 1,3
     210          96 :                 forc_a12(i) = forc_a12(i) + starsum(i)/sym%invarind(natrun)
     211             :              END DO
     212             :              !
     213             :              !  natrun loop end
     214             :           END DO
     215             :           !
     216             :           !
     217             :           !
     218             :           !     sum to existing forces
     219             :           !
     220             :           !  NOTE: force() is real and therefore takes only the
     221             :           !  real part of forc_a12(). in general, force must be
     222             :           !  real after k-star summation. Now, we put the proper
     223             :           !  operations into real space. Problem: what happens
     224             :           !  if in real space there is no inversion any more?
     225             :           !  But we have inversion in k-space due to time reversal
     226             :           !  symmetry, E(k)=E(-k)
     227             :           !  We argue that k-space inversion is automatically taken
     228             :           !  into account if force = (1/2)(forc_a12+conjg(forc_a12))
     229             :           !  because time reversal symmetry means that conjg(PSI)
     230             :           !  is also a solution of Schr. equ. if psi is one.
     231         168 :           DO i = 1,3
     232          72 :              results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a12(i))
     233          96 :              f_a12(i,n) = f_a12(i,n) + forc_a12(i)
     234             :           END DO
     235             :           !
     236             :           !     write result moved to force_a8
     237             :           !
     238             :        ENDIF
     239          36 :        natom = natom + atoms%neq(n)
     240             :     ENDDO
     241             : 
     242          12 :     CALL timestop("force_a12")
     243             : 
     244          12 :   END SUBROUTINE force_a12
     245             : END MODULE m_forcea12

Generated by: LCOV version 1.13