Line data Source code
1 : MODULE m_forcea12
2 : CONTAINS
3 60 : SUBROUTINE force_a12(atoms,nobd,sym,cell ,we,jsp,ne,usdus,eigVecCoeffs, &
4 60 : acoflo,bcoflo,e1cof,e2cof,f_a12,results)
5 : !--------------------------------------------------------------------------
6 : ! Pulay 1st term force contribution à la Rici et al.
7 : !
8 : ! Equation A12, Phys. Rev. B 43, 6411
9 : !--------------------------------------------------------------------------
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 :
17 : IMPLICIT NONE
18 :
19 : TYPE(t_atoms), INTENT(IN) :: atoms
20 : TYPE(t_sym), INTENT(IN) :: sym
21 : TYPE(t_cell), INTENT(IN) :: cell
22 :
23 : TYPE(t_usdus), INTENT(IN) :: usdus
24 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs
25 : TYPE(t_results), INTENT(INOUT) :: results
26 :
27 : INTEGER, INTENT(IN) :: nobd
28 : INTEGER, INTENT(IN) :: jsp, ne
29 :
30 : REAL, INTENT(IN) :: we(nobd)
31 : COMPLEX, INTENT(IN) :: acoflo(-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
32 : COMPLEX, INTENT(IN) :: bcoflo(-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
33 : COMPLEX, INTENT(IN) :: e1cof(ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
34 : COMPLEX, INTENT(IN) :: e2cof(ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
35 : COMPLEX, INTENT(INOUT) :: f_a12(3,atoms%ntype)
36 :
37 : ! Local scalars
38 : REAL, PARAMETER :: zero=0.0
39 : COMPLEX, PARAMETER :: czero=CMPLX(0.0,0.0)
40 : COMPLEX a12,cil1,cil2
41 : INTEGER i,ie,irinv,is,isinv,it,j,l,l1,l2,lm1,lm2 ,m1,m2,n,natom,natrun,ilo,m
42 :
43 : ! Local arrays
44 : COMPLEX forc_a12(3),gv(3)
45 60 : COMPLEX acof_flapw(nobd,0:atoms%lmaxd*(atoms%lmaxd+2)),bcof_flapw(nobd,0:atoms%lmaxd*(atoms%lmaxd+2))
46 : REAL aaa(2),bbb(2),ccc(2),ddd(2),eee(2),fff(2),gvint(3),starsum(3),vec(3),vecsum(3)
47 :
48 : ! Statement functions
49 : REAL alpha,beta,delta,epslon,gamma,phi
50 : INTEGER krondel
51 :
52 : ! Kronecker delta for arguments >=0 AND <0
53 : 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
54 : alpha(l,m) = (l+1)*0.5e0*SQRT(REAL((l-m)* (l-m-1))/ REAL((2*l-1)* (2*l+1)))
55 : beta(l,m) = l*0.5e0*SQRT(REAL((l+m+2)* (l+m+1))/ REAL((2*l+1)* (2*l+3)))
56 : GAMMA(l,m) = (l+1)*0.5e0*SQRT(REAL((l+m)* (l+m-1))/ REAL((2*l-1)* (2*l+1)))
57 : delta(l,m) = l*0.5e0*SQRT(REAL((l-m+2)* (l-m+1))/ REAL((2*l+1)* (2*l+3)))
58 : epslon(l,m) = (l+1)*SQRT(REAL((l-m)* (l+m))/ REAL((2*l-1)* (2*l+1)))
59 : phi(l,m) = l*SQRT(REAL((l-m+1)* (l+m+1))/REAL((2*l+1)* (2*l+3)))
60 :
61 60 : CALL timestart("force_a12")
62 :
63 182 : DO n = 1, atoms%ntype
64 122 : natom = atoms%firstAtom(n)
65 182 : IF (atoms%l_geo(n)) THEN
66 122 : forc_a12(:) = czero
67 :
68 314 : DO natrun = natom, natom + atoms%neq(n) - 1
69 :
70 192 : gv(:) = czero
71 :
72 : ! The local orbitals do not contribute to the term a12, because
73 : ! they vanish at the MT-boundary. Therefore, the LO-contribution
74 : ! to the a and b coefficients has to be subtracted before
75 : ! calculating a12.
76 :
77 1608 : DO l1 = 0, atoms%lmax(n)
78 12168 : DO m1 = -l1, l1
79 10560 : lm1 = l1*(l1+1) + m1
80 140016 : DO ie = 1, ne
81 128040 : acof_flapw(ie,lm1) = eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)
82 138600 : bcof_flapw(ie,lm1) = eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)
83 : END DO
84 : END DO
85 : END DO
86 :
87 208 : DO ilo = 1, atoms%nlo(n)
88 16 : l1 = atoms%llo(ilo,n)
89 240 : DO m1 = -l1, l1
90 32 : lm1 = l1*(l1+1) + m1
91 1680 : DO ie = 1, ne
92 1632 : acof_flapw(ie,lm1) = acof_flapw(ie,lm1) - acoflo(m1,ie,ilo,natrun)
93 1664 : bcof_flapw(ie,lm1) = bcof_flapw(ie,lm1) - bcoflo(m1,ie,ilo,natrun)
94 : END DO
95 : END DO
96 : END DO
97 :
98 1608 : DO l1 = 0, atoms%lmax(n)
99 1416 : cil1 = ImagUnit**l1
100 12168 : DO m1 = -l1, l1
101 10560 : lm1 = l1*(l1+1) + m1
102 91728 : DO l2 = 0, atoms%lmax(n)
103 79752 : cil2 = ImagUnit**l2
104 701064 : DO m2 = -l2, l2
105 610752 : lm2 = l2*(l2+1) + m2
106 :
107 610752 : a12 = czero
108 :
109 10492776 : DO ie = 1, ne
110 :
111 : a12 = a12 + CONJG(cil1*&
112 : (acof_flapw(ie,lm1)*usdus%us(l1,n,jsp) + bcof_flapw(ie,lm1)*usdus%uds(l1,n,jsp) ))*cil2*&
113 10492776 : (e1cof(ie,lm2,natrun)*usdus%us(l2,n,jsp)+ e2cof(ie,lm2,natrun)*usdus%uds(l2,n,jsp))*we(ie)
114 :
115 : END DO
116 :
117 610752 : aaa(1) = alpha(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1+1)
118 610752 : aaa(2) = alpha(l2,m2)*krondel(l1,l2-1)* krondel(m1,m2+1)
119 610752 : bbb(1) = beta(l1,m1)*krondel(l2,l1+1)* krondel(m2,m1+1)
120 610752 : bbb(2) = beta(l2,m2)*krondel(l1,l2+1)* krondel(m1,m2+1)
121 610752 : ccc(1) = GAMMA(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1-1)
122 610752 : ccc(2) = GAMMA(l2,m2)*krondel(l1,l2-1)* krondel(m1,m2-1)
123 610752 : ddd(1) = delta(l1,m1)*krondel(l2,l1+1)* krondel(m2,m1-1)
124 610752 : ddd(2) = delta(l2,m2)*krondel(l1,l2+1)* krondel(m1,m2-1)
125 610752 : eee(1) = epslon(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1)
126 610752 : eee(2) = epslon(l2,m2)*krondel(l1,l2-1)* krondel(m1,m2)
127 610752 : fff(1) = phi(l1,m1)*krondel(l2,l1+1)* krondel(m2,m1)
128 610752 : fff(2) = phi(l2,m2)*krondel(l1,l2+1)* krondel(m1,m2)
129 :
130 : gv(1) = gv(1) + (aaa(1)+bbb(1)-ccc(1)-ddd(1)+ &
131 : aaa(2)+bbb(2)-ccc(2)-ddd(2)) * &
132 610752 : 0.5* atoms%rmt(n)**2*a12
133 :
134 : gv(2) = gv(2) + ImagUnit* (aaa(1)+bbb(1)+ccc(1)+ddd(1) - &
135 : aaa(2)-bbb(2)-ccc(2)-ddd(2)) * &
136 610752 : 0.5* atoms%rmt(n)**2*a12
137 :
138 : gv(3) = gv(3) + (eee(1)+eee(2)-fff(1)-fff(2)) * &
139 690504 : 0.5*atoms%rmt(n)**2*a12
140 :
141 : END DO ! m2 (-l2:l2)
142 : END DO ! l2 (0:atoms%lmax(n))
143 : END DO ! m1 (-l1:l1)
144 : END DO ! l1 (0:atoms%lmax(n))
145 :
146 : ! To complete the k summation over the stars now sum over all
147 : ! operations which leave (k+G)*R(natrun)*taual(natrun) invariant.
148 : ! We sum over ALL these operations and not only the ones needed
149 : ! for the actual star of k. This should be ok if we divide properly
150 : ! by the number of operations.
151 : ! First, we find the operation S where RS=T.
152 : ! T (like R) leaves the above scalar product invariant (if S=1 then
153 : ! R=T). R is the operation which generates the position of an equivalent
154 : ! atom from the position of the representative.
155 :
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 the index of operation S: is
161 :
162 : ! Transform vector gv into internal coordinates:
163 :
164 768 : vec(:) = REAL(gv(:)) /atoms%neq(n)
165 :
166 3072 : gvint=MATMUL(cell%bmat,vec)/tpi_const
167 :
168 192 : vecsum(:) = zero
169 :
170 : !-gb2002
171 : ! irinv = invtab(ngopr(natrun))
172 : ! DO it = 1,invarind(natrun)
173 : ! is = multab(irinv,invarop(natrun,it))
174 : ! ! Note, actually we need the inverse of S but -in principle
175 : !c because {S} is agroup and we sum over all S- S should also
176 : !c work; to be lucid we take the inverse:
177 : ! isinv = invtab(is)
178 : !! isinv = is
179 : ! Rotation is alreadt done in to_pulay, here we work only in the
180 : ! coordinate system of the representative atom (natom):
181 : !!
182 :
183 704 : DO it = 1, sym%invarind(natom)
184 512 : is =sym%invarop(natom,it)
185 512 : isinv = sym%invtab(is)
186 : !-gb 2002
187 : ! Now we have the wanted index of operation with which we have
188 : ! to rotate gv. Note gv is given in cart. coordinates but mrot
189 : ! acts on internal ones.
190 2048 : DO i = 1, 3
191 1536 : vec(i) = zero
192 6656 : DO j = 1, 3
193 :
194 6144 : vec(i) = vec(i) + sym%mrot(i,j,isinv)*gvint(j)
195 :
196 : END DO
197 : END DO
198 :
199 2240 : DO i = 1, 3
200 2048 : vecsum(i) = vecsum(i) + vec(i)
201 : END DO
202 :
203 : END DO ! End of operator loop
204 :
205 : ! Transform from internal to cart. coordinates:
206 2496 : starsum=MATMUL(cell%amat,vecsum)
207 890 : DO i = 1, 3
208 768 : forc_a12(i) = forc_a12(i) + starsum(i)/sym%invarind(natrun)
209 : END DO
210 :
211 : END DO ! End of natrun loop
212 :
213 : ! Add onto existing forces.
214 :
215 : ! NOTE: results%force is real and therefore only the real part of
216 : ! forc_a12 is added. In general, force must be real after the k-star
217 : ! summation. Now, we put the proper operations into real space.
218 : ! Problem: What happens if in real space there is no inversion anymore?
219 : ! But we have inversion in k-space due to time reversal symmetry:
220 : ! E(k)=E(-k)
221 : ! We argue that k-space inversion is automatically taken into account
222 : ! if force = (1/2)(forc_a12+conjg(forc_a12)), because time reversal
223 : ! symmetry means that conjg(PSI) is also a solution of Schrödinger eq.
224 : ! if PSI is one.
225 :
226 488 : DO i = 1, 3
227 366 : results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a12(i))
228 488 : f_a12(i,n) = f_a12(i,n) + forc_a12(i)
229 : END DO
230 :
231 : END IF
232 : END DO
233 :
234 : ! The result is written in force_a8.
235 :
236 60 : CALL timestop("force_a12")
237 :
238 60 : END SUBROUTINE force_a12
239 : END MODULE m_forcea12
|
