force_a21_U.f90 4.66 KB
Newer Older
1 2
MODULE m_forcea21U
CONTAINS
3
  SUBROUTINE force_a21_U(atoms,i_u,itype,isp,we,ne,usdus,v_mmp,eigVecCoeffs,aveccof,bveccof,cveccof,a21)
4 5 6 7 8 9 10
    !
    !***********************************************************************
    ! This subroutine calculates the lda+U contribution to the HF forces, 
    ! similar to the A21 term, according to eqn. (22) of F. Tran et al.
    ! Comp.Phys.Comm. 179 (2008) 784-790
    !***********************************************************************
    !
11
    USE m_constants
12 13 14
    USE m_types_setup
    USE m_types_usdus
    USE m_types_cdnval
15 16
    IMPLICIT NONE

17 18 19
    TYPE(t_usdus),INTENT(IN)        :: usdus
    TYPE(t_atoms),INTENT(IN)        :: atoms
    TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
20 21
    !     ..
    !     .. Scalar Arguments ..
22
    INTEGER, INTENT (IN)    :: itype,isp,ne
Daniel Wortmann's avatar
Daniel Wortmann committed
23 24
    INTEGER, INTENT (INOUT) :: i_u ! on input: index for the first U for atom type "itype or higher"
                                   ! on exit: index for the first U for atom type "itype+1 or higher"
25 26
    !     ..
    !     .. Array Arguments ..
27 28 29 30 31 32
    REAL,    INTENT(IN)    :: we(ne) 
    COMPLEX, INTENT(IN)    :: v_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
    COMPLEX, INTENT(IN)    :: aveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
    COMPLEX, INTENT(IN)    :: bveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
    COMPLEX, INTENT(IN)    :: cveccof(3,-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
    REAL,    INTENT(INOUT) :: a21(3,atoms%nat)
33 34 35 36 37 38 39 40 41 42 43 44 45
    !     ..
    !     .. Local Scalars ..
    COMPLEX v_a,v_b,v_c,p1,p2,p3
    INTEGER lo,lop,l,lp ,mp,lm,lmp,iatom,ie,i,m
    !     ..
    !     ..
    !*************** ABBREVIATIONS *****************************************
    ! ccof       : coefficient of the local orbital function (u_lo*Y_lm)
    ! cveccof    : is defined equivalently to aveccof, but with the LO-fct.
    ! for information on nlo,llo,uulon,dulon, and uloulopn see
    ! comments in setlomap.
    !***********************************************************************

Daniel Wortmann's avatar
Daniel Wortmann committed
46 47 48 49 50 51
    IF (atoms%lda_u(i_u)%atomType.GT.itype) RETURN

    DO WHILE (atoms%lda_u(i_u)%atomType.EQ.itype)

       l = atoms%lda_u(i_u)%l

52 53 54 55 56 57 58 59
       !
       ! Add contribution for the regular LAPWs (like force_a21, but with
       ! the potential matrix, v_mmp, instead of the tuu, tdd ...)
       !
       DO m = -l,l
          lm = l* (l+1) + m
          DO mp = -l,l
             lmp = l* (l+1) + mp
Daniel Wortmann's avatar
Daniel Wortmann committed
60 61
             v_a = v_mmp(m,mp,i_u) 
             v_b = v_mmp(m,mp,i_u) * usdus%ddn(l,itype,isp) 
62 63 64
             DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype))
                DO ie = 1,ne
                   DO i = 1,3
65 66
                      p1 = (CONJG(eigVecCoeffs%acof(ie,lm,iatom,isp)) * v_a) * aveccof(i,ie,lmp,iatom)
                      p2 = (CONJG(eigVecCoeffs%bcof(ie,lm,iatom,isp)) * v_b) * bveccof(i,ie,lmp,iatom) 
Daniel Wortmann's avatar
Daniel Wortmann committed
67 68 69 70 71 72
                      a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2) * we(ie)/atoms%neq(itype)
                   END DO
                END DO
             END DO
          END DO ! mp
       END DO   ! m
73 74 75 76 77 78

       !
       ! If there are also LOs on this atom, with the same l as
       ! the one of LDA+U, add another few terms
       !
       DO lo = 1,atoms%nlo(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
79
          IF (l == atoms%llo(lo,itype)) THEN
80 81 82 83
             DO m = -l,l
                lm = l* (l+1) + m
                DO mp = -l,l
                   lmp = l* (l+1) + mp
Daniel Wortmann's avatar
Daniel Wortmann committed
84 85 86
                   v_a = v_mmp(m,mp,i_u)
                   v_b = v_mmp(m,mp,i_u) * usdus%uulon(lo,itype,isp)
                   v_c = v_mmp(m,mp,i_u) * usdus%dulon(lo,itype,isp)
87 88 89
                   DO iatom =  sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype))
                      DO ie = 1,ne
                         DO i = 1,3
90 91 92 93 94
                            p1 = v_a * (CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * cveccof(i,mp,ie,lo,iatom))
                            p2 = v_b * (CONJG(eigVecCoeffs%acof(ie,lm,iatom,isp)) * cveccof(i,mp,ie,lo,iatom) + &
                                        CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * aveccof(i,ie,lmp,iatom))
                            p3 = v_c * (CONJG(eigVecCoeffs%bcof(ie,lm,iatom,isp)) * cveccof(i,mp,ie,lo,iatom) + &
                                        CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * bveccof(i,ie,lmp,iatom))
Daniel Wortmann's avatar
Daniel Wortmann committed
95 96 97 98 99 100 101 102
                            a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2 + p3)*we(ie)/atoms%neq(itype)
                         END DO
                      END DO
                   END DO
                END DO
             END DO
          END IF   ! l == atoms%llo(lo,itype)
       END DO     ! lo = 1,atoms%nlo
103

Daniel Wortmann's avatar
Daniel Wortmann committed
104
       i_u = i_u + 1
105
       IF(i_u.GT.atoms%n_u) EXIT
Daniel Wortmann's avatar
Daniel Wortmann committed
106
    END DO
107 108 109

  END SUBROUTINE force_a21_U
END MODULE m_forcea21U