Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2020 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_force_a12_lv2
8 : USE m_constants
9 : USE m_ylm
10 : USE m_sphbes
11 : USE m_gaunt
12 : USE m_types_mat
13 : !------------------------------------------------------------------------------
14 : ! Calculates the surface integral due to kinetic energy terms according to
15 : ! the last two lines of equation (14) in
16 : !
17 : ! Klüppelberg et al., PRB 91 035105 (2015).
18 : !
19 : ! Klueppelberg (force level 2)
20 : !------------------------------------------------------------------------------
21 :
22 : IMPLICIT NONE
23 :
24 : CONTAINS
25 2 : SUBROUTINE force_a12_lv2(jsp,jspd,nobd,neigd,ntypd,ntype,natd,nbasfcn,nop,&
26 2 : nvd,lmaxd,omtil,nv,neq,k1,k2,k3,invarind,invarop,&
27 2 : invtab,mrot,ngopr,amat,bmat,eig,rmt,taual,we,bkpt,&
28 2 : zMat,f_a12,force )
29 :
30 : INTEGER, INTENT (IN) :: jsp,jspd,nobd,neigd,ntypd,ntype,natd
31 : INTEGER, INTENT (IN) :: nbasfcn,nop,nvd,lmaxd
32 : REAL , INTENT (IN) :: omtil
33 :
34 : INTEGER, INTENT (IN) :: nv(jspd),neq(ntypd)
35 : INTEGER, INTENT (IN) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
36 : INTEGER, INTENT (IN) :: invarind(natd),invarop(natd,nop)
37 : INTEGER, INTENT (IN) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
38 : REAL , INTENT (IN) :: amat(3,3),bmat(3,3),eig(neigd),rmt(ntypd)
39 : REAL , INTENT (IN) :: taual(3,natd),we(nobd),bkpt(3)
40 : TYPE(t_mat), INTENT(IN) :: zMat
41 : COMPLEX, INTENT (INOUT) :: f_a12(3,ntypd)
42 : REAL , INTENT (INOUT) :: force(3,ntypd,jspd)
43 :
44 : INTEGER :: lmax,kn,iband,iatom,itype,ieq,l,m,lm,t,lp,mp,lmp,it,is
45 : INTEGER :: isinv,i
46 : REAL :: normsq,r,r2vol
47 : COMPLEX :: img,noband,fpil
48 :
49 2 : COMPLEX :: force_a12(3,ntypd),gv(3),Ygaunt(3)
50 : COMPLEX :: coeff(3,-1:1),gvint(3),vecsum(3),starsum(3)
51 2 : REAL , ALLOCATABLE :: bsl(:,:,:),G(:,:),kG(:,:),kGreal(:,:)
52 2 : REAL , ALLOCATABLE :: kineticfactor(:,:)
53 2 : COMPLEX, ALLOCATABLE :: ylm(:,:),ppw(:,:),fpw(:,:),expf(:,:)
54 2 : REAL :: zr(nbasfcn,nobd)
55 2 : COMPLEX :: zc(nbasfcn,nobd)
56 :
57 2 : CALL timestart("Force level 2")
58 :
59 2 : IF (zMat%l_real) THEN
60 126329 : zr=zMat%data_r
61 : ELSE
62 0 : zc=zMat%data_c
63 : END IF
64 :
65 2 : lmax = 2*lmaxd
66 2 : img = cmplx(0.0,1.0)
67 :
68 16 : ALLOCATE ( bsl(0:lmax,nvd,ntypd),ylm((lmax+1)**2,nvd) )
69 12 : ALLOCATE ( ppw(nobd,(lmax+1)**2),fpw(nobd,(lmax+1)**2) )
70 16 : ALLOCATE ( G(3,nvd),kG(3,nvd),kGreal(3,nvd),expf(nvd,natd) )
71 8 : ALLOCATE ( kineticfactor(nobd,nvd) )
72 :
73 2 : coeff(:, :) = cmplx(0.0,0.0)
74 2 : coeff(1,-1) = sqrt(tpi_const/3.)
75 2 : coeff(1, 1) = -sqrt(tpi_const/3.)
76 2 : coeff(2,-1) = img*sqrt(tpi_const/3.)
77 2 : coeff(2, 1) = img*sqrt(tpi_const/3.)
78 2 : coeff(3, 0) = sqrt(2.*tpi_const/3.)
79 :
80 : ! Loop over plane waves
81 2445 : DO kn = 1,nv(jsp)
82 2443 : G(1,kn) = k1(kn,jsp)
83 2443 : G(2,kn) = k2(kn,jsp)
84 2443 : G(3,kn) = k3(kn,jsp)
85 :
86 9772 : kG(:,kn) = bkpt(:) + G(:,kn)
87 39088 : kGreal(:,kn) = matmul(kG(:,kn),bmat)
88 :
89 : ! This is only useable using the Laplacian for kinetic energy
90 9772 : normsq = dot_product(kGreal(:,kn),kGreal(:,kn))
91 127036 : DO iband = 1,nobd
92 127036 : kineticfactor(iband,kn) = 0.5*normsq - eig(iband)
93 : END DO ! iband
94 :
95 2443 : CALL ylm4(lmax,kGreal(:,kn),ylm(:,kn))
96 2443 : iatom = 1
97 9774 : DO itype = 1,ntype
98 7329 : r = sqrt(normsq)*rmt(itype)
99 :
100 7329 : CALL sphbes(lmax,r,bsl(:,kn,itype))
101 39088 : DO ieq = 1,neq(itype)
102 117264 : expf(kn,iatom) = exp(tpi_const*img*dot_product(kG(:,kn),taual(:,iatom)))
103 :
104 36645 : iatom = iatom + 1
105 : END DO ! ieq
106 : END DO ! itype
107 : END DO ! kn
108 :
109 26 : force_a12 = 0.0
110 :
111 2 : iatom = 1
112 8 : DO itype = 1,ntype
113 6 : r2vol = rmt(itype)**2/omtil
114 30 : DO ieq = 1,neq(itype)
115 :
116 24 : gv = 0.0
117 360696 : ppw = 0.0
118 360696 : fpw = 0.0
119 :
120 408 : DO l = 0,lmax-1 ! (arrays run to lmax, l+1 can only be supported til lmax-1)
121 384 : fpil = 2.*tpi_const*img**l
122 : ! Calculate ppw and fpw
123 469440 : DO kn = 1,nv(jsp)
124 7973952 : DO m = -l,l
125 7504896 : lm = l*(l+1) + m + 1
126 7504896 : noband = expf(kn,iatom)*conjg(ylm(lm,kn))*bsl(l,kn,itype) * fpil
127 390723648 : DO iband = 1,nobd
128 390254592 : IF (zMat%l_real) THEN
129 382749696 : ppw(iband,lm) = ppw(iband,lm) + zr(kn,iband) * noband
130 382749696 : IF (l.gt.0) CYCLE
131 : fpw(iband,lm) = fpw(iband,lm) + zr(kn,iband) * noband * &
132 1495116 : kineticfactor(iband,kn)
133 : ELSE
134 0 : ppw(iband,lm) = ppw(iband,lm) + zc(kn,iband) * noband
135 0 : IF (l.gt.0) CYCLE
136 : fpw(iband,lm) = fpw(iband,lm) + zc(kn,iband) * noband * &
137 0 : kineticfactor(iband,kn)
138 : END IF
139 : END DO ! iband
140 : END DO ! m
141 :
142 : ! If ppw is used with l, one needs fpw with l-1 (already calculated)
143 : ! and l+1 (calculated now)
144 8912448 : DO m = -l-1,l+1
145 8443008 : lm = (l+1)*(l+2) + m + 1
146 8443008 : noband =expf(kn,iatom)*conjg(ylm(lm,kn))*bsl(l+1,kn,itype) * fpil * img
147 439505472 : DO iband = 1,nobd
148 439036416 : IF (zMat%l_real) THEN
149 : fpw(iband,lm) = fpw(iband,lm) + zr(kn,iband) * noband * &
150 430593408 : kineticfactor(iband,kn)
151 : ELSE
152 : fpw(iband,lm) = fpw(iband,lm) + zc(kn,iband) * noband * &
153 0 : kineticfactor(iband,kn)
154 : END IF
155 : END DO ! iband
156 : END DO ! m
157 : END DO ! kn
158 :
159 : ! Calculate forces
160 6552 : DO m = -l,l
161 6144 : lm = l*(l+1) + m + 1
162 18792 : DO lp = abs(l-1),l+1,2
163 55200 : DO t = -1,1
164 36792 : mp = m-t
165 36792 : IF (lp.lt.abs(mp)) CYCLE
166 34632 : lmp = lp*(lp+1) + mp + 1
167 138528 : Ygaunt(:) = gaunt2(l,1,lp,m,t,mp,lmax)*coeff(:,t)
168 1813128 : DO iband = 1,nobd
169 : gv(:) = gv(:) + we(iband) * r2vol * Ygaunt(:) * &
170 7101720 : conjg(ppw(iband,lm)) * fpw(iband,lmp)
171 : END DO ! iband
172 : END DO ! t
173 : END DO ! lp
174 : END DO ! m
175 : END DO ! lmax
176 :
177 : ! starsumgedoens
178 : ! k summation is only on IBZ. Expansion to FBZ can be achieved by
179 : ! applying the rotation of the k-points to the atomic positions of
180 : ! equivalent atoms.
181 : ! Until now, the forces for one k-point are calculated for one
182 : ! specific atom, not considering such a rotation. To do so, the
183 : ! resulting forces have to be rotated back to the representative
184 : ! atom. This is done in the next 6 lines.
185 :
186 672 : vecsum = matmul(bmat,gv)/tpi_const/neq(itype)
187 :
188 600 : gvint = matmul(mrot(:,:,1),vecsum)
189 :
190 24 : vecsum = 0
191 :
192 48 : DO it = 1,invarind(iatom) ! loop over invariant operations
193 24 : is = invarop(iatom,it)
194 24 : isinv = invtab(is)
195 :
196 696 : vecsum = vecsum + matmul(mrot(:,:,isinv),gvint(:))
197 :
198 : END DO ! it
199 :
200 24 : is = ngopr(iatom)
201 24 : it = invtab(is)
202 :
203 600 : vecsum = matmul(mrot(:,:,is),vecsum)
204 :
205 600 : starsum = matmul(amat,vecsum)
206 96 : force_a12(:,itype) = force_a12(:,itype) + real(starsum(:))/invarind(iatom)
207 :
208 30 : iatom = iatom + 1
209 : END DO ! ieq
210 24 : force(:,itype,jsp) = force(:,itype,jsp) + force_a12(:,itype)
211 26 : f_a12(:,itype) = f_a12(:,itype) + force_a12(:,itype)
212 : END DO ! itype
213 :
214 2 : DEALLOCATE ( bsl,ylm,ppw,fpw )
215 2 : DEALLOCATE ( G,kG,kGreal,expf,kineticfactor )
216 :
217 2 : CALL timestop("Force level 2")
218 :
219 2 : END SUBROUTINE force_a12_lv2
220 : END MODULE m_force_a12_lv2
|
