Line data Source code
1 : MODULE m_force_a4
2 : CONTAINS
3 29 : SUBROUTINE force_a4(atoms,sym,sphhar,input,vr,force)
4 : !--------------------------------------------------------------------------
5 : ! Core density force contribution à la Rici et al.
6 : !
7 : ! Equation A4, Phys. Rev. B 43, 6411
8 : !--------------------------------------------------------------------------
9 : USE m_types
10 : USE m_constants
11 : USE m_intgr, ONLY : intgr0,intgr3
12 : USE m_differentiate, ONLY: difcub
13 : USE m_cdn_io
14 :
15 : IMPLICIT NONE
16 :
17 : TYPE(t_atoms), INTENT(IN) :: atoms
18 : TYPE(t_sym), INTENT(IN) :: sym
19 : TYPE(t_sphhar), INTENT(IN) :: sphhar
20 : TYPE(t_input), INTENT(IN) :: input
21 :
22 : REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
23 : REAL, INTENT (INOUT) :: force(3,atoms%ntype,input%jspins)
24 :
25 : ! Local scalars
26 : COMPLEX, PARAMETER :: czero=(0.0,0.0)
27 : REAL a4_1, a4_2, qcore, s13, s23, w, xi
28 : INTEGER i, ir, jsp, lh, lindex, mem, mindex, n, nd, na
29 :
30 : ! Local arrays
31 : COMPLEX forc_a4(3),gv(3),ycomp1(3,-1:1)
32 29 : REAL rhoaux(atoms%jmtd),rhoc(atoms%jmtd,atoms%ntype,input%jspins)
33 29 : REAL tec(atoms%ntype,input%jspins),qintc(atoms%ntype,input%jspins)
34 :
35 : ! Set Ylm related components.
36 29 : s13 = SQRT(1.0/3.0)
37 29 : s23 = SQRT(2.0/3.0)
38 29 : ycomp1(1,0) = czero
39 29 : ycomp1(2,0) = czero
40 29 : ycomp1(3,0) = CMPLX(2.0*s13,0.0)
41 29 : ycomp1(1,-1) = CMPLX(s23,0.0)
42 29 : ycomp1(2,-1) = CMPLX(0.0,-s23)
43 29 : ycomp1(3,-1) = czero
44 29 : ycomp1(1,1) = CMPLX(-s23,0.0)
45 29 : ycomp1(2,1) = CMPLX(0.0,-s23)
46 29 : ycomp1(3,1) = czero
47 :
48 29 : CALL timestart("force_a4")
49 :
50 : ! Read in core density.
51 29 : CALL readCoreDensity(input,atoms,rhoc,tec,qintc)
52 :
53 58 : DO jsp = 1, input%jspins
54 118 : DO n = 1, atoms%ntype
55 60 : na = atoms%firstAtom(n)
56 89 : IF (atoms%l_geo(n)) THEN
57 60 : nd = sym%ntypsy(na)
58 :
59 240 : DO i = 1, 3
60 240 : forc_a4(i) = czero
61 : END DO
62 :
63 : ! TODO: There is no output for this. Do we want some?
64 :
65 60 : CALL intgr0(rhoc(:,n,jsp),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qcore)
66 :
67 : 8000 FORMAT (' FORCE_A4: core charge=',1p,e16.8)
68 :
69 158 : DO lh = 1, sphhar%nlh(nd)
70 158 : lindex = sphhar%llh(lh,nd)
71 :
72 158 : IF (lindex.GT.1) EXIT
73 :
74 392 : DO i = 1, 3
75 392 : gv(i) = czero
76 : END DO
77 :
78 : ! Sum over all m for a particular lattice harmonic.
79 290 : DO mem = 1, sphhar%nmem(lh,nd)
80 192 : mindex = sphhar%mlh(mem,lh,nd)
81 866 : DO i = 1, 3
82 768 : gv(i) = gv(i) + sphhar%clnu(mem,lh,nd)*ycomp1(i,mindex)
83 : END DO
84 : END DO
85 :
86 : ! Construct the integrand rhocore*d/dr(v)*r**2.
87 : ! Note: rhoc is already multiplied by r**2 and sqrt(4*pi).
88 : ! difcub performs the analytic derivative of Lagrangian of 3rd order
89 :
90 98 : xi = atoms%rmsh(1,n)
91 98 : rhoaux(1) = difcub(atoms%rmsh(1,n),vr(1,lh,n,jsp),xi)*rhoc(1,n,jsp)
92 :
93 38318 : DO ir = 2, atoms%jri(n) - 2
94 38220 : xi = atoms%rmsh(ir,n)
95 38318 : rhoaux(ir) = difcub(atoms%rmsh(ir-1,n),vr(ir-1,lh,n,jsp),xi) * rhoc(ir,n,jsp)
96 : END DO
97 :
98 98 : ir = atoms%jri(n) - 1
99 98 : xi = atoms%rmsh(ir,n)
100 98 : rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
101 :
102 98 : ir = atoms%jri(n)
103 98 : xi = atoms%rmsh(ir,n)
104 98 : rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
105 :
106 98 : CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
107 :
108 98 : a4_1 = 0.5*w/sfp_const
109 :
110 : ! Construct the integrand rhocore*v*r.
111 : ! Note: rhocore is already multiplied by r**2 and sqrt(4*pi).
112 38612 : DO ir = 1, atoms%jri(n)
113 38612 : rhoaux(ir) = rhoc(ir,n,jsp)/atoms%rmsh(ir,n)*vr(ir,lh,n,jsp)
114 : END DO
115 :
116 98 : CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
117 :
118 98 : a4_2 = w/sfp_const
119 :
120 : ! Surface contribution from non-confined core-states
121 : ! Klueppelberg Sep'12 (force level 1)
122 98 : IF (input%ctail.AND.(input%f_level.GE.1)) THEN
123 9 : w = rhoc(atoms%jri(n),n,jsp)*vr(atoms%jri(n),lh,n,jsp)
124 9 : w = 0.5*w/sfp_const
125 : ELSE
126 89 : w = 0
127 : END IF
128 :
129 648 : DO i = 1, 3
130 392 : forc_a4(i) = forc_a4(i) - (a4_1+a4_2-w)*gv(i)
131 : END DO
132 : END DO ! lh (0:sphhar%nlh(nd))
133 :
134 : ! Add onto existing forces.
135 240 : DO i = 1, 3
136 240 : force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a4(i))
137 : END DO
138 :
139 : ! Write out result.
140 60 : WRITE (oUnit,FMT=8010) n
141 60 : WRITE (oUnit,FMT=8020) (forc_a4(i),i=1,3)
142 : 8010 FORMAT (' FORCES: EQUATION A4 FOR ATOM TYPE',i4)
143 : 8020 FORMAT (' FX_A4=',2f10.6,' FY_A4=',2f10.6,' FZ_A4=',2f10.6)
144 :
145 : END IF ! atoms%l_geo(n)
146 : END DO ! n (1:atoms%ntype)
147 : END DO ! jsp (1:input%jpins)
148 :
149 29 : CALL timestop("force_a4")
150 :
151 29 : END SUBROUTINE force_a4
152 : END MODULE m_force_a4
|
