LCOV - code coverage report
Current view: top level - force - force_a21.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 101 118 85.6 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_forcea21
       2             : CONTAINS
       3          12 :   SUBROUTINE force_a21(input,atoms,DIMENSION,sym,oneD,cell,&
       4          12 :                        we,jsp,epar,ne,eig,usdus,eigVecCoeffs,aveccof,bveccof,cveccof,f_a21,f_b4,results)
       5             : 
       6             :     ! ************************************************************
       7             :     ! Pulay 2nd and 3rd (A17+A20) term force contribution a la Rici
       8             :     ! combined
       9             :     ! NOTE: we do NOT include anymore  the i**l factors
      10             :     ! in the alm,blm coming from to_pulay. Therefore, we can
      11             :     ! use matrixelements from file 28,38 DIRECTLY
      12             :     ! note: present version only yields forces for
      13             :     ! highest energy window (=valence states)
      14             :     ! if also semicore forces are wanted the tmas and tmat files
      15             :     ! have to be saved, indexed and properly used here in force_a21
      16             :     ! 22/june/97: probably we found symmetrization error replacing
      17             :     ! now S^-1 by S (IS instead of isinv)
      18             :     ! ************************************************************
      19             :     !
      20             :     ! Force contribution B4 added following
      21             :     ! Madsen, Blaha, Schwarz, Sjostedt, Nordstrom
      22             :     ! GMadsen FZJ 20/3-01
      23             :     !
      24             :     USE m_forcea21lo
      25             :     USE m_forcea21U
      26             :     USE m_tlmplm_store
      27             :     USE m_types_setup
      28             :     USE m_types_misc
      29             :     USE m_types_usdus
      30             :     USE m_types_tlmplm
      31             :     USE m_types_cdnval
      32             :     USE m_constants
      33             :     USE m_juDFT
      34             :     IMPLICIT NONE
      35             :     TYPE(t_input),INTENT(IN)        :: input
      36             :     TYPE(t_results),INTENT(INOUT)   :: results
      37             :     TYPE(t_dimension),INTENT(IN)    :: DIMENSION
      38             :     TYPE(t_oneD),INTENT(IN)         :: oneD
      39             :     TYPE(t_sym),INTENT(IN)          :: sym
      40             :     TYPE(t_cell),INTENT(IN)         :: cell
      41             :     TYPE(t_atoms),INTENT(IN)        :: atoms
      42             :     TYPE(t_usdus),INTENT(IN)        :: usdus
      43             :     TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
      44             :     !     ..
      45             :     !     .. Scalar Arguments ..
      46             :     INTEGER, INTENT (IN) :: ne,jsp
      47             :     !     ..
      48             :     !     .. Array Arguments ..
      49             :     REAL,    INTENT(IN)    :: we(ne),epar(0:atoms%lmaxd,atoms%ntype)
      50             :     REAL,    INTENT(IN)    :: eig(DIMENSION%neigd)
      51             :     COMPLEX, INTENT(IN)    :: aveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
      52             :     COMPLEX, INTENT(IN)    :: bveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
      53             :     COMPLEX, INTENT(IN)    :: cveccof(3,-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
      54             :     COMPLEX, INTENT(INOUT) :: f_a21(3,atoms%ntype)
      55             :     COMPLEX, INTENT(INOUT) :: f_b4(3,atoms%ntype)
      56             :     !     ..
      57             :     !     .. Local Scalars ..
      58             :     COMPLEX dtd,dtu,utd,utu
      59             :     INTEGER lo, mlotot, mlolotot, mlot_d, mlolot_d
      60             :     INTEGER i,ie,im,in,l1,l2,ll1,ll2,lm1,lm2,m1,m2,n,natom,m,i_u
      61             :     INTEGER natrun,is,isinv,j,irinv,it
      62             :     REAL   ,PARAMETER:: zero=0.0
      63             :     COMPLEX,PARAMETER:: czero=CMPLX(0.,0.)
      64             :     !     ..
      65             :     !     .. Local Arrays ..
      66          12 :     COMPLEX, ALLOCATABLE :: v_mmp(:,:,:)
      67          12 :     REAL,    ALLOCATABLE :: a21(:,:),b4(:,:)
      68             :     COMPLEX forc_a21(3),forc_b4(3)
      69             :     REAL starsum(3),starsum2(3),gvint(3),gvint2(3)
      70             :     REAL vec(3),vec2(3),vecsum(3),vecsum2(3)
      71             : 
      72          12 :     TYPE(t_tlmplm)::tlmplm
      73             : 
      74          12 :     CALL timestart("force_a21")
      75             : 
      76             :     !dimension%lmplmd = (dimension%lmd* (dimension%lmd+3))/2
      77          12 :     mlotot = 0 ; mlolotot = 0
      78          36 :     DO n = 1, atoms%ntype
      79          24 :        mlotot = mlotot + atoms%nlo(n)
      80          36 :        mlolotot = mlolotot + atoms%nlo(n)*(atoms%nlo(n)+1)/2
      81             :     ENDDO
      82          12 :     mlot_d = MAX(mlotot,1)
      83          12 :     mlolot_d = MAX(mlolotot,1)
      84             :     ALLOCATE ( tlmplm%tdd(0:DIMENSION%lmplmd,atoms%ntype,1),tlmplm%tuu(0:DIMENSION%lmplmd,atoms%ntype,1),&
      85             :          tlmplm%tdu(0:DIMENSION%lmplmd,atoms%ntype,1),tlmplm%tud(0:DIMENSION%lmplmd,atoms%ntype,1),&
      86             :          tlmplm%tuulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,1),&
      87             :          tlmplm%tdulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,1),&
      88             :          tlmplm%tuloulo(-atoms%llod:atoms%llod,-atoms%llod:atoms%llod,mlolot_d,1),&
      89          12 :          a21(3,atoms%nat),b4(3,atoms%nat),tlmplm%ind(0:DIMENSION%lmd,0:DIMENSION%lmd,atoms%ntype,1) )
      90             :     !
      91          12 :     IF(atoms%n_u.GT.0) THEN
      92          12 :        ALLOCATE(v_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u))
      93          60 :        v_mmp = CMPLX(0.0,0.0)
      94          12 :        CALL read_tlmplm_vs_mmp(jsp,atoms%n_u,v_mmp)
      95             :     END IF
      96             : 
      97          12 :     i_u = 1
      98          12 :     natom = 1
      99          36 :     DO  n = 1,atoms%ntype
     100          24 :        IF (atoms%l_geo(n)) THEN
     101          24 :           forc_a21(:) = czero
     102          24 :           forc_b4(:) = czero
     103             : 
     104             : 
     105             :           CALL read_tlmplm(n,jsp,atoms%nlo,&
     106             :                tlmplm%tuu(:,n,1),tlmplm%tud(:,n,1),tlmplm%tdu(:,n,1),tlmplm%tdd(:,n,1),&
     107          24 :                tlmplm%ind(:,:,n,1),tlmplm%tuulo(:,:,:,1),tlmplm%tuloulo(:,:,:,1),tlmplm%tdulo(:,:,:,1))
     108             : 
     109          48 :           DO natrun = natom,natom + atoms%neq(n) - 1
     110          24 :              a21(:,natrun) = zero
     111         120 :              b4(:,natrun) = zero
     112             : 
     113             :           END DO
     114             :           !
     115         360 :           DO ie = 1,ne
     116             :              !
     117             :              !
     118         360 :              DO l1 = 0,atoms%lmax(n)
     119        3024 :                 ll1 = l1* (l1+1)
     120       30240 :                 DO m1 = -l1,l1
     121       27216 :                    lm1 = ll1 + m1
     122       27216 :                    DO l2 = 0,atoms%lmax(n)
     123             :                       !
     124      244944 :                       ll2 = l2* (l2+1)
     125     2449440 :                       DO m2 = -l2,l2
     126     2204496 :                          lm2 = ll2 + m2
     127     4653936 :                          DO natrun = natom,natom + atoms%neq(n) - 1
     128     2204496 :                             in = tlmplm%ind(lm1,lm2,n,1)
     129     4408992 :                             IF (in.NE.-9999) THEN
     130     1925616 :                                IF (in.GE.0) THEN
     131             :                                   !
     132             :                                   ! ATTENTION: the matrix elements tuu,tdu,tud,tdd
     133             :                                   ! as calculated in tlmplm are the COMPLEX CONJUGATE
     134             :                                   ! of the non-spherical matrix elements because in the
     135             :                                   ! matrix building routine hssphn (or similar routines)
     136             :                                   ! the COMPLEX CONJUGATE of alm,blm is calculated (to
     137             :                                   ! save complex operations presumably)
     138             :                                   ! Her, A20 is formulated in the usual way therefore
     139             :                                   ! we have to take the COMPLEX CONJUGATE versions
     140             :                                   ! of tuu,tdu,tud,tdd as compared to hssphn!
     141             :                                   !
     142      975744 :                                   utu = tlmplm%tuu(in,n,1)
     143      975744 :                                   dtu = tlmplm%tdu(in,n,1)
     144      975744 :                                   utd = tlmplm%tud(in,n,1)
     145      975744 :                                   dtd = tlmplm%tdd(in,n,1)
     146             :                                ELSE
     147      949872 :                                   im = -in
     148      949872 :                                   utu = CONJG(tlmplm%tuu(im,n,1))
     149      949872 :                                   dtd = CONJG(tlmplm%tdd(im,n,1))
     150      949872 :                                   utd = CONJG(tlmplm%tdu(im,n,1))
     151      949872 :                                   dtu = CONJG(tlmplm%tud(im,n,1))
     152             :                                END IF
     153     7702464 :                                DO i = 1,3
     154             :                                   a21(i,natrun) = a21(i,natrun) + 2.0*&
     155             :                                        AIMAG( CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utu*aveccof(i,ie,lm2,natrun)&
     156             :                                        +CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utd*bveccof(i,ie,lm2,natrun)&
     157             :                                        +CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtu*aveccof(i,ie,lm2,natrun)&
     158     7702464 :                                        +CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtd*bveccof(i,ie,lm2,natrun))*we(ie)/atoms%neq(n)
     159             :                                   !   END i loop
     160             :                                END DO
     161             :                             END IF
     162             :                             !   END natrun
     163             :                          END DO
     164             :                          !
     165             :                          !   END m2 loop
     166             :                       END DO
     167             :                       !   END l2 loop
     168             :                    END DO
     169             :                    !+gu 20.11.97
     170       27216 :                    in = tlmplm%ind(lm1,lm1,n,1)
     171       27216 :                    IF (in.NE.-9999) THEN
     172       25872 :                       utu = -eig(ie)
     173       25872 :                       utd = 0.0
     174       25872 :                       dtu = 0.0
     175       25872 :                       dtd = utu*usdus%ddn(l1,n,jsp)
     176             :                    ELSE
     177        1344 :                       utu = epar(l1,n)-eig(ie)
     178        1344 :                       utd = 0.5
     179        1344 :                       dtu = 0.5
     180        1344 :                       dtd = utu*usdus%ddn(l1,n,jsp)
     181             :                    END IF
     182      111888 :                    DO i = 1,3
     183      272160 :                       DO natrun = natom,natom + atoms%neq(n) - 1
     184             :                          a21(i,natrun) = a21(i,natrun) + 2.0*&
     185             :                               AIMAG(CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utu*aveccof(i,ie,lm1,natrun)&
     186             :                               +CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utd*bveccof(i,ie,lm1,natrun)&
     187             :                               +CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtu*aveccof(i,ie,lm1,natrun)&
     188             :                               +CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtd*bveccof(i,ie,lm1,natrun)&
     189      163296 :                               )*we(ie) /atoms%neq(n)
     190             :                       END DO
     191             :                       !
     192             :                       !-gu
     193             :                       ! END  i loop
     194             :                    END DO
     195             :                    !   END m1 loop
     196             :                 END DO
     197             :                 !   END l1 loop
     198             :              END DO
     199             :              !   END ie loop
     200             :           END DO
     201             :           !
     202             :           !--->    add the local orbital and U contribution to a21
     203             :           !
     204          24 :           CALL force_a21_lo(atoms,jsp,n,we,eig,ne,eigVecCoeffs,aveccof,bveccof,cveccof,tlmplm,usdus,a21)
     205             : 
     206          24 :           IF ((atoms%n_u.GT.0).AND.(i_u.LE.atoms%n_u)) THEN
     207          24 :              CALL force_a21_U(atoms,i_u,n,jsp,we,ne,usdus,v_mmp,eigVecCoeffs,aveccof,bveccof,cveccof,a21)
     208             :           END IF
     209          24 :           IF (input%l_useapw) THEN
     210             :              ! -> B4 force
     211           0 :              DO ie = 1,ne
     212           0 :                 DO l1 = 0,atoms%lmax(n)
     213           0 :                    ll1 = l1* (l1+1)
     214           0 :                    DO m1 = -l1,l1
     215           0 :                       lm1 = ll1 + m1
     216           0 :                       DO i = 1,3
     217           0 :                          DO natrun = natom,natom + atoms%neq(n) - 1
     218             :                             b4(i,natrun) = b4(i,natrun) + 0.5 *&
     219             :                                  we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
     220             :                                  CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)*usdus%us(l1,n,jsp)&
     221             :                                  +eigVecCoeffs%bcof(ie,lm1,natrun,jsp)*usdus%uds(l1,n,jsp))*&
     222             :                                  (aveccof(i,ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
     223             :                                  +bveccof(i,ie,lm1,natrun)*usdus%duds(l1,n,jsp) )&
     224             :                                  -CONJG(aveccof(i,ie,lm1,natrun)*usdus%us(l1,n,jsp)&
     225             :                                  +bveccof(i,ie,lm1,natrun)*usdus%uds(l1,n,jsp) )*&
     226             :                                  (eigVecCoeffs%acof(ie,lm1,natrun,jsp)*usdus%dus(l1,n,jsp)&
     227           0 :                                  +eigVecCoeffs%bcof(ie,lm1,natrun,jsp)*usdus%duds(l1,n,jsp)) )
     228             :                          END DO
     229             :                       END DO
     230             :                    END DO
     231             :                 END DO
     232           0 :                 DO lo = 1,atoms%nlo(n)
     233           0 :                    l1 = atoms%llo(lo,n)
     234           0 :                    DO m = -l1,l1
     235           0 :                       lm1 = l1* (l1+1) + m
     236           0 :                       DO i=1,3
     237           0 :                          DO natrun = natom,natom + atoms%neq(n) - 1
     238             :                             b4(i,natrun) = b4(i,natrun) + 0.5 *&
     239             :                                  we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
     240             :                                  CONJG( eigVecCoeffs%acof(ie,lm1,natrun,jsp)* usdus%us(l1,n,jsp)&
     241             :                                  + eigVecCoeffs%bcof(ie,lm1,natrun,jsp)* usdus%uds(l1,n,jsp) ) *&
     242             :                                  cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp)&
     243             :                                  + CONJG(eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)*usdus%ulos(lo,n,jsp)) *&
     244             :                                  ( aveccof(i,ie,lm1,natrun)* usdus%dus(l1,n,jsp)&
     245             :                                  + bveccof(i,ie,lm1,natrun)* usdus%duds(l1,n,jsp)&
     246             :                                  + cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) )  &
     247             :                                  - (CONJG( aveccof(i,ie,lm1,natrun) *usdus%us(l1,n,jsp)&
     248             :                                  + bveccof(i,ie,lm1,natrun) *usdus%uds(l1,n,jsp) ) *&
     249             :                                  eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)  *usdus%dulos(lo,n,jsp)&
     250             :                                  + CONJG(cveccof(i,m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
     251             :                                  ( eigVecCoeffs%acof(ie,lm1,natrun,jsp)*usdus%dus(l1,n,jsp)&
     252             :                                  + eigVecCoeffs%bcof(ie,lm1,natrun,jsp)*usdus%duds(l1,n,jsp)&
     253           0 :                                  + eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)*usdus%dulos(lo,n,jsp) ) ) )  
     254             :                          END DO
     255             :                       ENDDO
     256             :                    ENDDO
     257             :                 ENDDO
     258             :              END DO
     259             :           ENDIF
     260             :           !
     261          48 :           DO natrun = natom,natom + atoms%neq(n) - 1
     262             :              !
     263             :              !  to complete summation over stars of k now sum
     264             :              !  over all operations which leave (k+G)*R(natrun)*taual(natrun)
     265             :              !  invariant. We sum over ALL these operations and not only
     266             :              !  the ones needed for the actual star of k. Should be
     267             :              !  ok if we divide properly by the number of operations
     268             :              !  First, we find operation S where RS=T. T -like R- leaves
     269             :              !  the above scalar product invariant (if S=1 then R=T).
     270             :              !  R is the operation which generates position of equivalent atom
     271             :              !  out of position of representative
     272             :              !  S=R^(-1) T
     273             :              !  number of ops which leave (k+G)*op*taual invariant: invarind
     274             :              !  index of inverse operation of R: irinv
     275             :              !  index of operation T: invarop
     276             :              !  now, we calculate index of operation S: is
     277             :              !
     278             :              !  note, that vector in expression A17,A20 + A21 is a
     279             :              !  reciprocal lattice vector! other transformation rules
     280             :              !
     281             :              !  transform recip vector g-g' into internal coordinates
     282             : 
     283          24 :              vec(:) = a21(:,natrun)
     284          96 :              vec2(:) = b4(:,natrun)
     285             : 
     286          24 :              gvint=MATMUL(cell%bmat,vec)/tpi_const
     287          24 :              gvint2=MATMUL(cell%bmat,vec2)/tpi_const
     288          24 :              vecsum(:) = zero
     289          24 :              vecsum2(:) = zero
     290             : 
     291             :              !-gb2002
     292             :              !            irinv = invtab(ngopr(natrun))
     293             :              !            DO it = 1,invarind(natrun)
     294             :              !               is = multab(irinv,invarop(natrun,it))
     295             :              !c  note, actually we need the inverse of S but -in principle
     296             :              !c  because {S} is a group and we sum over all S- S should also
     297             :              !c  work; to be lucid we take the inverse:
     298             :              !                isinv = invtab(is)
     299             :              !!               isinv = is
     300             :              ! Rotation is alreadt done in to_pulay, here we work only in the
     301             :              ! coordinate system of the representative atom (natom):
     302             :              !!        
     303         168 :              DO it = 1,sym%invarind(natom)
     304         144 :                 is =sym%invarop(natom,it)
     305         144 :                 isinv = sym%invtab(is)
     306         144 :                 IF (oneD%odi%d1) isinv = oneD%ods%ngopr(natom)
     307             :                 !-gb 2002
     308             :                 !  now we have the wanted index of operation with which we have
     309             :                 !  to rotate gv. Note gv is given in cart. coordinates but
     310             :                 !  mrot acts on internal ones
     311         576 :                 DO i = 1,3
     312         432 :                    vec(i) = zero
     313         432 :                    vec2(i) = zero
     314        1872 :                    DO j = 1,3
     315        1728 :                       IF (.NOT.oneD%odi%d1) THEN
     316        1296 :                          vec(i) = vec(i) + sym%mrot(i,j,isinv)*gvint(j)
     317        1296 :                          vec2(i) = vec2(i) + sym%mrot(i,j,isinv)*gvint2(j)
     318             :                       ELSE
     319           0 :                          vec(i) = vec(i) + oneD%ods%mrot(i,j,isinv)*gvint(j)
     320           0 :                          vec2(i) = vec2(i) + oneD%ods%mrot(i,j,isinv)*gvint2(j)
     321             :                       END IF
     322             :                    END DO
     323             :                 END DO
     324        1032 :                 DO i = 1,3
     325         432 :                    vecsum(i) = vecsum(i) + vec(i)
     326         576 :                    vecsum2(i) = vecsum2(i) + vec2(i)
     327             :                 END DO
     328             :                 !   end operator loop
     329             :              END DO
     330             :              !
     331             :              !   transform from internal to cart. coordinates
     332          24 :              starsum=MATMUL(cell%amat,vecsum)
     333          24 :              starsum2=MATMUL(cell%amat,vecsum2)
     334         192 :              DO i = 1,3
     335          72 :                 forc_a21(i) = forc_a21(i) + starsum(i)/sym%invarind(natrun)
     336          96 :                 forc_b4(i) = forc_b4(i) + starsum2(i)/sym%invarind(natrun)
     337             :              END DO
     338             :              !
     339             :              !  natrun loop end
     340             :           END DO
     341             :           !
     342             :           !     sum to existing forces
     343             :           !
     344             :           !  NOTE: force() IS REAL AND THEREFORE TAKES ONLY THE
     345             :           !  REAL PART OF forc_a21(). IN GENERAL, FORCE MUST BE
     346             :           !  REAL AFTER k-STAR SUMMATION. NOW, WE PUT THE PROPER
     347             :           !  OPERATIONS INTO REAL SPACE. PROBLEM: WHAT HAPPENS
     348             :           !  IF IN REAL SPACE THERE IS NO INVERSION ANY MORE?
     349             :           !  BUT WE HAVE INVERSION IN k-SPACE DUE TO TIME REVERSAL
     350             :           !  SYMMETRY, E(k)=E(-k)
     351             :           !  WE ARGUE THAT k-SPACE INVERSION IS AUTOMATICALLY TAKEN
     352             :           !  INTO ACCOUNT IF FORCE = (1/2)(forc_a21+conjg(forc_a21))
     353             :           !  BECAUSE TIME REVERSAL SYMMETRY MEANS THAT conjg(PSI)
     354             :           !  IS ALSO A SOLUTION OF SCHR. EQU. IF PSI IS ONE.
     355         168 :           DO i = 1,3
     356          72 :              results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a21(i) + forc_b4(i))
     357          72 :              f_a21(i,n)     = f_a21(i,n)     + forc_a21(i)
     358          96 :              f_b4(i,n)      = f_b4(i,n)      + forc_b4(i)
     359             :           END DO
     360             :           !
     361             :           !     write result moved to force_a8
     362             :           !
     363             :           !         write(*,*) a21(:,n) 
     364             :        ENDIF                                            !  IF (atoms%l_geo(n)) ...
     365          36 :        natom = natom + atoms%neq(n)
     366             :     ENDDO
     367             : 
     368          12 :     CALL timestop("force_a21")
     369             : 
     370          12 :   END SUBROUTINE force_a21
     371             : END MODULE m_forcea21

Generated by: LCOV version 1.13