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