force_a4.f90 5.58 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
MODULE m_force_a4
CONTAINS
  SUBROUTINE force_a4(atoms,sphhar,input,&
       &                    vr,&
       &                    force)
    !
    ! ************************************************************
    ! rhocore force contribution a la Rici et al.
    !
    ! ************************************************************
    !
    USE m_intgr, ONLY : intgr0,intgr3
    USE m_constants, ONLY : sfp_const
    USE m_differentiate,ONLY: difcub
    USE m_types
    IMPLICIT NONE
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_sphhar),INTENT(IN)  :: sphhar
    TYPE(t_atoms),INTENT(IN)   :: atoms
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,input%jspins)
    REAL,    INTENT (INOUT) :: force(3,atoms%ntypd,input%jspins)
    !     ..
    !     .. Local Scalars ..
    REAL a4_1,a4_2,qcore,s13,s23,w,xi
    INTEGER i,ir,jsp,lh,lindex,mem,mindex,n,nd,na
    !     ..
    !     .. Local Arrays ..
    COMPLEX forc_a4(3),gv(3),ycomp1(3,-1:1)
    REAL rhoaux(atoms%jmtd),rhoc(atoms%jmtd)
    !     ..
    !     .. Data statements ..
    COMPLEX,PARAMETER:: czero=(0.000,0.000)
    !     ..
    !
    !     set ylm related components
    !
    s13 = SQRT(1.0/3.0)
    s23 = SQRT(2.0/3.0)
    ycomp1(1,0) = czero
    ycomp1(2,0) = czero
    ycomp1(3,0) = CMPLX(2.0*s13,0.0)
    ycomp1(1,-1) = CMPLX(s23,0.0)
    ycomp1(2,-1) = CMPLX(0.0,-s23)
    ycomp1(3,-1) = czero
    ycomp1(1,1) = CMPLX(-s23,0.0)
    ycomp1(2,1) = CMPLX(0.0,-s23)
    ycomp1(3,1) = czero
    !
    OPEN (17,file='cdnc',form='unformatted',status='old')
    DO  jsp = 1,input%jspins
       IF (jsp.EQ.1) THEN
          REWIND 17
       END IF
       na = 1
       DO  n = 1,atoms%ntype
          !     --->    read in core density
          READ (17) (rhoc(ir),ir=1,atoms%jri(n))
          !     --->   core density is multiplied by r**2*sqrt(4*pi)
          !     --->   such that intrg0 gives correct core charge
          !
          !     --->    read in kinetic enrgy of the core
          READ (17)
          !
          IF (atoms%l_geo(n)) THEN
             DO i = 1,3
                forc_a4(i) = czero
             END DO
             !
             CALL intgr0(rhoc,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qcore)
             !     write(16,1616) qcore
8000         FORMAT (' FORCE_A4: core charge=',1p,e16.8)
             !
             !
             nd = atoms%ntypsy(na)
             !
             !
             lh_loop: DO  lh = 1,sphhar%nlh(nd)
                lindex = sphhar%llh(lh,nd)
                IF (lindex.GT.1) EXIT lh_loop
                DO i = 1,3
                   gv(i) = czero
                END DO
                !
                !    sum over all m for particular symm. harmonic
                DO mem = 1,sphhar%nmem(lh,nd)
                   mindex = sphhar%mlh(mem,lh,nd)
                   DO i = 1,3
                      gv(i) = gv(i) + sphhar%clnu(mem,lh,nd)*ycomp1(i,mindex)
                   END DO
                END DO
                !
                !
                !     construct integrand rhocore*d/dr(v)*r**2
                !     note: rhocore is already multiplied by r**2 and srt(4.*pi)
                !     difcub performs analytic derivative of Lagrangian of 3rd order
                !
                xi = atoms%rmsh(1,n)
                rhoaux(1) = difcub(atoms%rmsh(1,n),vr(1,lh,n,jsp),xi)*rhoc(1)
                DO ir = 2,atoms%jri(n) - 2
                   xi = atoms%rmsh(ir,n)
                   rhoaux(ir) = difcub(atoms%rmsh(ir-1,n),&
                        &                                vr(ir-1,lh,n,jsp),xi) * rhoc(ir)
                END DO
                !
                ir = atoms%jri(n) - 1
                xi = atoms%rmsh(ir,n)
                rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
                     &                              vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir)
                !
                ir = atoms%jri(n)
                xi = atoms%rmsh(ir,n)
                rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
                     &                              vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir)
                CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
                a4_1 = 0.5*w/sfp_const
                !
                !     construct integrand rhocore*v*r
                !     note: rhocore is already multiplied by r**2 and srt(4.*pi)
                !
                DO ir = 1,atoms%jri(n)
                   rhoaux(ir) = rhoc(ir)/atoms%rmsh(ir,n)*vr(ir,lh,n,jsp)
                END DO
                !
                CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
                a4_2 = w/sfp_const
                !
                DO i = 1,3
                   forc_a4(i) = forc_a4(i) - (a4_1+a4_2)*gv(i)
                END DO
                !
                !  lh loop ends
             ENDDO lh_loop
             !
             !     sum to existing forces
             !
             DO i = 1,3
                force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a4(i))
             END DO
             !
             !     write result
             !
             WRITE (6,FMT=8010) n
             WRITE (16,FMT=8010) n
             WRITE (6,FMT=8020) (forc_a4(i),i=1,3)
             WRITE (16,FMT=8020) (forc_a4(i),i=1,3)
8010         FORMAT (' FORCES: EQUATION A4 FOR ATOM TYPE',i4)
8020         FORMAT (' FX_A4=',2f10.6,' FY_A4=',2f10.6,' FZ_A4=',2f10.6)
             ! type loop ends
          ENDIF
          na = na + atoms%neq(n)
       ENDDO
       READ (17)
       ! spin loop ends
    ENDDO
    !
    CLOSE (17)
    RETURN
  END SUBROUTINE force_a4
END MODULE m_force_a4