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 : !!***
|
