Branch data Line data Source code
1 : : !!****m* ABINIT/m_nonlop_pl
2 : : !! NAME
3 : : !! nonlop_pl
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, GZ, MT, FF, DRH)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! PARENTS
14 : : !! m_nonlop
15 : : !!
16 : : !! CHILDREN
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : module m_nonlop_pl
27 : :
28 : : use defs_basis
29 : : use m_errors
30 : : use m_abicore
31 : : use m_xmpi
32 : : use m_contistr01
33 : : use m_contistr03
34 : : use m_contistr12
35 : : use m_contstr21
36 : : use m_contstr23
37 : : use m_contstr25
38 : : use m_contstr25a
39 : : use m_contstr26
40 : : use m_metstr
41 : : use m_opernl
42 : :
43 : : use defs_abitypes, only : MPI_type
44 : : use m_geometry, only : strconv
45 : : use m_kg, only : ph1d3d
46 : : use m_contract, only : cont22cso, cont22so, cont24, cont33cso, cont33so, cont35, cont22, cont3, cont13, &
47 : : metcon, metcon_so, metric_so
48 : : use m_numeric_tools, only : geteuler
49 : : implicit none
50 : :
51 : : private
52 : : !!***
53 : :
54 : : public :: nonlop_pl
55 : : !!***
56 : :
57 : : contains
58 : : !!***
59 : :
60 : : !!****f* ABINIT/nonlop_pl
61 : : !! NAME
62 : : !! nonlop_pl
63 : : !!
64 : : !! FUNCTION
65 : : !! * Compute application of a nonlocal operator Vnl in order to get:
66 : : !! - contracted elements (energy, forces, stresses, ...), if signs=1
67 : : !! - a function in reciprocal space (|out> = Vnl|in>), if signs=2
68 : : !! Operator Vnl, as the following form:
69 : : !! $Vnl=sum_{R,lmn,l''m''n''} {|P_{Rlmn}> Enl^{R}_{lmn,l''m''n''} <P_{Rl''m''n''}|}$
70 : : !! Operator Vnl is -- in the typical case -- the nonlocal potential.
71 : : !! - With norm-conserving pseudopots, $Enl^{R}_{lmn,l''m''n''}$ is the
72 : : !! Kleinmann-Bylander energy $Ekb^{R}_{ln}$.
73 : : !! - The |P_{Rlmn}> are the projector functions.
74 : : !! * This routine uses Legendre polynomials Pl to express Vnl.
75 : : !!
76 : : !! INPUTS
77 : : !! choice: chooses possible output:
78 : : !! choice=1 => a non-local energy contribution
79 : : !! =2 => a gradient with respect to atomic position(s)
80 : : !! =3 => a gradient with respect to strain(s)
81 : : !! =23=> a gradient with respect to atm. pos. and strain(s)
82 : : !! =4 => a gradient and 2nd derivative with respect to atomic pos.
83 : : !! =5 => a gradient with respect to k wavevector
84 : : !! =6 => 2nd derivatives with respect to strain
85 : : !! dimekb1,dimekb2=dimensions of ekb (see ekb)
86 : : !! dimffnlin=second dimension of ffnlin (1+number of derivatives)
87 : : !! dimffnlout=second dimension of ffnlout (1+number of derivatives)
88 : : !! ekb(dimekb1,dimekb2,nspinortot**2)= (Real) Kleinman-Bylander energies (hartree)
89 : : !! dimekb1=lmnmax - dimekp2=ntypat
90 : : !! ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used
91 : : !! for the application of the nonlocal operator to the |in> vector
92 : : !! ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used
93 : : !! for the application of the nonlocal operator to the |out> vector
94 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
95 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
96 : : !! (bohr^-1)
97 : : !! idir=direction of the - atom to be moved in the case (choice=2,signs=2),
98 : : !! - k point direction in the case (choice=5,signs=2)
99 : : !! - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
100 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
101 : : !! istwf_k=option parameter that describes the storage of wfs
102 : : !! kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector
103 : : !! kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector
104 : : !! kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector
105 : : !! kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector
106 : : !! kptin(3)=k point in terms of recip. translations, for the |in> vector
107 : : !! kptout(3)=k point in terms of recip. translations, for the |out> vector
108 : : !! lmnmax=max. number of (l,m,n) components over all types of atoms
109 : : !! matblk=dimension of the arrays ph3din and ph3dout
110 : : !! mgfft=maximum size of 1D FFTs
111 : : !! mpi_enreg=information about MPI parallelization
112 : : !! natom=number of atoms in cell
113 : : !! nattyp(ntypat)=number of atoms of each type
114 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
115 : : !! nkpgin,nkpgout=second sizes of arrays kpgin/kpgout
116 : : !! nloalg(3)=governs the choice of the algorithm for nonlocal operator
117 : : !! nnlout=dimension of enlout: choice=1=>nnlout=1 choice=2=>nnlout=3*natom
118 : : !! choice=3=>nnlout=6 choice=4=>nnlout=6*natom
119 : : !! choice=5=>nnlout=1 choice=6=>nnlout=6*(3*natom+6)
120 : : !! choice=23=>nnlout=6+3*natom
121 : : !! npwin=number of planewaves for given k point, for the |in> vector
122 : : !! npwout=number of planewaves for given k point, for the |out> vector
123 : : !! nspinor=number of spinorial components of the wavefunctions on current proc
124 : : !! nspinortot=total number of spinorial components of the wavefunctions
125 : : !! ntypat=number of types of atoms in cell
126 : : !! only_SO=flag to calculate only the SO part in nonlop
127 : : !! phkxredin(2,natom)=phase factors exp(2 pi kptin.xred)
128 : : !! phkxredout(2,natom)=phase factors exp(2 pi kptout.xred)
129 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
130 : : !! ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in)
131 : : !! ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out)
132 : : !! --- pspso removed in beautification because it was unused ---
133 : : !! pspso(ntypat)=spin-orbit characteristic for each atom type
134 : : !! -------------------------------------------------------------
135 : : !! signs= if 1, get contracted elements (energy, forces, stress, ...)
136 : : !! if 2, applies the non-local operator to a function in reciprocal space
137 : : !! ucvol=unit cell volume (bohr^3)
138 : : !! vectin(2,nspinor*npwin)=input cmplx wavefunction coefficients <G|Cnk>
139 : : !! use_gbt= if 1, no spin-orbit coupling (scaler-Relativistic only);
140 : : !! if 2, include only the \sigma_z component of the spin-orbit coupling.
141 : : !!
142 : : !! OUTPUT
143 : : !! ==== if (signs==1) ====
144 : : !! enlout(nnlout)= contribution of this state to the nl part
145 : : !! of the following properties:
146 : : !! if choice=1 : enlout(1) -> the energy
147 : : !! if choice=2 : enlout(1:3*natom) -> the forces
148 : : !! if choice=3 : enlout(1:6) -> the stresses
149 : : !! if choice=23: enlout(1:6+3*natom) -> the forces and the stresses
150 : : !! if choice=4 : enlout(1:6*natom) -> the frozen wf part of dynam. matrix
151 : : !! if choice=6 : enlout(1:6*(3*natom+6)) -> the frozen wf part of elastic tensor
152 : : !! ==== if (signs==2) ====
153 : : !! vectout(2,nspinor*npwout)= result of the aplication of the nl operator
154 : : !! or one of its derivative to the input vect.:
155 : : !! if choice=1 : Vnl |vectin>
156 : : !! if choice=2 : dVnl/d(xred(idir,iatom) |vectin> (xred=reduced atm. pos.)
157 : : !! if choice=3 : dVnl/d(strain(idir)) |vectin> (symmetric strain =>idir=1...6)
158 : : !! if choice=5 : dVnl/dk(idir) |vectin> (k wavevector)
159 : : !!
160 : : !! NOTES
161 : : !! In the case signs=1, the array vectout is not used, nor modified
162 : : !! so that the same array as vectin can be used as a dummy argument;
163 : : !! the same is true for the pairs npwin-npwout, ffnlin-ffnlout,
164 : : !! kgin-kgout, ph3din-ph3dout, phkredin-phkxredout).
165 : : !!
166 : : !! Calculation includes contributions to grads of Etot wrt coord and
167 : : !! wrt strains for l=0,1,2,3.
168 : : !!
169 : : !! WARNINGS
170 : : !! - Warning 1: This routine is in a transient state, during the
171 : : !! time of the implementation of the spin-orbit coupling...
172 : : !! In particular, the OMP parallelisation is still missing,
173 : : !! but it matters here only when nspinor==2.
174 : : !! - Warning 2: the order of atoms is governed by atindx
175 : : !!
176 : : !! PARENTS
177 : : !! nonlop
178 : : !!
179 : : !! CHILDREN
180 : : !! cont13,cont22,cont22cso,cont22so,cont24,cont3,cont33cso,cont33so,cont35
181 : : !! contistr01,contistr03,contistr12,contstr21,contstr23,contstr25
182 : : !! contstr25a,contstr26,ddkten,metcon,metcon_so,metric_so,metstr,opernl2
183 : : !! opernl3,opernl4a,opernl4b,ph1d3d,scalewf_nonlop,strconv,strsocv,trace2
184 : : !! xmpi_sum
185 : : !!
186 : : !! SOURCE
187 : :
188 : 88372944 : subroutine nonlop_pl(choice,dimekb1,dimekb2,dimffnlin,dimffnlout,ekb,enlout,&
189 : 58915296 : & ffnlin,ffnlout,gmet,gprimd,idir,indlmn,ispin_gbt,istwf_k,kgin,kgout,kpgin,kpgout,&
190 : : & kptin,kptout,lmnmax,matblk,mgfft,mpi_enreg,mpsang,mpssoang,&
191 [ - + ]: 29457648 : & natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,npwin,npwout,nspinor,nspinortot,&
192 : 29457648 : & ntypat,only_SO,phkxredin,phkxredout,ph1d,ph3din,ph3dout,signs,&
193 [ - + - + ]: 29457648 : & spinaxis,ucvol,use_gbt,vectin,vectout)
194 : :
195 : : !Arguments ------------------------------------
196 : : !This type is defined in defs_mpi
197 : : !The (inout) classification below is misleading; mpi_enreg is temporarily
198 : : ! changed but reset to its initial condition before exiting.
199 : : !scalars
200 : : integer,intent(in) :: choice,dimekb1,dimekb2,dimffnlin,dimffnlout,idir,istwf_k,ispin_gbt
201 : : integer,intent(in) :: lmnmax,matblk,mgfft,mpsang,mpssoang,natom,nkpgin,nkpgout
202 : : integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,only_SO,signs,use_gbt
203 : : real(dp),intent(in) :: ucvol
204 : : type(MPI_type),intent(in) :: mpi_enreg
205 : : !arrays
206 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),kgin(3,npwin),kgout(3,npwout)
207 : : integer,intent(in) :: nattyp(ntypat),ngfft(18),nloalg(3) !,pspso(ntypat) UNUSED
208 : : real(dp),intent(in) :: ekb(dimekb1,dimekb2,nspinortot**2)
209 : : real(dp),intent(in) :: ffnlin(npwin,dimffnlin,lmnmax,ntypat)
210 : : real(dp),intent(in) :: ffnlout(npwout,dimffnlout,lmnmax,ntypat),gmet(3,3)
211 : : real(dp),intent(in) :: gprimd(3,3),kpgin(npwin,nkpgin),kpgout(npwout,nkpgout)
212 : : !real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom) !vz_d
213 : : real(dp),intent(in) :: kptin(3),kptout(3) !vz_d
214 : : real(dp),intent(in) :: ph1d(2,*),spinaxis(3) !vz_d
215 : : real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom)
216 : : real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk)
217 : : real(dp),intent(inout) :: vectin(:,:)
218 : : real(dp),intent(out) :: enlout(:) !vz_i
219 : : real(dp),intent(inout) :: vectout(:,:) !vz_i
220 : :
221 : : !Local variables-------------------------------
222 : : !mlang is the maximum number of different angular momenta
223 : : !(mlang=4 means s,p,d,f)
224 : : ! Note : in a future version, one should adjust mlang to mpsang.
225 : : !mlang2 is the maximum number of unique tensor components for a tensor
226 : : !of rank (mlang-1) with index range 1-3
227 : : !mlang3 is the maximum number of unique tensor components summed over
228 : : !all tensors of rank 0 through mlang-1.
229 : : !mlang4 is the total number of additional unique tensor components
230 : : !related to strain gradients, ranks 2 through mlang+1.
231 : : !mlang6 is the total number of additional unique tensor components
232 : : !related to strain 2nd derivaives, ranks 4 through mlang+3.
233 : : !mlang1 is the total number of certain additional unique tensor components
234 : : !related to internal strain, ranks 1 through mlang
235 : : !mlang5 is the total number of other additional unique tensor components
236 : : !related to internal strain, ranks 1 through mlang
237 : : !scalars
238 : : integer,parameter :: mlang=4
239 : : ! MG: I tried to use parameters instead of saved variables but [tutorespfn][trf2_1] gets stuck on milou_g95_snofbfpe
240 : : integer,save :: mlang1=((mlang+1)*(mlang+2)*(mlang+3))/6-1
241 : : !integer,save :: mlang2=(mlang*(mlang+1))/2 ! Unused
242 : : integer,save :: mlang3=(mlang*(mlang+1)*(mlang+2))/6
243 : : integer,save :: mlang4=((mlang+2)*(mlang+3)*(mlang+4))/6-4
244 : : integer,save :: mlang5=((mlang+3)*(mlang+4)*(mlang+5))/6-10
245 : : integer,save :: mlang6=((mlang+4)*(mlang+5)*(mlang+6))/6-20
246 : : integer :: compact,ia,ia1,ia2,ia3,ia4,ia5,ierr,iest,ig,ii,ilang,ilang2,ilmn
247 : : integer :: iln,iln0,indx,iproj,ipsang,ishift,isp,ispin,ispinor,ispinor_index,ispinp,ispinor_ekb
248 : : integer :: istr,istr1,istr2,iterm,itypat,jj,jjk,jjs,jjs1,jjs2,jjs3,jjs4,jjstr,jspin
249 : : integer :: mincat,mproj,mu,mumax,n1,n2,n3,ndgxdt,ndgxdtfac,nincat,nlang
250 : : integer :: nproj,nspinso,rank
251 : : integer :: sign,spaceComm, isft
252 : : real(dp) :: alpha,beta,e2nl,e2nldd,enlk
253 : : character(len=500) :: msg
254 : : !arrays
255 : 29457648 : integer,allocatable :: indlmn_s(:,:,:),jproj(:)
256 : : real(dp) :: amet(2,3,3,2,2),amet_lo(3,3),e2nl_tmp(6),eisnl(3),rank2(6)
257 : : real(dp) :: rank2c(2,6),soc_weight(3),strsnl(6),strsnl_out(6),strsso(6,3),strssoc(6),trace(2)!,tsec(2)
258 : 29457648 : real(dp),allocatable :: d2gxdis(:,:,:,:,:),d2gxdis_s(:,:,:,:)
259 : 29457648 : real(dp),allocatable :: d2gxds2(:,:,:,:,:),d2gxds2_s(:,:,:,:)
260 : 29457648 : real(dp),allocatable :: dgxdis(:,:,:,:,:),dgxdis_s(:,:,:,:),dgxds(:,:,:,:,:)
261 : 29457648 : real(dp),allocatable :: dgxds_s(:,:,:,:),dgxdsfac(:,:,:,:,:)
262 : 29457648 : real(dp),allocatable :: dgxdt(:,:,:,:,:,:),dgxdt_s(:,:,:,:,:)
263 : 29457648 : real(dp),allocatable :: dgxdtfac(:,:,:,:,:),ekb_s(:,:),gxa(:,:,:,:,:)
264 : 29457648 : real(dp),allocatable :: gxa_s(:,:,:,:),gxafac(:,:,:,:),pauli(:,:,:,:)
265 : 29457648 : real(dp),allocatable :: temp(:,:),tmpfac(:,:),vectin_s(:,:),vectout_s(:,:)
266 : 29457648 : real(dp),allocatable :: wt(:,:)
267 : :
268 : : ! **********************************************************************
269 : :
270 : : ABI_UNUSED(mgfft)
271 : :
272 : : !Test: spin orbit not allowed for choice=5,6
273 [ + + - + ]: 29457648 : if (nspinortot==2 .and. choice==6 ) then
274 : 0 : ABI_BUG('For nspinortot=2, choice=6 is not yet allowed.')
275 : : end if
276 : :
277 [ - + ]: 29457648 : if ((choice<1 .or. choice>6) .and. choice/=23 ) then
278 : 0 : write(msg,'(a,i0)')'Does not presently support this choice=',choice
279 : 0 : ABI_BUG(msg)
280 : : end if
281 : :
282 : : !Test: choice 51 and 52 only allowed with nonlop_ylm
283 : : !JWZ, 01-Sep-08
284 [ - + ]: 29457648 : if (choice==51 .or. choice==52) then
285 : 0 : ABI_BUG('choice 51 or 52 is not yet allowed.')
286 : : end if
287 : :
288 : : !Define dimension of work arrays.
289 [ + + ]: 71836647 : mincat=min(NLO_MINCAT,maxval(nattyp))
290 [ + + + + ]: 270290592 : mproj=maxval(indlmn(3,:,:))
291 [ - + + - : 88372944 : ABI_MALLOC(temp,(2,mlang4))
+ - + - -
+ - + ]
292 [ - + + - : 58915296 : ABI_MALLOC(tmpfac,(2,mlang4))
- + - + ]
293 [ + + + - : 88372944 : ABI_MALLOC(wt,(mlang,mproj))
+ - + + -
+ - + ]
294 [ - + ]: 29457648 : ABI_MALLOC(jproj,(mlang))
295 : 29457648 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
296 [ + + + - : 58915296 : ABI_MALLOC(ekb_s,(mlang,mproj))
- + - + ]
297 [ - + + - : 117830592 : ABI_MALLOC(indlmn_s,(6,lmnmax,ntypat))
- + + - +
- + - - +
- + ]
298 : :
299 : : !Eventually compute the spin-orbit metric tensor:
300 [ + + ]: 29457648 : if (mpssoang>mpsang) then
301 [ + - ]: 2339994 : ABI_MALLOC(pauli,(2,2,2,3))
302 [ + + ]: 9359976 : soc_weight = one
303 : : ! GBT: keep only the z-component, optionally flip its sign
304 [ + + ]: 2339994 : if (use_gbt == 2) then
305 [ + + ]: 47388 : soc_weight(1:2) = 0
306 [ + + ]: 15796 : if (ispin_gbt == 2) soc_weight(3) = -1
307 : : end if
308 : 2339994 : call geteuler(spinaxis,alpha,beta)
309 : 2339994 : call metric_so(amet,soc_weight,gprimd,pauli,alpha,beta)
310 : : end if
311 : :
312 : : !Allocate array gxa (contains projected scalars).
313 [ - + + - : 176745888 : ABI_MALLOC(gxa,(2,mlang3,mincat,mproj,nspinortot))
- + + - +
+ + - - +
+ - + - +
+ - + -
+ ]
314 [ + + ]: 29457648 : if(nspinor==2) then
315 [ - + + - : 7760870 : ABI_MALLOC(gxa_s,(2,mlang3,mincat,mproj))
- + + - +
+ + - + -
+ + - + -
+ ]
316 : : else
317 [ - + ]: 27905474 : ABI_MALLOC(gxa_s,(0,0,0,0))
318 : : end if
319 : :
320 [ - + + - : 147288240 : ABI_MALLOC(gxafac,(2,mlang3,mincat,mproj))
- + + - +
+ + - + -
+ + - + +
- ]
321 [ + + + + : 5618103892 : gxa(:,:,:,:,:)=zero
+ + + + +
+ ]
322 : :
323 : : !If choice==2 : first-order atomic displacements
324 : : !If signs==2, only one direction is considered
325 : : !If signs==1, the three directions are considered
326 : : !If choice==4 and signs==1 : second-order atomic displacements,
327 : : !the nine components are considered
328 : : !If choice==5 and signs==2 : ddk
329 : : !component 1 -> from ffnl(:,2,...)
330 : : !component 2 -> from ffnl(:,1,...) (too much space is booked for this
331 : : !case, since the number of angular momenta is smaller than mlang3, but it is easier)
332 : 29457648 : ndgxdt=0
333 [ + + + + ]: 29457648 : if(signs==2 .and. choice==2) ndgxdt=1
334 [ + + + + ]: 29457648 : if(signs==1 .and. (choice==2.or.choice==23)) ndgxdt=3
335 [ + + ]: 29457648 : if(choice==4) ndgxdt=9
336 [ + + ]: 29457648 : if(choice==5) ndgxdt=2
337 : : !Allocate dgxdt (contains derivatives of gxa with respect to atomic displacements or ddk).
338 [ + + + - : 206203536 : ABI_MALLOC(dgxdt,(2,ndgxdt,mlang3,mincat,mproj,nspinortot))
- + + - -
+ + - + +
+ - - + +
- + - + +
- + - + ]
339 [ + + + + : 3068408272 : dgxdt(:,:,:,:,:,:)=zero
+ + + + +
+ + + ]
340 [ + + ]: 29457648 : if(nspinor==2)then
341 [ + + + - : 9313044 : ABI_MALLOC(dgxdt_s,(2,ndgxdt,mlang3,mincat,mproj))
- + + - -
+ + - + +
+ - + - +
+ - + +
- ]
342 [ + + + + : 432512190 : dgxdt_s(:,:,:,:,:)=zero
+ + + + +
+ ]
343 : : else
344 [ - + ]: 27905474 : ABI_MALLOC(dgxdt_s,(0,0,0,0,0))
345 : : end if
346 : 29457648 : ndgxdtfac=0
347 [ + + + + ]: 29457648 : if(signs==2 .and. choice==2) ndgxdtfac=1
348 [ + + ]: 29457648 : if(choice==4) ndgxdtfac=3
349 [ + + ]: 29457648 : if(choice==5) ndgxdtfac=2
350 [ + + + - : 176745888 : ABI_MALLOC(dgxdtfac,(2,ndgxdtfac,mlang3,mincat,mproj))
- + + - -
+ + - + +
+ - + - +
+ - + -
+ ]
351 : :
352 : : !Allocate dgxds (contains derivatives of gxa with respect to strains).
353 [ - + + - : 176745888 : ABI_MALLOC(dgxds,(2,mlang4,mincat,mproj,nspinor))
- + + - +
+ + - - +
+ - + - +
+ - + -
+ ]
354 [ + + + + : 14275266580 : dgxds(:,:,:,:,:)=zero
+ + + + +
+ ]
355 [ - + + - : 147288240 : ABI_MALLOC(dgxdsfac,(2,mlang4,mincat,mproj,nspinor))
- + + - +
+ + - - +
+ - - + -
+ ]
356 [ + + ]: 29457648 : if(choice==6) then
357 [ - + + - : 126192 : ABI_MALLOC(dgxdis,(2,mlang1,mincat,mproj,nspinor))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
358 [ - + + - : 126192 : ABI_MALLOC(d2gxdis,(2,mlang5,mincat,mproj,nspinor))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
359 [ - + + - : 126192 : ABI_MALLOC(d2gxds2,(2,mlang6,mincat,mproj,nspinor))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
360 : : else
361 [ - + ]: 29436616 : ABI_MALLOC(dgxdis ,(0,0,0,0,0))
362 [ - + ]: 29436616 : ABI_MALLOC(d2gxdis,(0,0,0,0,0))
363 [ - + ]: 29436616 : ABI_MALLOC(d2gxds2,(0,0,0,0,0))
364 : : end if
365 [ - + ]: 29457648 : ABI_MALLOC(dgxds_s ,(0,0,0,0))
366 [ - + ]: 29457648 : ABI_MALLOC(dgxdis_s ,(0,0,0,0))
367 [ - + ]: 29457648 : ABI_MALLOC(d2gxdis_s,(0,0,0,0))
368 [ - + ]: 29457648 : ABI_MALLOC(d2gxds2_s,(0,0,0,0))
369 [ + + ]: 29457648 : if(nspinor==2)then
370 [ - + ]: 1552174 : ABI_FREE(dgxds_s)
371 [ - + + - : 7760870 : ABI_MALLOC(dgxds_s,(2,mlang4,mincat,mproj))
- + + - +
+ + - + -
+ + - + +
- ]
372 [ + + + + : 2410320282 : dgxds_s(:,:,:,:)=zero
+ + + + ]
373 [ - + ]: 1552174 : if(choice==6) then
374 [ # # ]: 0 : ABI_FREE(dgxdis_s)
375 [ # # ]: 0 : ABI_FREE(d2gxdis_s)
376 [ # # ]: 0 : ABI_FREE(d2gxds2_s)
377 [ # # # # : 0 : ABI_MALLOC(dgxdis_s,(2,mlang1,mincat,mproj))
# # # # #
# # # # #
# # # # #
# ]
378 [ # # # # : 0 : ABI_MALLOC(d2gxdis_s,(2,mlang5,mincat,mproj))
# # # # #
# # # # #
# # # # #
# ]
379 [ # # # # : 0 : ABI_MALLOC(d2gxds2_s,(2,mlang6,mincat,mproj))
# # # # #
# # # # #
# # # # #
# ]
380 : : else
381 : : end if
382 : : end if
383 : :
384 : : !Zero out some arrays
385 [ + + + + ]: 47247716 : if(choice==2 .or. choice==4 .or. choice==5 .or. choice==6 .or. choice==23) enlout(:)=0.0d0
386 : :
387 [ + + + + : 28490720499 : if(signs==2) vectout(:,:)=zero
+ + ]
388 : :
389 : : !if(choice==3.or.choice==23) then
390 : : ! strsnl(:)=zero
391 : : ! if(mpssoang>mpsang) strsso(:,:)=zero
392 : : !end if
393 : 29457648 : enlk=zero
394 : 29457648 : strsso = zero
395 : 29457648 : strsnl = zero
396 : :
397 : : !In the case vectin is a spinor, split its second part.
398 : : !Also, eventually take into account the storage format of the wavefunction
399 : : !(the original vector will be restored before leaving the routine,
400 : : !except for the vectin(2,1) component with istwf_k==2,
401 : : !that should vanish)
402 : : !In sequential, treat the second spinor part first
403 [ + + ]: 29457648 : if (nspinor==2)then
404 [ - + + - : 4656522 : ABI_MALLOC(vectin_s,(2,npwin))
+ - + - -
+ - + ]
405 [ - + + - : 4656522 : ABI_MALLOC(vectout_s,(2,npwout))
+ - + - -
+ - + ]
406 : :
407 [ - + ]: 1552174 : isft = npwin;if (mpi_enreg%nproc_spinor>1) isft=0
408 : :
409 : : ! Initialize it
410 : : !$OMP PARALLEL DO
411 [ + + ]: 329881518 : do ig=1,npwin
412 : 328329344 : vectin_s(1,ig)=vectin(1,ig+isft)
413 : 329881518 : vectin_s(2,ig)=vectin(2,ig+isft)
414 : : end do
415 : :
416 : : ! Take into account the storage
417 [ - + ]: 1552174 : if(istwf_k/=1)then
418 : 0 : call scalewf_nonlop(istwf_k,mpi_enreg,npwin,1,vectin_s)
419 : : end if
420 : : end if ! nspinortot==2
421 : :
422 : : !Treat the first spinor part now
423 [ + + ]: 29457648 : if(istwf_k/=1) then
424 [ - + - - : 904765 : call scalewf_nonlop(istwf_k,mpi_enreg,npwin,1,vectin)
- - - - -
+ - - -
- ]
425 : : end if
426 : :
427 : :
428 : : !Big loop on atom types.
429 : 29457648 : ia1=1
430 [ + + ]: 71836647 : do itypat=1,ntypat
431 : :
432 : : ! Get atom loop indices for different types:
433 : 42378999 : ia2=ia1+nattyp(itypat)-1
434 : :
435 : : ! Cut the sum on different atoms in blocks, to allow memory saving.
436 : : ! Inner summations on atoms will be done from ia3 to ia4.
437 : : ! Note that the maximum range from ia3 to ia4 is mincat (maximum
438 : : ! increment of atoms).
439 [ - + - - : 85163057 : do ia3=ia1,ia2,mincat
+ + + + ]
440 : 42784058 : ia4=min(ia2,ia3+mincat-1)
441 : : ! Give the increment of number of atoms in this subset.
442 : 42784058 : nincat=ia4-ia3+1
443 : :
444 : : ! Prepare the phase factors for the atoms between ia3 and ia4 :
445 : : ! For nloalg(2)<=0, they were not prepared previously,it is needed to
446 : : ! compute them again.
447 [ + + ]: 42784058 : if(nloalg(2)<=0)then
448 : : ! For nloalg(2)==0, it is needed to compute the phase factors.
449 [ - + ]: 2243368 : if(mincat>matblk)then
450 : : write(msg,'(a,a,a,i4,a,i4,a)')&
451 : 0 : 'With nloc_mem<=0, mincat must be less than matblk.',ch10,&
452 : 0 : 'Their value is ',mincat,' and ',matblk,'.'
453 : 0 : ABI_BUG(msg)
454 : : end if
455 : 2243368 : call ph1d3d(ia3,ia4,kgin,matblk,natom,npwin,n1,n2,n3,phkxredin,ph1d,ph3din)
456 : : end if
457 : :
458 : : ! Here begins the different treatment for the scalar-relativistic
459 : : ! part(s) and the spin-orbit part.
460 : : ! Loop on ispinor : 1 for scalar-Relativistic, 2 for spin-orbit
461 [ + + ]: 42784058 : nspinso=1;if (mpssoang>mpsang) nspinso=2
462 : :
463 : : ! Change nspinso if collinear run or if nspinor == 2 and SOC is not wanted.
464 : : ! TODO: The last check requires pspso
465 [ + + + + ]: 42784058 : if (nspinortot == 1 .and. use_gbt /= 2) nspinso = 1
466 : :
467 [ + + ]: 130403329 : do ispinor=1,nspinso
468 : 45240272 : ispinor_index=ispinor
469 [ + + ]: 45240272 : if (mpi_enreg%paral_spinor==1) ispinor_index=mpi_enreg%me_spinor+1
470 : :
471 : : !
472 : : ! mjv 25 6 2008: if only_SO == 1 or 2 skip the scalar relativistic terms
473 : : ! only output the spin orbit ones
474 : : !
475 [ + + + - ]: 45240272 : if (ispinor==1 .and. only_SO>0) cycle
476 : :
477 : : ! Select scalar-relativistic or spin-orbit KB-energies:
478 [ + + + + : 901451054 : ekb_s(:,:)=zero ; wt(:,:)=zero
+ + + + ]
479 : : ! Loop over (l,n) values (loop over l,m,n and test on l,n)
480 [ + + ]: 226201360 : iln0=0 ; jproj(:)=0 ; nlang=0
481 [ + + + + ]: 1675401805 : indlmn_s(:,:,itypat)=0
482 [ + + ]: 278120491 : do ilmn=1,lmnmax
483 [ + + ]: 232880219 : if(ispinor/=indlmn(6,ilmn,itypat))cycle
484 [ + + ]: 1218826266 : indlmn_s(:,ilmn,itypat)=indlmn(:,ilmn,itypat)
485 : 174118038 : iln =indlmn(5,ilmn,itypat)
486 [ - + ]: 219358310 : if (iln>iln0) then
487 : 174118038 : iln0=iln
488 : 174118038 : ipsang=indlmn(1,ilmn,itypat)+1
489 : : ! DEBUG
490 : : ! write(std_out,*)' nonlop : ipsang,ilmn,itypat,ispinor=',ipsang,ilmn,itypat,ispinor
491 : : ! ENDDEBUG
492 : 174118038 : iproj=indlmn(3,ilmn,itypat)
493 : : ! This shift is not needed anymore
494 : : ! if (ispinor==2) ipsang=indlmn(1,ilmn,itypat)-mpsang+2
495 : : ! ekb_s(ipsang,iproj)=ekb(iln,itypat,ispinor)
496 [ + - ]: 174118038 : ispinor_ekb=min(ispinor,ubound(ekb,3))
497 : 174118038 : ekb_s(ipsang,iproj)=ekb(iln,itypat,ispinor_ekb)
498 : 174118038 : wt(ipsang,iproj)=4.d0*pi/ucvol*dble(2*ipsang-1)*ekb_s(ipsang,iproj)
499 : : !
500 : : ! mjv 27 6 2008: if only_SO == 2 remove the factor of l in the operator
501 : : !
502 [ - + ]: 174118038 : if (only_SO == 2) then
503 : 0 : wt(ipsang,iproj)=4.d0*pi/ucvol*ekb_s(ipsang,iproj)
504 : : end if
505 : 174118038 : jproj(ipsang)=max(jproj(ipsang),iproj)
506 [ + - ]: 174118038 : if(iproj>0)nlang=max(nlang,ipsang)
507 : : end if
508 : : end do ! ilmn
509 : :
510 : :
511 : : ! If nlang is not 0, then some non-local part is to be applied for that type of atom.
512 [ + + ]: 88024330 : if (nlang/=0) then
513 : : ! Operate with the non-local potential on the wavefunction, in order
514 : : ! to get projected quantities. Call different routines opernl2,
515 : : ! opernl3, opernl4x which corresponds to different writings of the
516 : : ! same numerical operations. There is still optimisation left for
517 : : ! the case istwf_k/=1 (up to a factor 2 in speed).
518 : : ! call timab(74+choice,1,tsec)
519 : 43565098 : sign=1
520 [ + + + + : 12632044519 : gxa(:,:,:,:,:)=zero
+ + + + +
+ ]
521 [ + + + + : 6143043699 : dgxdt(:,:,:,:,:,:)=zero
+ + + + +
+ + + ]
522 [ + + + + : 32197922983 : dgxds(:,:,:,:,:)=zero
+ + + + +
+ ]
523 : :
524 : : ! Only the first spinorial component of vectin is taken into account first
525 [ + + ]: 43565098 : ispin=1;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index
526 [ + + ]: 43565098 : if(nloalg(1)==2) then
527 : : call opernl2(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
528 : : & ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
529 : : & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
530 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
531 : 10474 : & ntypat,ph3din,sign,vectin)
532 [ + + ]: 43554624 : else if(nloalg(1)==3) then
533 : : call opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
534 : : & ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
535 : : & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
536 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
537 : 126990 : & ntypat,ph3din,sign,vectin)
538 [ + - ]: 43427634 : else if(nloalg(1)==4) then
539 : : call opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
540 : : & ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
541 : : & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
542 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
543 : 43427634 : & ntypat,ph3din,vectin)
544 : : end if
545 : : ! This duplication of the opernl calls is needed to avoid copying
546 : : ! vectin, with a detrimental effect on speed.
547 [ + + ]: 43565098 : if (nspinor==2)then
548 : 4901725 : ispin=2
549 [ - + ]: 4901725 : if(nloalg(1)==2) then
550 : : call opernl2(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,&
551 : : & ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
552 : : & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
553 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
554 : 0 : & ntypat,ph3din,sign,vectin_s)
555 [ - + ]: 4901725 : else if(nloalg(1)==3) then
556 : : call opernl3(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,&
557 : : & ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
558 : : & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
559 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
560 : 0 : & ntypat,ph3din,sign,vectin_s)
561 [ + - ]: 4901725 : else if(nloalg(1)==4) then
562 : : call opernl4a(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,&
563 : : & ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
564 : : & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
565 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
566 : 4901725 : & ntypat,ph3din,vectin_s)
567 : : end if
568 [ + + + + : 9104999873 : dgxds(:,:,:,:,ispin)=dgxds_s(:,:,:,:)
+ + + + ]
569 [ + + + + : 1541116165 : dgxdt(:,:,:,:,:,ispin)=dgxdt_s(:,:,:,:,:)
+ + + + +
+ ]
570 [ + + + + : 3552064865 : gxa(:,:,:,:,ispin)=gxa_s(:,:,:,:)
+ + + + ]
571 : : end if
572 : :
573 : : ! Parallelism stuff
574 : 43565098 : spaceComm=mpi_enreg%comm_fft
575 : 43565098 : call xmpi_sum(dgxds,spaceComm,ierr)
576 [ + + ]: 43565098 : if (mpi_enreg%paral_spinor==1) then
577 : 19712 : spaceComm=mpi_enreg%comm_spinorfft
578 : 19712 : jspin=3-ispin
579 [ + + + + : 4908288 : gxa(:,:,:,:,ispin)=gxa(:,:,:,:,1)
+ + + + ]
580 [ + + + + : 4908288 : gxa(:,:,:,:,jspin)=zero
+ + + + ]
581 [ + + ]: 19712 : if ( ndgxdt>0)then
582 [ + + + + : 155328 : dgxdt(:,:,:,:,:,ispin)=dgxdt(:,:,:,:,:,1)
+ + + + +
+ ]
583 [ + + + + : 155328 : dgxdt(:,:,:,:,:,jspin)=zero
+ + + + +
+ ]
584 : : end if
585 : : end if
586 : :
587 : 43565098 : call xmpi_sum(gxa,spaceComm,ierr)
588 [ + + ]: 43565098 : if(ndgxdt>0) then
589 : 6998360 : call xmpi_sum(dgxdt,spaceComm,ierr)
590 : : end if
591 : :
592 : : ! XG030513 : MPIWF, at this place, one should perform the reduction
593 : : ! and spread of data of gxa, dgxdt and dgxds
594 : :
595 : :
596 : :
597 : : ! BUG FIX START
598 : : ! MS100725: First loop over spins, ilang, proj to take care of the ddk
599 : : ! decompaction PRIOR to entering the main loop. This fixes a subtle bug
600 : : ! in the calculation of the velocity operator with SOC
601 [ + + ]: 92031921 : do isp=1,nspinor
602 [ + + ]: 48466823 : ispin=isp;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index
603 : :
604 [ + + ]: 167556628 : do ia=1,nincat
605 [ + + ]: 315037219 : do ilang=1,nlang
606 : 191045689 : nproj=jproj(ilang)
607 [ + + ]: 266570396 : if(nproj/=0) then
608 : 176785943 : ilang2=(ilang*(ilang+1))/2
609 [ + + ]: 506152007 : do iproj=1,nproj
610 : :
611 : : ! The rank of the tensor gxa equals l:
612 : 329366064 : rank=ilang-1
613 : : ! jjs gives the starting address of the relevant components
614 : 329366064 : jjs=1+((ilang-1)*ilang*(ilang+1))/6
615 [ - + ]: 329366064 : if (ilang>4) then
616 : 0 : write(msg,'(a,i0)')' ilang must fall in range [1..4] but value is ',ilang
617 : 0 : ABI_BUG(msg)
618 : : end if
619 : :
620 : : ! Eventual tensorial decompaction for ddk perturbation:
621 [ + + + + ]: 506152007 : if(choice==5 .and. ilang>=2) then
622 : 4330797 : jjk=1+((ilang-2)*(ilang-1)*ilang)/6
623 : 4330797 : compact=-1
624 : : temp(:,1:(rank*(rank+1))/2)= &
625 [ + + + + ]: 28351860 : & dgxdt(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj,ispin)
626 : 4330797 : call ddkten(compact,idir,rank,temp,tmpfac)
627 : : dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)=&
628 : : & dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)&
629 [ + + + + ]: 59845908 : & +tmpfac(:,1:((rank+1)*(rank+2))/2)
630 : : end if
631 : :
632 : : ! End loop over iproj:
633 : : end do
634 : : ! End condition of existence of a reference state:
635 : : end if
636 : :
637 : : ! End loop over ilang:
638 : : end do
639 : :
640 : : ! End loop over ia:
641 : : end do
642 : :
643 : : ! End loop over isp
644 : : end do
645 : : ! BUG FIX END
646 : :
647 : :
648 : : ! Main loop over spins:
649 [ + + ]: 92031921 : do isp=1,nspinor
650 [ + + ]: 48466823 : ispin=isp;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index
651 : :
652 : : ! Perform contractions for the various tensors (d)gx?, producing the
653 : : ! contracted tensors (d)gx?fac to be passed back to opernl:
654 [ + + ]: 123991530 : do ia=1,nincat
655 [ + + ]: 315037219 : do ilang=1,nlang
656 : 191045689 : nproj=jproj(ilang)
657 [ + + ]: 266570396 : if(nproj/=0) then
658 : 176785943 : ilang2=(ilang*(ilang+1))/2
659 [ + + ]: 506152007 : do iproj=1,nproj
660 : :
661 : : ! The rank of the tensor gxa equals l:
662 : 329366064 : rank=ilang-1
663 : : ! jjs gives the starting address of the relevant components
664 : 329366064 : jjs=1+((ilang-1)*ilang*(ilang+1))/6
665 [ - + ]: 329366064 : if (ilang>4) then
666 : 0 : write(msg,'(a,i0)')' ilang must fall in range [1..4] but value is ',ilang
667 : 0 : ABI_BUG(msg)
668 : : end if
669 : :
670 : : ! Metric & spinorial contraction from gxa to gxafac. The treatment
671 : : ! is different for the scalar-relativistic and spin-orbit parts.
672 [ + + ]: 329366064 : if(ispinor==1) then
673 : : ! ------ Scalar-Relativistic ------
674 : : temp(:,1:((rank+1)*(rank+2))/2)= &
675 [ + + + + ]: 2592208547 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
676 : 264173618 : call metcon(rank,gmet,temp,tmpfac)
677 : : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
678 [ + + + + ]: 2592208547 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
679 : : else
680 : : ! ------ Spin-orbit ------
681 [ + + + + ]: 924560782 : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero
682 : : ! Contraction over spins:
683 [ + + ]: 195529950 : do ispinp=1,nspinortot
684 : : ! => Imaginary part (multiplying by i, then by the Im of amet):
685 : : temp(1,1:((rank+1)*(rank+2))/2)= &
686 [ + + ]: 703060176 : & -gxa(2,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
687 : : temp(2,1:((rank+1)*(rank+2))/2)= &
688 [ + + ]: 703060176 : & gxa(1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
689 [ + + + + ]: 1694387552 : amet_lo(:,:)=amet(2,:,:,ispin,ispinp)
690 : 130337504 : call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
691 : : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
692 : : & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
693 [ + + + + ]: 1848505520 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
694 : : ! => Real part:
695 : : temp(:,1:((rank+1)*(rank+2))/2)= &
696 [ + + + + ]: 1848505520 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
697 [ + + + + ]: 1694387552 : amet_lo(:,:)=amet(1,:,:,ispin,ispinp)
698 : 130337504 : call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
699 : : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
700 : : & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
701 [ + + + + ]: 1913697966 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
702 : : end do
703 : : end if
704 : :
705 : : ! Eventual tensorial compaction of gxafac for ddk
706 : : ! perturbation:
707 [ + + + + ]: 329366064 : if(choice==5 .and. ilang>=2) then
708 : 4330797 : jjk=1+((ilang-2)*(ilang-1)*ilang)/6
709 : : ! BUG FIX START
710 : : ! MS100725: Moved this chunk of code to a separate preliminary loop (see above),
711 : : ! to fix the k-derivative of the SOC Hamiltonian
712 : : ! compact=-1
713 : : ! temp(:,1:(rank*(rank+1))/2)= &
714 : : !& dgxdt(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj,ispin)
715 : : ! call ddkten(compact,idir,rank,temp,tmpfac)
716 : : ! dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)= &
717 : : !& dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)&
718 : : !& +tmpfac(:,1:((rank+1)*(rank+2))/2)
719 : : ! BUG FIX END
720 : 4330797 : compact=1
721 : : tmpfac(:,1:((rank+1)*(rank+2))/2)= &
722 [ + + + + ]: 59845908 : & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)
723 : 4330797 : call ddkten(compact,idir,rank,temp,tmpfac)
724 : : dgxdtfac(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj)= &
725 [ + + + + ]: 28351860 : & temp(:,1:(rank*(rank+1))/2)
726 : : end if
727 : :
728 : : ! Section for strain perturbation
729 : : ! no spin-orbit yet
730 : :
731 [ + + + + ]: 329366064 : if(choice==3 .and. signs==2) then
732 : 4843316 : istr=idir
733 [ + - ]: 4843316 : if(ispinor==1) then
734 : : ! ------ Scalar-Relativistic ------
735 : : ! jjstr is the starting address for dgxds and dgxdsfac
736 : 4843316 : jjstr=-3+((ilang+1)*(ilang+2)*(ilang+3))/6
737 : : ! diagonal enlk contribution
738 : : ! note sign change (12/05/02)
739 [ + + ]: 4843316 : if(istr>3) then
740 [ + + + + ]: 27251968 : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero
741 : : else
742 : : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
743 [ + + + + ]: 25344412 : & -gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)
744 : : end if
745 [ + + + + ]: 153764492 : dgxdsfac(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)=zero
746 : 4843316 : iterm=1
747 : : temp(:,1:((rank+1)*(rank+2))/2)= &
748 [ + + + + ]: 52596380 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
749 : 4843316 : call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac)
750 : : dgxdsfac(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)= &
751 [ + + + + ]: 153764492 : & wt(ilang,iproj)*tmpfac(:,1:((rank+3)*(rank+4))/2)
752 : 4843316 : iterm=2
753 : : temp(:,1:((rank+1)*(rank+2))/2)= &
754 [ + + + + ]: 52596380 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
755 : 4843316 : call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac)
756 : : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
757 : : & +gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
758 [ + + + + ]: 52596380 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
759 : 4843316 : iterm=3
760 : : temp(:,1:((rank+3)*(rank+4))/2)= &
761 [ + + + + ]: 153764492 : dgxds(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)
762 : 4843316 : call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac)
763 : : gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
764 : : & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
765 [ + + + + ]: 52596380 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
766 : : end if
767 : : ! end section for strain perturbation
768 : : end if
769 : :
770 : : ! Eventual metric & spinorial contraction from dgxdt to dgxdtfac.
771 : : ! either for the dynamical matrix, or for the application of the
772 : : ! gradient of the operator. The treatment is different for the
773 : : ! scalar-relativistic and spin-orbit parts.
774 : : if ((choice==2.and.signs==2).or. &
775 [ + + + + : 329366064 : & (choice==5.and.signs==2).or. &
+ + ]
776 : : & (choice==4)) then
777 [ + + ]: 30000616 : mumax=ndgxdtfac;if (choice==5) mumax=1
778 [ + + ]: 61955096 : do mu=1,mumax
779 [ + + ]: 361320544 : if(ispinor==1) then
780 : : ! ------ Scalar-Relativistic ------
781 : : temp(:,1:((rank+1)*(rank+2))/2)= &
782 [ + + + + ]: 264791350 : & dgxdt(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
783 : 27972016 : call metcon(rank,gmet,temp,tmpfac)
784 : : dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
785 [ + + + + ]: 264791350 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
786 : : else
787 : : ! ------ Spin-orbit ------
788 [ + + + + ]: 56174976 : dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero
789 : : ! Contraction over spins:
790 [ + + ]: 11947392 : do ispinp=1,nspinortot
791 : : ! => Imaginary part (multiplying by i, then by the Im of amet):
792 : : temp(1,1:((rank+1)*(rank+2))/2)= &
793 [ + + ]: 42759936 : & -dgxdt(2,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
794 : : temp(2,1:((rank+1)*(rank+2))/2)= &
795 [ + + ]: 42759936 : & dgxdt(1,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
796 [ + + + + ]: 103544064 : amet_lo(:,:)=amet(2,:,:,ispin,ispinp)
797 : 7964928 : call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
798 : : dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
799 : : & dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+&
800 [ + + + + ]: 112349952 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
801 : : ! => Real part:
802 : : temp(:,1:((rank+1)*(rank+2))/2)= &
803 [ + + + + ]: 112349952 : & dgxdt(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
804 [ + + + + ]: 103544064 : amet_lo(:,:)=amet(1,:,:,ispin,ispinp)
805 : 7964928 : call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
806 : : dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
807 : : & dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+&
808 [ + + + + ]: 116332416 : & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
809 : : end do
810 : : end if
811 : : end do
812 : : end if
813 : :
814 : :
815 : : ! ---- Accumulate the nonlocal energy.
816 [ + + ]: 1391833819 : do ii=1,ilang2
817 : 1062467755 : jj=ii-1+jjs
818 : : enlk=enlk+(gxafac(1,jj,ia,iproj)*gxa(1,jj,ia,iproj,ispin)&
819 : 1391833819 : & +gxafac(2,jj,ia,iproj)*gxa(2,jj,ia,iproj,ispin) )
820 : : end do
821 : :
822 : : ! ---- Accumulate force contributions if requested.
823 : : ! Note that the contraction of gxa with dgxdt can be done like
824 : : ! a cartesian dot product now because the symmetrically-related
825 : : ! terms are accounted for with weights in gxa.
826 [ + + + + ]: 329366064 : if ((choice==2.or.choice==23) .and. signs==1) then
827 [ + + ]: 3319319 : ishift=0;if (choice==23) ishift=6
828 : 3319319 : ia5=ia+ia3-1
829 [ + + ]: 14224555 : do ii=1,ilang2
830 : 10905236 : jj=ii-1+jjs
831 [ + + ]: 46940263 : do mu=1,3
832 : : ! (includes also factor of 2 from "2*Re[stuff]")
833 : 32715708 : indx=mu+3*(ia5-1)+ishift
834 : : enlout(indx)=enlout(indx)+two*&
835 : : & ( gxafac(1,jj,ia,iproj)*dgxdt(1,mu,jj,ia,iproj,ispin)&
836 : 43620944 : & +gxafac(2,jj,ia,iproj)*dgxdt(2,mu,jj,ia,iproj,ispin))
837 : : end do
838 : : end do
839 : : end if
840 : :
841 : : ! ---- Accumulate stress tensor contributions if requested.
842 [ + + + + ]: 329366064 : if ((choice==3.or.choice==23).and.signs==1) then
843 : : ! 1- Compute contractions involving gxa and dgxds:
844 : : ! --- Same formula for Scalar-relativistic and Spin-orbit ---
845 [ + + ]: 2269630 : if (ilang==1) then
846 [ + + ]: 4888702 : do ii=1,6
847 : : rank2(ii)=2.d0*&
848 : : & (gxafac(1,1,ia,iproj)*dgxds(1,ii,ia,iproj,isp)+ &
849 : 4888702 : & gxafac(2,1,ia,iproj)*dgxds(2,ii,ia,iproj,isp) )
850 : : end do
851 [ + + ]: 1571244 : else if (ilang==2) then
852 : : call cont13(gxafac(:,jjs:jjs+2,ia,iproj),&
853 : 867129 : & dgxds(:, 7:16,ia,iproj,isp),rank2)
854 [ + + ]: 704115 : else if (ilang==3) then
855 : : call cont24(gxafac(:,jjs:jjs+5,ia,iproj),&
856 : 634069 : & dgxds(:,17:31,ia,iproj,isp),rank2)
857 [ - + ]: 70046 : else if (ilang==4) then
858 : : call cont35(gxafac(:,jjs:jjs+9,ia,iproj),&
859 : 70046 : & dgxds(:,32:52,ia,iproj,isp),rank2)
860 : : end if
861 : : ! In all cases add rank2 term into stress tensor
862 [ + + ]: 15887410 : strsnl(:)=strsnl(:)-rank2(:)
863 : : ! 2- Compute contractions involving gxa and gxa:
864 [ + + ]: 2269630 : if(ispinor==1) then
865 : : ! 2a ------ Scalar-Relativistic ------
866 [ + + ]: 1923722 : if (ilang==2) then
867 : : strsnl(1)=strsnl(1)-wt(ilang,iproj)*2.d0*&
868 : : & (gxa(1,jjs ,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispin)+&
869 : 696849 : & gxa(2,jjs ,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispin))
870 : : strsnl(2)=strsnl(2)-wt(ilang,iproj)*2.d0*&
871 : : & (gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispin)+&
872 : 696849 : & gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispin))
873 : : strsnl(3)=strsnl(3)-wt(ilang,iproj)*2.d0*&
874 : : & (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+2,ia,iproj,ispin)+&
875 : 696849 : & gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+2,ia,iproj,ispin))
876 : : strsnl(4)=strsnl(4)-wt(ilang,iproj)*2.d0*&
877 : : & (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispin)+&
878 : 696849 : & gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispin))
879 : : strsnl(5)=strsnl(5)-wt(ilang,iproj)*2.d0*&
880 : : & (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispin)+&
881 : 696849 : & gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispin))
882 : : strsnl(6)=strsnl(6)-wt(ilang,iproj)*2.d0*&
883 : : & (gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispin)+&
884 : 696849 : & gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispin))
885 [ + + ]: 1226873 : else if (ilang==3) then
886 : 491121 : call trace2(gxa(:,jjs:jjs+5,ia,iproj,ispin),gmet,trace)
887 : 491121 : call cont22(gxa(:,jjs:jjs+5,ia,iproj,ispin),gmet,rank2)
888 [ + + ]: 3437847 : do ii=1,6
889 : : strsnl(ii)=strsnl(ii)+wt(ilang,iproj)*&
890 : : & (2.d0*(trace(1)*gxa(1,jjs+ii-1,ia,iproj,ispin)+&
891 : 3437847 : & trace(2)*gxa(2,jjs+ii-1,ia,iproj,ispin))-3.d0*rank2(ii))
892 : : end do
893 [ + + ]: 735752 : else if (ilang==4) then
894 : 37366 : call cont3(gxa(:,jjs:jjs+9,ia,iproj,ispin),gmet,rank2)
895 [ + + ]: 261562 : strsnl(:)=strsnl(:)-wt(ilang,iproj)*rank2(:)
896 : : end if
897 : : else
898 : : ! 2b ------ Spin-orbit ------
899 [ + + ]: 1037724 : do ispinp=1,nspinortot
900 [ + + ]: 1037724 : if (ilang==3) then
901 : : call cont22so(gxa(:,jjs:jjs+5,ia,iproj,ispin),&
902 : : & gxa(:,jjs:jjs+5,ia,iproj,ispinp),&
903 : 285896 : & amet(:,:,:,ispin,ispinp),rank2)
904 [ + + ]: 2001272 : strsnl(:)=strsnl(:)-wt(ilang,iproj)*3.d0*rank2(:)
905 [ + + ]: 405920 : else if (ilang==4) then
906 : : call cont33so(gxa(:,jjs:jjs+9,ia,iproj,ispin),&
907 : : & gxa(:,jjs:jjs+9,ia,iproj,ispinp),&
908 : 65360 : & gmet,amet(:,:,:,ispin,ispinp),rank2)
909 [ + + ]: 457520 : strsnl(:)=strsnl(:)-wt(ilang,iproj)*rank2(:)
910 : : end if
911 : : end do
912 : : end if
913 : : ! 3- Compute contractions involving gxa and gxa due to
914 : : ! gradients of antisymmetric tensor (amet):
915 : : ! --- Only in case of Spin-orbit ---
916 [ + + ]: 2269630 : if(ispinor==2) then
917 [ + + ]: 1037724 : do ispinp=1,nspinortot
918 : : ! Be carefull: no need to compute rank2c(:,1:3) !
919 [ + + ]: 691816 : if (ilang==2) then
920 : : rank2c(1,4)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispinp)&
921 : 340560 : & +gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispinp)
922 : : rank2c(2,4)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispinp)&
923 : 340560 : & -gxa(2,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispinp)
924 : : rank2c(1,5)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp)&
925 : 340560 : & +gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp)
926 : : rank2c(2,5)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp)&
927 : 340560 : & -gxa(2,jjs+2,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp)
928 : : rank2c(1,6)=gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp)&
929 : 340560 : & +gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp)
930 : : rank2c(2,6)=gxa(1,jjs+1,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp)&
931 : 340560 : & -gxa(2,jjs+1,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp)
932 [ + + ]: 351256 : else if (ilang==3) then
933 : : call cont22cso(gxa(:,jjs:jjs+5,ia,iproj,ispin),&
934 : : & gxa(:,jjs:jjs+5,ia,iproj,ispinp),&
935 : 285896 : & gmet,rank2c)
936 [ + - ]: 65360 : else if (ilang==4) then
937 : : call cont33cso(gxa(:,jjs:jjs+9,ia,iproj,ispin),&
938 : : & gxa(:,jjs:jjs+9,ia,iproj,ispinp),&
939 : 65360 : & gmet,rank2c)
940 : : end if
941 [ + - ]: 1037724 : if (ilang>1) then
942 [ + + ]: 2767264 : do jj=1,3
943 [ + + ]: 8993608 : do ii=4,6
944 : : strsso(ii,jj)=strsso(ii,jj)-2.d0*wt(ilang,iproj)*&
945 : : & (pauli(1,ispin,ispinp,jj)*rank2c(2,ii)+&
946 : 8301792 : & pauli(2,ispin,ispinp,jj)*rank2c(1,ii))
947 : : end do
948 : : end do
949 : : end if
950 : : end do
951 : : end if
952 : : ! Enf if (choice==3)
953 : : end if
954 : :
955 : : ! ---- Accumulate dynamical matrix contributions if requested.
956 [ + + ]: 329366064 : if (choice==4) then
957 : 976932 : ia5=ia+ia3-1
958 [ + + ]: 4010446 : do ii=1,ilang2
959 : 3033514 : jj=ii-1+jjs
960 [ + + ]: 21234598 : do mu=1,6
961 : : ! (includes also factor of 2 from "2*Re[stuff]")
962 : : enlout(mu+6*(ia5-1))=enlout(mu+6*(ia5-1))+two*&
963 : : & (gxafac(1,jj,ia,iproj)*dgxdt(1,mu+3,jj,ia,iproj,ispin)&
964 : 21234598 : & +gxafac(2,jj,ia,iproj)*dgxdt(2,mu+3,jj,ia,iproj,ispin))
965 : : end do
966 [ + + ]: 12134056 : do mu=1,3
967 : : enlout(mu+6*(ia5-1))=enlout(mu+6*(ia5-1))+two*&
968 : : & (dgxdtfac(1,mu,jj,ia,iproj)*dgxdt(1,mu,jj,ia,iproj,ispin)&
969 : 12134056 : & +dgxdtfac(2,mu,jj,ia,iproj)*dgxdt(2,mu,jj,ia,iproj,ispin))
970 : : end do
971 : : enlout(4+6*(ia5-1))=enlout(4+6*(ia5-1)) +two*&
972 : : & (dgxdtfac(1,2,jj,ia,iproj)*dgxdt(1,3,jj,ia,iproj,ispin)&
973 : 3033514 : & +dgxdtfac(2,2,jj,ia,iproj)*dgxdt(2,3,jj,ia,iproj,ispin))
974 : : enlout(5+6*(ia5-1))=enlout(5+6*(ia5-1)) +two*&
975 : : & (dgxdtfac(1,1,jj,ia,iproj)*dgxdt(1,3,jj,ia,iproj,ispin)&
976 : 3033514 : & +dgxdtfac(2,1,jj,ia,iproj)*dgxdt(2,3,jj,ia,iproj,ispin))
977 : : enlout(6+6*(ia5-1))=enlout(6+6*(ia5-1)) +two*&
978 : : & (dgxdtfac(1,1,jj,ia,iproj)*dgxdt(1,2,jj,ia,iproj,ispin)&
979 : 4010446 : & +dgxdtfac(2,1,jj,ia,iproj)*dgxdt(2,2,jj,ia,iproj,ispin))
980 : : end do
981 : : end if
982 : :
983 : : ! ---- Accumulate elastic tensor contributions if requested.
984 : :
985 [ + + ]: 506152007 : if(choice==6) then
986 : : ! XG 081121 : msg to the person who has introduced this CPP option (sorry, I did not have to time to locate who did this ...)
987 : : ! This section of ABINIT should be allowed by default to the user. I have found that on the contrary, the build
988 : : ! system defaults are such that this section is forbidden by default. You might restore this flag if at the same time,
989 : : ! you modify the build system in such a way that by default this section is included, and if the user wants, it can disable it.
990 : : ! #if defined USE_NLSTRAIN_LEGENDRE
991 : 164504 : ia5=ia+ia3-1
992 : 164504 : jjs1=((ilang)*(ilang+1)*(ilang+2))/6
993 : 164504 : jjs2=-3+((ilang+1)*(ilang+2)*(ilang+3))/6
994 : 164504 : jjs3=-9+((ilang+2)*(ilang+3)*(ilang+4))/6
995 : 164504 : jjs4=-19+((ilang+3)*(ilang+4)*(ilang+5))/6
996 : :
997 : : ! Contribution for two diagonal strains (basically, the nonlocal
998 : : ! enlk)
999 : : temp(:,1:((rank+1)*(rank+2))/2)= &
1000 [ + + + + ]: 1781360 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
1001 : 164504 : call metcon(rank,gmet,temp,tmpfac)
1002 : 164504 : e2nldd=zero
1003 [ + + ]: 703456 : do ii=1,((rank+1)*(rank+2))/2
1004 : : e2nldd=e2nldd+&
1005 : : & (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+&
1006 : 703456 : & gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii))
1007 : : end do
1008 : :
1009 : : ! Terms involving one ucvol derivative (diagonal strain only)
1010 : : ! and one derivative of nonlocal operator
1011 : : ! Loop over strain index
1012 [ + + ]: 1151528 : do istr2=1,6
1013 : :
1014 : : ! rank->rank+2
1015 : 987024 : iterm=1
1016 : : temp(:,1:((rank+1)*(rank+2))/2)= &
1017 [ + + + + ]: 10688160 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
1018 : 987024 : call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac)
1019 : 987024 : e2nl_tmp(istr2)=0.d0
1020 [ + + ]: 11084976 : do ii=1,((rank+3)*(rank+4))/2
1021 : : e2nl_tmp(istr2)=e2nl_tmp(istr2)-2.d0*&
1022 : : & (dgxds(1,jjs2-1+ii,ia,iproj,isp)*tmpfac(1,ii)+&
1023 : 11084976 : & dgxds(2,jjs2-1+ii,ia,iproj,isp)*tmpfac(2,ii))
1024 : : end do
1025 : : ! rank->rank
1026 : :
1027 : 987024 : iterm=2
1028 : : temp(:,1:((rank+1)*(rank+2))/2)= &
1029 [ + + + + ]: 10688160 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
1030 : 987024 : call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac)
1031 [ + + ]: 4385240 : do ii=1,((rank+1)*(rank+2))/2
1032 : : e2nl_tmp(istr2)=e2nl_tmp(istr2)-&
1033 : : & (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+&
1034 : 4220736 : & gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii))
1035 : : end do
1036 : : ! DEBUG
1037 : : ! This and subsequent similar debug sections evaluate the
1038 : : ! hermitial conjugate of the contraction immeditely above
1039 : : ! and the comparison was useful for development purposes.
1040 : : ! rank+2->rank
1041 : : ! iterm=3
1042 : : ! temp(:,1:((rank+3)*(rank+4))/2)= &
1043 : : ! dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin)
1044 : : ! call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac)
1045 : : ! e2nl_tmp(istr2)=0.d0
1046 : : ! do ii=1,((rank+1)*(rank+2))/2
1047 : : ! e2nl_tmp(istr2)=e2nl_tmp(istr2)-&
1048 : : ! & (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+&
1049 : : ! & gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii))
1050 : : ! end do
1051 : : ! ENDDEBUG
1052 : : end do
1053 : :
1054 : : ! Terms involving two derivatives of the nonlocal operator
1055 : : ! Loop over 2nd strain index
1056 [ + + ]: 1151528 : do istr2=1,6
1057 : : ! Loop over 1st strain index, upper triangle only
1058 [ + + ]: 6909168 : do istr1=1,6
1059 : 5922144 : iest=istr1+(3*natom+6)*(istr2-1)
1060 : :
1061 : : ! Accumulate terms corresponding to two derivatives of ucvol
1062 : : ! (simply the nonlocal energy contributin) for both indices
1063 : : ! corresponding to diagonal strains
1064 : :
1065 [ + + + + ]: 5922144 : if(istr1<=3 .and. istr2<=3) then
1066 : 1480536 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nldd
1067 : : end if
1068 : :
1069 : : ! Accumulate terms computed above from 1st derivatives
1070 : : ! when one or more indices corresponds to diagonal strain
1071 [ + + ]: 5922144 : if(istr2<=3) then
1072 : 2961072 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl_tmp(istr1)
1073 : : end if
1074 [ + + ]: 5922144 : if(istr1<=3) then
1075 : 2961072 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl_tmp(istr2)
1076 : : end if
1077 : :
1078 : : ! rank->rank+4
1079 : : call contstr21(istr2,istr1,rank,gmet,gprimd,e2nl,&
1080 : : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1081 : 5922144 : & d2gxds2(:,jjs4:jjs4-1+((rank+5)*(rank+6))/2,ia,iproj,isp))
1082 : 5922144 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1083 : :
1084 : : ! DEBUG
1085 : : ! rank+4->rank
1086 : : ! call contstr22(istr2,istr1,rank,gmet,gprimd,e2nl,&
1087 : : ! & d2gxds2(:,jjs4:jjs4-1+((rank+5)*(rank+6))/2,ia,iproj,ispin),&
1088 : : ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1089 : : ! enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1090 : : ! ENDDEBUG
1091 : :
1092 : : ! rank->rank+2
1093 : : call contstr23(istr2,istr1,rank,gmet,gprimd,e2nl,&
1094 : : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1095 : 5922144 : & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1096 : 5922144 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1097 : : ! DEBUG
1098 : :
1099 : : ! rank+2->rank
1100 : : ! call contstr24(istr2,istr1,rank,gmet,gprimd,e2nl,&
1101 : : ! & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin),&
1102 : : ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1103 : : ! enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1104 : : ! ENDDEBUG
1105 : :
1106 : : ! rank+2->rank+2
1107 [ + + ]: 5922144 : if(rank<=2) then
1108 : : call contstr25(istr2,istr1,rank,gmet,gprimd,e2nl,&
1109 : : & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp),&
1110 : 5921568 : & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1111 : : else
1112 : : call contstr25a(istr2,istr1,rank,gmet,gprimd,e2nl,&
1113 : : & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp),&
1114 : 576 : & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1115 : : end if
1116 : 5922144 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1117 : :
1118 : : ! rank->rank
1119 : : call contstr26(istr2,istr1,rank,gmet,gprimd,e2nl,&
1120 : : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1121 : 5922144 : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1122 : 6909168 : enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1123 : :
1124 : : end do !istr1
1125 : :
1126 : : ! Contributions to internal strain (one cartesian strain and one
1127 : : ! reduced-coordinate atomic displacement derivative).
1128 : 987024 : iest=7+3*(ia5-1)+(3*natom+6)*(istr2-1)
1129 : :
1130 : : ! rank->rank+3
1131 : : call contistr03(istr2,rank,gmet,gprimd,eisnl,&
1132 : : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1133 : 987024 : & d2gxdis(:,jjs3:jjs3-1+((rank+4)*(rank+5))/2,ia,iproj,isp))
1134 : : enlout(iest:iest+2)= enlout(iest:iest+2)&
1135 [ + + ]: 3948096 : & +wt(ilang,iproj)*eisnl(1:3)
1136 : :
1137 : : ! DEBUG
1138 : : ! rank+3->rank
1139 : : ! call contistr30(istr2,rank,gmet,gprimd,eisnl,&
1140 : : ! & d2gxdis(:,jjs3:jjs3-1+((rank+4)*(rank+5))/2,ia,iproj,ispin),&
1141 : : ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1142 : : ! enlout(iest:iest+2)= enlout(iest:iest+2)&
1143 : : ! & +wt(ilang,iproj)*eisnl(1:3)
1144 : : ! ENDDEBUG
1145 : :
1146 : : ! rank+1->rank+2
1147 : : call contistr12(istr2,rank,gmet,gprimd,eisnl,&
1148 : : & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,isp),&
1149 : 987024 : & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1150 : : enlout(iest:iest+2)= enlout(iest:iest+2)&
1151 [ + + ]: 3948096 : & +wt(ilang,iproj)*eisnl(1:3)
1152 : :
1153 : : ! DEBUG
1154 : : ! rank+2->rank+1
1155 : : ! call contistr21(istr2,rank,gmet,gprimd,eisnl,&
1156 : : ! & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin),&
1157 : : ! & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,ispin))
1158 : : ! enlout(iest:iest+2)= enlout(iest:iest+2)&
1159 : : ! & +wt(ilang,iproj)*eisnl(1:3)
1160 : : ! ENDDEBUG
1161 : :
1162 : : ! rank->rank+1
1163 : : call contistr01(istr2,rank,gmet,gprimd,eisnl,&
1164 : : & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1165 : 987024 : & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,isp))
1166 : : enlout(iest:iest+2)= enlout(iest:iest+2)&
1167 [ + + ]: 4112600 : & +wt(ilang,iproj)*eisnl(1:3)
1168 : : !
1169 : : ! DEBUG
1170 : : ! rank+1->rank
1171 : : ! call contistr10(istr2,rank,gmet,gprimd,eisnl,&
1172 : : ! & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,ispin),&
1173 : : ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1174 : : ! enlout(iest:iest+2)= enlout(iest:iest+2)&
1175 : : ! & +wt(ilang,iproj)*eisnl(1:3)
1176 : : ! ENDDEBUG
1177 : :
1178 : : end do !istr2
1179 : : end if !choice==6
1180 : :
1181 : : ! End loop over iproj:
1182 : : end do
1183 : : ! End condition of existence of a reference state:
1184 : : end if
1185 : :
1186 : : ! End loop over ilang:
1187 : : end do
1188 : :
1189 : : ! End loop over ia:
1190 : : end do
1191 : :
1192 : : ! Operate with the non-local potential on the projected scalars,
1193 : : ! in order to get matrix element [NOT needed if only force or stress
1194 : : ! or dynamical matrix is being computed].
1195 : :
1196 [ + + ]: 92031921 : if (signs==2) then
1197 [ + + + - ]: 47778907 : if(nloalg(2)<=0 .and. choice==2)then
1198 : : ! Prepare the phase k+q factors for the atoms between ia3 and ia4:
1199 : : ! they were not prepared previously for nloalg(2)<=0 and choice==2.
1200 : : call ph1d3d(ia3,ia4,kgout,matblk,natom,npwout,&
1201 : 2835862 : & n1,n2,n3,phkxredout,ph1d,ph3dout)
1202 : : end if
1203 : :
1204 : : ! call timab(74+choice,1,tsec)
1205 : 47778907 : sign=-1
1206 : : ! The duplication of the opernl calls has the purpose to avoid
1207 : : ! copying vectout/vectout_s
1208 [ + + + + ]: 47778907 : if(ispin==1.or.nspinor==1)then
1209 [ + + ]: 42915664 : if( nloalg(1)==2) then
1210 : : call opernl2(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1211 : : & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1212 : : & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1213 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1214 : 10394 : & ntypat,ph3dout,sign,vectout)
1215 [ + + ]: 42905270 : else if( nloalg(1)==3) then
1216 : : call opernl3(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1217 : : & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1218 : : & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1219 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1220 : 125702 : & ntypat,ph3dout,sign,vectout)
1221 [ + - ]: 42779568 : else if( nloalg(1)==4) then
1222 : : call opernl4b(choice,dgxdsfac,dgxdtfac,ffnlout,gmet,gxafac,&
1223 : : & ia3,idir,indlmn_s,ispinor,itypat,jproj,kgout,kpgout,kptout,&
1224 : : & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,&
1225 : 42779568 : & dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,ntypat,ph3dout,vectout)
1226 : : end if
1227 : : else ! if ispin == 2
1228 [ + + + + ]: 2752832490 : vectout_s(:,:)=zero
1229 [ - + ]: 4863243 : if( nloalg(1)==2) then
1230 : : call opernl2(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1231 : : & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1232 : : & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1233 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1234 : 0 : & ntypat,ph3dout,sign,vectout_s)
1235 [ - + ]: 4863243 : else if( nloalg(1)==3) then
1236 : : call opernl3(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1237 : : & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1238 : : & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1239 : : & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1240 : 0 : & ntypat,ph3dout,sign,vectout_s)
1241 [ + - ]: 4863243 : else if( nloalg(1)==4) then
1242 : : call opernl4b(choice,dgxds,dgxdtfac,ffnlout,gmet,gxafac,&
1243 : : & ia3,idir,indlmn_s,ispinor,itypat,jproj,kgout,kpgout,kptout,&
1244 : : & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,&
1245 : 4863243 : & dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,ntypat,ph3dout,vectout_s)
1246 : : end if
1247 : : vectout(1,1+npwout:2*npwout)=&
1248 [ + + ]: 920852992 : & vectout(1,1+npwout:2*npwout)+vectout_s(1,1:npwout)
1249 : : vectout(2,1+npwout:2*npwout)=&
1250 [ + + ]: 920852992 : & vectout(2,1+npwout:2*npwout)+vectout_s(2,1:npwout)
1251 : :
1252 : : end if ! end ispin if
1253 : : end if ! end signs==2 if
1254 : :
1255 : : ! End loops over spins:
1256 : : end do
1257 : :
1258 : : ! End condition of existence of a non-local part for that type of atom:
1259 : : end if
1260 : :
1261 : : ! End loop over ispinor:
1262 : : end do
1263 : :
1264 : : ! End sum on atom subset loop, over ia3:
1265 : : end do
1266 : :
1267 : : ! End atom type loop, over itypat:
1268 : 71836647 : ia1=ia2+1
1269 : : end do
1270 : :
1271 : : !De-allocate temporary space.
1272 [ - + ]: 29457648 : ABI_FREE(ekb_s)
1273 [ - + ]: 29457648 : ABI_FREE(gxa)
1274 [ - + ]: 29457648 : ABI_FREE(gxafac)
1275 [ - + ]: 29457648 : ABI_FREE(dgxds)
1276 [ - + ]: 29457648 : ABI_FREE(dgxdt)
1277 [ - + ]: 29457648 : ABI_FREE(dgxdtfac)
1278 [ - + ]: 29457648 : ABI_FREE(wt)
1279 [ - + ]: 29457648 : ABI_FREE(jproj)
1280 [ - + ]: 29457648 : ABI_FREE(temp)
1281 [ - + ]: 29457648 : ABI_FREE(tmpfac)
1282 [ - + ]: 29457648 : ABI_FREE(dgxdsfac)
1283 [ - + ]: 29457648 : ABI_FREE(indlmn_s)
1284 : : !if(choice==6) then
1285 [ - + ]: 29457648 : ABI_FREE(dgxdis)
1286 [ - + ]: 29457648 : ABI_FREE(d2gxdis)
1287 [ - + ]: 29457648 : ABI_FREE(d2gxds2)
1288 : : !end if
1289 : : !if(nspinor==2) then
1290 [ - + ]: 29457648 : ABI_FREE(dgxds_s)
1291 [ - + ]: 29457648 : ABI_FREE(dgxdt_s)
1292 [ - + ]: 29457648 : ABI_FREE(gxa_s)
1293 : : !end if
1294 : : !if(nspinor==2.and.choice==6) then
1295 [ - + ]: 29457648 : ABI_FREE(dgxdis_s)
1296 [ - + ]: 29457648 : ABI_FREE(d2gxdis_s)
1297 [ - + ]: 29457648 : ABI_FREE(d2gxds2_s)
1298 : : !end if
1299 [ + + ]: 29457648 : if (mpssoang>mpsang) then
1300 [ - + ]: 2339994 : ABI_FREE(pauli)
1301 : : end if
1302 : :
1303 : : !Restore the original content of the vectin array.
1304 : : !Note that only the first part was modified
1305 [ + + ]: 29457648 : if(istwf_k/=1) then
1306 [ - + - - : 904765 : call scalewf_nonlop(istwf_k,mpi_enreg,npwin,2,vectin)
- - - - -
+ - - -
- ]
1307 : : end if
1308 : :
1309 [ + + ]: 29457648 : if (nspinor==2) then
1310 [ - + ]: 1552174 : ABI_FREE(vectin_s)
1311 [ - + ]: 1552174 : ABI_FREE(vectout_s)
1312 : : end if
1313 : :
1314 [ + + ]: 29457648 : if (mpi_enreg%paral_spinor==1) then
1315 [ + - ]: 4928 : if (size(enlout)>0) call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr)
1316 : 4928 : call xmpi_sum(strsnl,mpi_enreg%comm_spinor,ierr)
1317 : 4928 : call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr)
1318 : 4928 : call xmpi_sum(strsso,mpi_enreg%comm_spinor,ierr)
1319 : : end if
1320 : :
1321 : : !Save non-local energy
1322 [ + + + - ]: 29457648 : if((choice==1).and.size(enlout)>0) enlout(1)=enlk ! on test v4/93 size(enlout) can be zero (PMA)
1323 : :
1324 : : !Do final manipulations to produce strain gradients for
1325 : : !stress tensor, in cartesian coordinates
1326 [ + + + + : 29564436 : if ((choice==3.or.choice==23) .and. signs==1) then
- + - + -
+ ]
1327 : : ! Convert strsnl from reduced to cartesian coordinates
1328 : 106788 : strsnl_out(:)=0.d0
1329 : 106788 : call strconv(strsnl,gprimd,strsnl_out)
1330 : 106788 : strsnl(:) = strsnl_out(:)
1331 : : ! Add diagonal part (fill up first 6 components of enlout with
1332 : : ! these gradients; elements 7,8,9 of enlout are not used)
1333 : 106788 : enlout(1)=strsnl(1)-enlk
1334 : 106788 : enlout(2)=strsnl(2)-enlk
1335 : 106788 : enlout(3)=strsnl(3)-enlk
1336 : 106788 : enlout(4)=strsnl(4)
1337 : 106788 : enlout(5)=strsnl(5)
1338 : 106788 : enlout(6)=strsnl(6)
1339 : : ! Eventually, add spin-orbit part due to gradients of
1340 : : ! antisymmetric tensor
1341 [ + + ]: 106788 : if (mpssoang>mpsang) then
1342 : 17537 : call strsocv(strsso,gprimd,strssoc)
1343 [ + + ]: 122759 : enlout(1:6)=enlout(1:6)+strssoc(:)
1344 : : end if
1345 : : end if
1346 : :
1347 : : !DEBUG
1348 : : !write(std_out,*)' nonlop_pl: exit '
1349 : : !ENDDEBUG
1350 : :
1351 : : contains
1352 : : !!***
1353 : :
1354 : : !!****f* ABINIT/trace2
1355 : : !! NAME
1356 : : !! trace2
1357 : : !!
1358 : : !! FUNCTION
1359 : : !! Sum indices to compute trace of rank 2 tensor gxa related to l=2
1360 : : !! $trace=sum_{i,j} {gxa(i,j) gmet(i,j)}$
1361 : : !!
1362 : : !! INPUTS
1363 : : !! gxa(2,6) = $sum_{G} {e^(2 \pi i G cdot t} {{f_2(|k+G|)} \over {|k+G|^2}} (k+G) cdot (k+G) C(G_{nk})}$
1364 : : !! gmet(3,3)=(symmetric) metric tensor in reciprocal space (bohr^-2)
1365 : : !!
1366 : : !! OUTPUT
1367 : : !! trace(2)=sum_{i,j} {gxa(i,j) gmet(i,j)}$ (Re and Im).
1368 : : !!
1369 : : !! NOTES
1370 : : !! Here index 6 refers to vector components
1371 : : !! of (k+G) but note tensor is symmetric=>only 6 components.
1372 : : !! The components are given in the order 11 22 33 32 31 21.
1373 : : !! The initial 2 handles the Re and Im parts.
1374 : : !!
1375 : : !! PARENTS
1376 : : !! m_nonlop_pl
1377 : : !!
1378 : : !! CHILDREN
1379 : : !!
1380 : : !! SOURCE
1381 : :
1382 : 491121 : subroutine trace2(gxa,gmet,trace)
1383 : :
1384 : : !Arguments ------------------------------------
1385 : : !arrays
1386 : : real(dp),intent(in) :: gmet(3,3),gxa(2,6)
1387 : : real(dp),intent(out) :: trace(2)
1388 : :
1389 : : !Local variables-------------------------------
1390 : : !scalars
1391 : : integer :: reim
1392 : :
1393 : : ! *************************************************************************
1394 : :
1395 : : !Write out index summation, Re and Im parts
1396 [ + + ]: 1473363 : do reim=1,2
1397 : : trace(reim)=gxa(reim,1)*gmet(1,1)+gxa(reim,2)*gmet(2,2)+&
1398 : : & gxa(reim,3)*gmet(3,3)+&
1399 : : & 2.0d0*(gxa(reim,4)*gmet(3,2)+gxa(reim,5)*gmet(3,1)+&
1400 : 1473363 : & gxa(reim,6)*gmet(2,1))
1401 : : end do
1402 : :
1403 : 491121 : end subroutine trace2
1404 : : !!***
1405 : :
1406 : : !!****f* ABINIT/strsocv
1407 : : !! NAME
1408 : : !! strsocv
1409 : : !!
1410 : : !! FUNCTION
1411 : : !! Convert from antisymmetric storage mode 3x3x3 rank3 tensor in reduced
1412 : : !! coordinates "red" to symmetric storage mode 3x3 rank2 tensor in
1413 : : !! cartesian coordinates "cart", using metric tensor "gprimd".
1414 : : !!
1415 : : !! INPUTS
1416 : : !! red(6,3)=3x3x3 tensor in antisymmetric storage mode,
1417 : : !! reduced coordinates
1418 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
1419 : : !!
1420 : : !! OUTPUT
1421 : : !! cart(6)=3x3 tensor in symmetric storage mode,
1422 : : !! cartesian coordinates
1423 : : !!
1424 : : !! NOTES
1425 : : !! This routine is used to compute spin-orbit stress tensor.
1426 : : !!
1427 : : !! red is antisymmetric in first two indices and stored
1428 : : !! as 11 22 33 32 31 21.
1429 : : !! cart is symmetric and stored as 11 22 33 32 31 21.
1430 : : !!
1431 : : !!{{\ \begin{eqnarray}
1432 : : !! cart(1,1) & = & & red(i,j,2) G(3,i) G(1,j) + red(i,j,3) G(1,i) G(2,j) \nonumber
1433 : : !! cart(2,2) & = & & red(i,j,1) G(2,i) G(3,j) + red(i,j,3) G(1,i) G(2,j) \nonumber
1434 : : !! cart(3,3) & = & & red(i,j,1) G(2,i) G(3,j) + red(i,j,2) G(3,i) G(1,j) \nonumber
1435 : : !! cart(3,2) & = & 0.5 ( & red(i,j,3) G(1,i) G(3,j) + red(i,j,2) G(2,i) G(1,j)) \nonumber
1436 : : !! cart(3,1) & = & 0.5 ( & red(i,j,3) G(3,i) G(2,j) + red(i,j,1) G(2,i) G(1,j)) \nonumber
1437 : : !! cart(2,1) & = & 0.5 ( & red(i,j,2) G(3,i) G(2,j) + red(i,j,1) G(1,i) G(3,j))
1438 : : !! \end{eqnarray} }}
1439 : : !!
1440 : : !! PARENTS
1441 : : !! m_nonlop_pl
1442 : : !!
1443 : : !! CHILDREN
1444 : : !!
1445 : : !! SOURCE
1446 : :
1447 : 17537 : subroutine strsocv(red,gprimd,cart)
1448 : :
1449 : : !Arguments ------------------------------------
1450 : : !arrays
1451 : : real(dp),intent(in) :: gprimd(3,3),red(6,3)
1452 : : real(dp),intent(out) :: cart(6)
1453 : :
1454 : : !Local variables-------------------------------
1455 : : !scalars
1456 : : integer :: ii,jj
1457 : : !arrays
1458 : : real(dp) :: work(3,3,3)
1459 : :
1460 : : ! *************************************************************************
1461 : :
1462 [ + + ]: 70148 : do ii=1,3
1463 : 52611 : work(1,1,ii)=0.d0
1464 : 52611 : work(2,2,ii)=0.d0
1465 : 52611 : work(3,3,ii)=0.d0
1466 : 52611 : work(3,2,ii)=red(4,ii) ; work(2,3,ii)=-red(4,ii)
1467 : 52611 : work(3,1,ii)=red(5,ii) ; work(1,3,ii)=-red(5,ii)
1468 : 70148 : work(2,1,ii)=red(6,ii) ; work(1,2,ii)=-red(6,ii)
1469 : : end do
1470 : :
1471 : 17537 : cart(:)=0.d0
1472 [ + + ]: 70148 : do jj=1,3
1473 [ + + ]: 227981 : do ii=1,3
1474 : : cart(1)=cart(1)+work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)&
1475 : 157833 : & +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj)
1476 : : cart(2)=cart(2)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)&
1477 : 157833 : & +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj)
1478 : : cart(3)=cart(3)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)&
1479 : 157833 : & +work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)
1480 : : cart(4)=cart(4)+work(ii,jj,3)*gprimd(1,ii)*gprimd(3,jj)&
1481 : 157833 : & +work(ii,jj,2)*gprimd(2,ii)*gprimd(1,jj)
1482 : : cart(5)=cart(5)+work(ii,jj,3)*gprimd(3,ii)*gprimd(2,jj)&
1483 : 157833 : & +work(ii,jj,1)*gprimd(2,ii)*gprimd(1,jj)
1484 : : cart(6)=cart(6)+work(ii,jj,2)*gprimd(3,ii)*gprimd(2,jj)&
1485 : 210444 : & +work(ii,jj,1)*gprimd(1,ii)*gprimd(3,jj)
1486 : : end do
1487 : : end do
1488 : 17537 : cart(4)=0.5d0*cart(4)
1489 : 17537 : cart(5)=0.5d0*cart(5)
1490 : 17537 : cart(6)=0.5d0*cart(6)
1491 : :
1492 : 17537 : end subroutine strsocv
1493 : : !!***
1494 : :
1495 : : !!****f* ABINIT/scalewf_nonlop
1496 : : !! NAME
1497 : : !! scalewf_nonlop
1498 : : !!
1499 : : !! FUNCTION
1500 : : !! At the start of nonlop (or similar routines), as well as its end,
1501 : : !! the wavefunctions, when stored with istwfk/=2,
1502 : : !! need to be scaled (by a factor of 2 or 1/2),
1503 : : !! except for the G=0 component.
1504 : : !! Only the first spinor component is to be modified.
1505 : : !!
1506 : : !! INPUTS
1507 : : !! istwf_k=storage mode of the vector
1508 : : !! mpi_enreg=information about MPI parallelization
1509 : : !! npw=number of planewaves
1510 : : !! option=1 multiply by 2
1511 : : !! =2 multiply by 1/2
1512 : : !!
1513 : : !! OUTPUT
1514 : : !! (see side effects)
1515 : : !!
1516 : : !! SIDE EFFECTS
1517 : : !! vect(2,npw)=vector that is rescaled
1518 : : !!
1519 : : !! NOTES
1520 : : !! XG030513 : MPIWF One should pay attention to the
1521 : : !! G=0 component, that will be only one one proc...
1522 : : !!
1523 : : !! PARENTS
1524 : : !! m_nonlop_pl
1525 : : !!
1526 : : !! CHILDREN
1527 : : !!
1528 : : !! SOURCE
1529 : :
1530 : 1809530 : subroutine scalewf_nonlop(istwf_k,mpi_enreg,npw,option,vect)
1531 : :
1532 : : !Arguments ------------------------------------
1533 : : !This type is defined in defs_mpi
1534 : : !scalars
1535 : : integer,intent(in) :: istwf_k,npw,option
1536 : : type(MPI_type),intent(in) :: mpi_enreg
1537 : : !arrays
1538 : : real(dp),intent(inout) :: vect(2,npw)
1539 : :
1540 : : !Local variables-------------------------------
1541 : : !scalars
1542 : : integer :: ipw
1543 : : real(dp) :: scale
1544 : : character(len=500) :: msg
1545 : :
1546 : : ! *************************************************************************
1547 : :
1548 : : DBG_ENTER("COLL")
1549 : :
1550 [ + - ]: 1809530 : if(istwf_k/=1)then
1551 : :
1552 [ - + ]: 1809530 : if(option/=1 .and. option/=2)then
1553 : : write(msg,'(a,a,a,i0)')&
1554 : 0 : 'The argument option should be 1 or 2,',ch10,&
1555 : 0 : 'however, option=',option
1556 : 0 : ABI_BUG(msg)
1557 : : end if
1558 : :
1559 : 1809530 : scale=two
1560 [ + + ]: 1809530 : if(option==2)scale=half
1561 : :
1562 : : ! Storage for the Gamma point. The component of the G=0 vector
1563 : : ! should not be scaled, and no G=0 imaginary part is allowed.
1564 [ + + ]: 1809530 : if(istwf_k==2)then
1565 [ + + ]: 993540 : if (mpi_enreg%me_g0_fft==1) then
1566 : 991160 : vect(2,1)=zero
1567 : : !$OMP PARALLEL DO
1568 [ + + ]: 632878994 : do ipw=2,npw
1569 : 631887834 : vect(1,ipw)=scale*vect(1,ipw)
1570 : 632878994 : vect(2,ipw)=scale*vect(2,ipw)
1571 : : end do
1572 : : !$OMP END PARALLEL DO
1573 : : else
1574 : : !$OMP PARALLEL DO
1575 [ + + ]: 247520 : do ipw=1,npw
1576 : 245140 : vect(1,ipw)=scale*vect(1,ipw)
1577 : 247520 : vect(2,ipw)=scale*vect(2,ipw)
1578 : : end do
1579 : : !$OMP END PARALLEL DO
1580 : : end if
1581 : : end if
1582 : :
1583 : : ! Other storage modes, for k points with time-reversal symmetry.
1584 : : ! All components should be scaled.
1585 [ + + ]: 1809530 : if(istwf_k>2)then
1586 : : !$OMP PARALLEL DO
1587 [ + + ]: 262493056 : do ipw=1,npw
1588 : 261677066 : vect(1,ipw)=scale*vect(1,ipw)
1589 : 262493056 : vect(2,ipw)=scale*vect(2,ipw)
1590 : : end do
1591 : : !$OMP END PARALLEL DO
1592 : : end if
1593 : :
1594 : : end if ! istwf_k/=1
1595 : :
1596 : : DBG_EXIT("COLL")
1597 : :
1598 : 1809530 : end subroutine scalewf_nonlop
1599 : : !!***
1600 : :
1601 : : !!****f* ABINIT/ddkten
1602 : : !! NAME
1603 : : !! ddkten
1604 : : !!
1605 : : !! FUNCTION
1606 : : !! Compact or decompact the tensors related to the ffnl(:,1,...)
1607 : : !! part of the ddk operator, taking into account the direction
1608 : : !! of the ddk perturbation.
1609 : : !!
1610 : : !! INPUTS
1611 : : !! compact= if 1, compact from tmpfac
1612 : : !! idir=direction of the ddk perturbation
1613 : : !! rank=0,1,2, or 3 = rank of tmpfac tensor, also angular momentum (=l)
1614 : : !!
1615 : : !! OUTPUT
1616 : : !! (see side effects)
1617 : : !!
1618 : : !! SIDE EFFECTS
1619 : : !! Input/Output:
1620 : : !! temp(2,(rank*(rank+1))/2)=compacted tensor
1621 : : !! for l=1, just a scalar
1622 : : !! for l=2, a vector
1623 : : !! tmpfac(2,(rank+1)*(rank+2)/2)=decompacted tensor
1624 : : !! for l=1, a vector
1625 : : !! for l=2, a symmetric matrix, stored as
1626 : : !! (1 . .)
1627 : : !! (6 2 .)
1628 : : !! (5 4 3)
1629 : : !!
1630 : : !! NOTES
1631 : : !! For l=0, there is no contribution.
1632 : : !!
1633 : : !! PARENTS
1634 : : !! m_nonlop_pl
1635 : : !!
1636 : : !! CHILDREN
1637 : : !!
1638 : : !! SOURCE
1639 : :
1640 : 8661594 : subroutine ddkten(compact,idir,rank,temp,tmpfac)
1641 : :
1642 : : !Arguments ------------------------------------
1643 : : !scalars
1644 : : integer,intent(in) :: compact,idir,rank
1645 : : !arrays
1646 : : real(dp),intent(inout) :: temp(2,(rank*(rank+1))/2)
1647 : : real(dp),intent(inout) :: tmpfac(2,((rank+1)*(rank+2))/2)
1648 : :
1649 : : !Local variables-------------------------------
1650 : : !scalars
1651 : : character(len=500) :: msg
1652 : :
1653 : : ! *************************************************************************
1654 : :
1655 [ - + ]: 8661594 : if(rank/=1 .and. rank/=2 .and. rank/=3)then
1656 : : write(msg, '(a,i10,a,a,a)' )&
1657 : 0 : 'Input rank=',rank,' not allowed.',ch10,&
1658 : 0 : 'Possible values are 1,2,3 only.'
1659 : 0 : ABI_BUG(msg)
1660 : : end if
1661 : :
1662 : : !Take care of p angular momentum
1663 [ + + ]: 8661594 : if(rank==1)then
1664 : :
1665 : : ! Compaction tmpfac -> temp
1666 [ + + ]: 4995510 : if(compact==1)then
1667 [ + + ]: 7493265 : temp(:,1)=tmpfac(:,idir)
1668 : :
1669 : : ! Decompaction temp -> tmpfac
1670 : : else
1671 [ + + + + ]: 24977550 : tmpfac(:,1:3)=0.0d0
1672 [ + + ]: 7493265 : tmpfac(:,idir)=temp(:,1)
1673 : : end if
1674 : :
1675 : : ! Take care of d angular momentum
1676 : : ! rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1677 : :
1678 [ + + ]: 3666084 : else if(rank==2)then
1679 : :
1680 : : ! Compaction tmpfac -> temp
1681 [ + + ]: 3659324 : if(compact==1)then
1682 [ + + ]: 1829662 : if(idir==1)then
1683 : : ! Count the number of non-zero derivatives with respect to k(idir)
1684 : : ! The factor of 2 on the diagonal comes from the derivative with
1685 : : ! respect to the first K then to the second K
1686 [ + + + + : 7451982 : temp(:,1)=2.0d0*tmpfac(:,1); temp(:,2)=tmpfac(:,6); temp(:,3)=tmpfac(:,5)
+ + ]
1687 [ + + ]: 1001664 : else if(idir==2)then
1688 [ + + + + : 4527936 : temp(:,2)=2.0d0*tmpfac(:,2); temp(:,1)=tmpfac(:,6); temp(:,3)=tmpfac(:,4)
+ + ]
1689 [ + - ]: 498560 : else if(idir==3)then
1690 [ + + + + : 4487040 : temp(:,3)=2.0d0*tmpfac(:,3); temp(:,1)=tmpfac(:,5); temp(:,2)=tmpfac(:,4)
+ + ]
1691 : : end if
1692 : : ! Decompaction temp -> tmpfac
1693 : : else
1694 [ + + + + ]: 34763578 : tmpfac(:,1:6)=0.0d0
1695 [ + + ]: 5488986 : tmpfac(:,idir)=2.0d0*temp(:,idir)
1696 [ + + ]: 1829662 : if(idir==1)then
1697 [ + + + + ]: 4967988 : tmpfac(:,5)=temp(:,3); tmpfac(:,6)=temp(:,2)
1698 [ + + ]: 1001664 : else if(idir==2)then
1699 [ + + + + ]: 3018624 : tmpfac(:,4)=temp(:,3); tmpfac(:,6)=temp(:,1)
1700 [ + - ]: 498560 : else if(idir==3)then
1701 [ + + + + ]: 2991360 : tmpfac(:,4)=temp(:,2); tmpfac(:,5)=temp(:,1)
1702 : : end if
1703 : : end if
1704 : :
1705 : : ! Take care of f angular momentum
1706 [ + - ]: 6760 : else if(rank==3)then
1707 : : ! rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
1708 : : ! rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1709 : :
1710 : : ! Compaction tmpfac -> temp
1711 [ + + ]: 6760 : if(compact==1)then
1712 [ + - ]: 3380 : if(idir==1)then
1713 : : ! Count the number of non-zero derivatives with respect to k(idir)
1714 [ + + ]: 10140 : temp(:,1)=3.0d0*tmpfac(:,1)
1715 [ + + + + ]: 33800 : temp(:,2:4)=tmpfac(:,2:4)
1716 [ + + + + ]: 23660 : temp(:,5:6)=2.0d0*tmpfac(:,5:6)
1717 [ # # ]: 0 : else if(idir==2)then
1718 [ # # ]: 0 : temp(:,6)=2.0d0*tmpfac(:,2)
1719 [ # # ]: 0 : temp(:,4)=2.0d0*tmpfac(:,9)
1720 [ # # ]: 0 : temp(:,5)=tmpfac(:,4)
1721 [ # # ]: 0 : temp(:,1)=tmpfac(:,6)
1722 [ # # ]: 0 : temp(:,3)=tmpfac(:,8)
1723 [ # # ]: 0 : temp(:,2)=3.0d0*tmpfac(:,7)
1724 [ # # ]: 0 : else if(idir==3)then
1725 [ # # ]: 0 : temp(:,3)=3.0d0*tmpfac(:,10)
1726 [ # # ]: 0 : temp(:,5)=2.0d0*tmpfac(:,3)
1727 [ # # ]: 0 : temp(:,4)=2.0d0*tmpfac(:,8)
1728 [ # # ]: 0 : temp(:,6)=tmpfac(:,4)
1729 [ # # ]: 0 : temp(:,1)=tmpfac(:,5)
1730 [ # # ]: 0 : temp(:,2)=tmpfac(:,9)
1731 : : end if
1732 : : ! Decompaction temp -> tmpfac
1733 : : else
1734 [ + + + + ]: 104780 : tmpfac(:,1:10)=0.0d0
1735 [ + - ]: 3380 : if(idir==1)then
1736 [ + + ]: 10140 : tmpfac(:,1)=3.0d0*temp(:,1)
1737 [ + + + + ]: 33800 : tmpfac(:,2:4)=temp(:,2:4)
1738 [ + + + + ]: 23660 : tmpfac(:,5:6)=2.0d0*temp(:,5:6)
1739 [ # # ]: 0 : else if(idir==2)then
1740 [ # # ]: 0 : tmpfac(:,2)=2.0d0*temp(:,6)
1741 [ # # ]: 0 : tmpfac(:,9)=2.0d0*temp(:,4)
1742 [ # # ]: 0 : tmpfac(:,4)=temp(:,5)
1743 [ # # ]: 0 : tmpfac(:,6)=temp(:,1)
1744 [ # # ]: 0 : tmpfac(:,8)=temp(:,3)
1745 [ # # ]: 0 : tmpfac(:,7)=3.0d0*temp(:,2)
1746 [ # # ]: 0 : else if(idir==3)then
1747 [ # # ]: 0 : tmpfac(:,10)=3.0d0*temp(:,3)
1748 [ # # ]: 0 : tmpfac(:,3)=2.0d0*temp(:,5)
1749 [ # # ]: 0 : tmpfac(:,8)=2.0d0*temp(:,4)
1750 [ # # ]: 0 : tmpfac(:,4)=temp(:,6)
1751 [ # # ]: 0 : tmpfac(:,5)=temp(:,1)
1752 [ # # ]: 0 : tmpfac(:,9)=temp(:,2)
1753 : : end if
1754 : : end if
1755 : :
1756 : : end if
1757 : :
1758 : 8661594 : end subroutine ddkten
1759 : : !!***
1760 : :
1761 : : end subroutine nonlop_pl
1762 : : !!***
1763 : :
1764 : : end module m_nonlop_pl
1765 : : !!***
|
