Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfpt_nstwf
2 : : !! NAME
3 : : !! m_dfpt_nstwf
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group ()
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 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : ! nvtx related macro definition
22 : : #include "nvtx_macros.h"
23 : :
24 : : module m_dfpt_nstwf
25 : :
26 : : use defs_basis
27 : : use m_xmpi
28 : : use m_errors
29 : : use m_abicore
30 : : use m_wfk
31 : : use m_hamiltonian
32 : : use m_cgtools
33 : : use m_nctk
34 : : use m_dtset
35 : : use m_dtfil
36 : : use m_gputk
37 : : use m_abi_linalg
38 : :
39 : : use defs_datatypes, only : pseudopotential_type
40 : : use defs_abitypes, only : MPI_type
41 : : use m_time, only : timab
42 : : use m_io_tools, only : file_exists
43 : : use m_fourier_interpol, only : transgrid
44 : : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, gemm_nonlop_use_gemm
45 : : use m_geometry, only : stresssym
46 : : use m_dynmat, only : dfpt_sygra
47 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle, proc_distrb_band, proc_distrb_nband
48 : : use m_hdr, only : hdr_skip
49 : : use m_occ, only : occeig
50 : : use m_pawang, only : pawang_type
51 : : use m_pawrad, only : pawrad_type
52 : : use m_pawtab, only : pawtab_type
53 : : use m_paw_an, only : paw_an_type, paw_an_reset_flags
54 : : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
55 : : use m_pawfgrtab,only : pawfgrtab_type
56 : : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, pawrhoij_nullify, &
57 : : pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim
58 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_copy, pawcprj_output
59 : : use m_pawdij, only : pawdijfr
60 : : use m_pawfgr, only : pawfgr_type
61 : : use m_paw_mkrho,only : pawmkrho
62 : : use m_paw_nhat, only : pawnhatfr
63 : : use m_paw_dfpt, only : pawdfptenergy
64 : : use m_kg, only : mkkin, kpgstr, mkkpg
65 : : use m_fft, only : fftpac
66 : : use m_spacepar, only : hartrestr, symrhg
67 : : use m_initylmg, only : initylmg
68 : : use m_mkffnl, only : mkffnl
69 : : use m_getgh1c, only : getgh1c, getdc1
70 : : use m_dfpt_mkrho, only : dfpt_accrho
71 : : use m_atm2fft, only : dfpt_atm2fft
72 : : use m_mkcore, only : dfpt_mkcore
73 : : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll
74 : : use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
75 : : use m_mklocl, only : dfpt_vlocal, vlocalstr
76 : : use m_cgprj, only : getcprj
77 : :
78 : : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
79 : :
80 : : #if defined(HAVE_GPU_MARKERS)
81 : : use m_nvtx_data
82 : : #endif
83 : :
84 : : implicit none
85 : :
86 : : private
87 : : !!***
88 : :
89 : : public :: dfpt_nstpaw
90 : : public :: dfpt_nstwf
91 : : public :: gaugetransfo
92 : : !!***
93 : :
94 : : contains
95 : : !!***
96 : :
97 : : !!****f* ABINIT/dfpt_nstpaw
98 : : !! NAME
99 : : !! dfpt_nstpaw
100 : : !!
101 : : !! FUNCTION
102 : : !! Initially designed for PAW approach, but works also for NCPP.
103 : : !! This routine compute the non-stationary expression for the
104 : : !! second derivative of the total energy, for a whole row of
105 : : !! mixed derivatives (including diagonal terms contributing
106 : : !! to non-stationary 2nd-order total energy).
107 : : !! Compared with NC-pseudopotentials, PAW contributions include:
108 : : !! - changes of the overlap between 0-order wave-functions,
109 : : !! - on-site contributions.
110 : : !!
111 : : !! COPYRIGHT
112 : : !! Copyright (C) 2010-2026 ABINIT group (MT, AM)
113 : : !! This file is distributed under the terms of the
114 : : !! GNU General Public License, see ~abinit/COPYING
115 : : !! or http://www.gnu.org/copyleft/gpl.txt .
116 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
117 : : !!
118 : : !! INPUTS
119 : : !! cg (2,mpw *nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions at k
120 : : !! cgq(2,mpw1*nspinor*mband_mem*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
121 : : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
122 : : !! cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
123 : : !! cprj(natom,nspinor*mband_mem*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors
124 : : !! cprjq(natom,nspinor*mband_mem*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors
125 : : !! docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
126 : : !! doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
127 : : !! dtfil <type(datafiles_type)>=variables related to files
128 : : !! dtset <type(dataset_type)>=all input variables for this dataset
129 : : !! eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
130 : : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
131 : : !! eigen1(2*mband*mband*nkpt_rbz*nsppol)=1st-order eigenvalues at k,q (hartree)
132 : : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
133 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
134 : : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
135 : : !! idir=direction of the perturbation
136 : : !! indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points
137 : : !! indsy1(4,nsym1,natom)=indirect indexing array for atom labels
138 : : !! ipert=type of the perturbation
139 : : !! irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for RF symmetries
140 : : !! istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
141 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
142 : : !! kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
143 : : !! kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ
144 : : !! kxc(nfftf,nkxc)=exchange and correlation kernel
145 : : !! mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90)
146 : : !! mkmem =number of k points treated by this node.
147 : : !! mkqmem =number of k+q points treated by this node (GS data).
148 : : !! mk1mem =number of k points treated by this node (RF data)
149 : : !! mpert =maximum number of ipert
150 : : !! mpi_enreg=information about MPI parallelization
151 : : !! mpw=maximum dimensioned size of npw or wfs at k
152 : : !! mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
153 : : !! nattyp(ntypat)= # atoms of each type.
154 : : !! nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
155 : : !! mband_mem=number of bands per processor
156 : : !! ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
157 : : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid
158 : : !! ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid
159 : : !! nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
160 : : !! nhat1(cplex*nfftf,nspden*usepaw)=1st-order compensation charge density (PAW)
161 : : !! nkpt_rbz=number of k points in the reduced BZ for this perturbation
162 : : !! nkxc=second dimension of the kxc array
163 : : !! npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
164 : : !! npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
165 : : !! nspden=number of spin-density components
166 : : !! nspinor=number of spinorial components of the wavefunctions
167 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
168 : : !! nsym1=number of symmetry elements in space group consistent with i perturbation
169 : : !! n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, cplex*nfftf
170 : : !! occkq(mband*nkpt_rbz*nsppol)=occupation number for each band at each k+q point of the reduced BZ
171 : : !! occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced BZ
172 : : !! paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS
173 : : !! paw_an1(natom) <type(paw_an_type)>=1st-order paw arrays given on angular mesh for the perturbation (j1)
174 : : !! paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
175 : : !! paw_ij1(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels
176 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
177 : : !! pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
178 : : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
179 : : !! pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
180 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
181 : : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
182 : : !! pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
183 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
184 : : !! phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic transl. phases, for RF symmetries
185 : : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
186 : : !! ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information for the "fine" grid
187 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
188 : : !! rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
189 : : !! rhor1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3.
190 : : !! rmet(3,3)=real space metric (bohr**2)
191 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
192 : : !! symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations
193 : : !! symrc1(3,3,nsym1)=symmetry operations in reciprocal space
194 : : !! symrl1(3,3,nsym1)=symmetry operations in real space in terms
195 : : !! ucvol=unit cell volume in bohr**3.
196 : : !! usecprj= 1 if cprj, cprjq arrays are stored in memory
197 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
198 : : !! usevxctau=1 if if XC functional depends on kinetic energy density
199 : : !! usexcnhat= -PAW only- flag controlling use of compensation density in Vxc
200 : : !! useylmgr1= 1 if ylmgr1 array is allocated
201 : : !! vectornd(with_vectornd*nfftf,3)=nuclear dipole moment vector potential
202 : : !! vhartr1(cplex*nfft)=1-order Hartree potential
203 : : !! vpsp1(cplex*nfftf)=first-order derivative of the ionic potential
204 : : !! vtrial(nfftf,nspden)=GS potential (Hartree).
205 : : !! vtrial1(cplex*nfftf,nspden)= RF 1st-order potential (Hartree).
206 : : !! vxc(nfftf,nspden)=XC GS potential
207 : : !! vxctau(nfftf,nspden,4*usevxctau)=derivative of e_xc with respect to kinetic energy density, for mGGA
208 : : !! with_vectornd = 1 if vectornd allocated
209 : : !! wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
210 : : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
211 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
212 : : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
213 : : !! ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
214 : : !! ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q
215 : : !!
216 : : !! OUTPUT
217 : : !! blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
218 : : !! d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
219 : : !! d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
220 : : !! d2ovl(2,3,mpert,3,mpert*usepaw)=overlap contributions to the 2DTEs (PAW only)
221 : : !! eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
222 : : !! PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
223 : : !!
224 : : !! NOTES
225 : : !! We perform here the computation of
226 : : !! delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
227 : : !! see PRB 78, 035105 (2008), Eq. (42) [[cite:Audouze2008]]
228 : : !!
229 : : !! SOURCE
230 : :
231 : 1083 : subroutine dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,&
232 : 722 : & eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,&
233 : 361 : & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,&
234 : 1083 : & mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,mband_mem_rbz,ncpgr,nfftf,ngfftf,nhat,nhat1,&
235 : 722 : & nkpt_rbz,nkxc,npwarr,npwar1,nspden,nspinor,nsppol,nsym1,n3xccc,occkq,occ_rbz,&
236 [ - + - + ]: 361 : & paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,&
237 [ - + ]: 361 : & pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,&
238 [ - + ]: 361 : & ucvol,usecprj,usepaw,usevxctau,usexcnhat,useylmgr1,vectornd,vhartr1,vpsp1,vtrial,vtrial1,vxc,vxctau,&
239 [ - + - + : 361 : & with_vectornd,wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr1)
- + ]
240 : :
241 : : !Arguments -------------------------------
242 : : !scalars
243 : : integer,intent(in) :: cplex,idir,ipert,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpw,mpw1
244 : : integer,intent(in) :: ncpgr,nfftf,nkpt_rbz,nkxc,nspden,nspinor,nsppol,nsym1
245 : : integer,intent(in) :: n3xccc,usecprj,usepaw,usevxctau,usexcnhat,useylmgr1
246 : : integer,intent(in) :: mband_mem_rbz,with_vectornd
247 : : real(dp),intent(in) :: gsqcut,ucvol
248 : : real(dp),intent(out) :: eovl1
249 : : type(datafiles_type),intent(in) :: dtfil
250 : : type(dataset_type),intent(in) :: dtset
251 : : type(MPI_type),intent(in) :: mpi_enreg
252 : : type(pawang_type),intent(in) :: pawang,pawang1
253 : : type(pawfgr_type),intent(in) :: pawfgr
254 : : type(pseudopotential_type),intent(in) :: psps
255 : : !arrays
256 : : integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol)
257 : : integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom)
258 : : integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))
259 : : integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
260 : : integer,intent(in) :: ngfftf(18),npwarr(nkpt_rbz),npwar1(nkpt_rbz)
261 : : integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1)
262 : : integer,intent(inout) :: blkflg(3,mpert,3,mpert)
263 : : real(dp),intent(in),target :: cg(2,mpw*nspinor*mband_mem_rbz*mkmem*nsppol)
264 : : real(dp),intent(in),target :: cgq(2,mpw1*nspinor*mband_mem_rbz*mkqmem*nsppol)
265 : : real(dp),intent(in),target :: cg1(2,mpw1*nspinor*mband_mem_rbz*mk1mem*nsppol)
266 : : real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*nsppol)
267 : : real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*nsppol)
268 : : real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*nsppol),eigen0(dtset%mband*nkpt_rbz*nsppol)
269 : : real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol)
270 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt_rbz(3,nkpt_rbz)
271 : : real(dp),intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden),nhat1(cplex*nfftf,nspden*usepaw)
272 : : real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol)
273 : : real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*nsppol)
274 : : real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))
275 : : real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
276 : : real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),rhog(2,nfftf)
277 : : real(dp),intent(in) :: rhor(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden),rmet(3,3),rprimd(3,3)
278 : : real(dp),intent(in) :: tnons1(3,nsym1),vhartr1(cplex*nfftf),vtrial1(cplex*nfftf,nspden),vxc(nfftf,nspden)
279 : : real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom)
280 : : real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
281 : : real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
282 : : real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
283 : : real(dp),target,intent(in) :: vpsp1(cplex*nfftf),vtrial(nfftf,nspden),xccc3d1(cplex*n3xccc)
284 : : real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert)
285 : : real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw)
286 : : real(dp),intent(inout) :: vectornd(with_vectornd*nfftf,3)
287 : : real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*usevxctau)
288 : : type(pawcprj_type),intent(in) :: cprj(dtset%natom,nspinor*mband_mem_rbz*mkmem*nsppol*usecprj)
289 : : type(pawcprj_type),intent(in) :: cprjq(dtset%natom,nspinor*mband_mem_rbz*mkqmem*nsppol*usecprj)
290 : : type(paw_an_type),intent(in) :: paw_an(:)
291 : : type(paw_an_type),intent(inout) :: paw_an1(:)
292 : : type(paw_ij_type),intent(in) :: paw_ij(:)
293 : : type(paw_ij_type),intent(inout) :: paw_ij1(:)
294 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(:)
295 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*usepaw)
296 : : type(pawrhoij_type),intent(in) :: pawrhoij(:)
297 : : type(pawrhoij_type),intent(in) :: pawrhoij1(:)
298 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw)
299 : :
300 : : !Local variables-------------------------------
301 : : !scalars
302 : : integer,parameter :: tim_fourwf=18,tim_getgh1c=2,tim_projbd=3,formeig1=1
303 : : integer :: bd2tot_index,bdtot_index,berryopt,bufsz,choice,cpopt,cplex_rhoij,ddkcase
304 : : integer :: dimffnl,dimffnl1,dimffnl1_idir1,dimylmgr1,g0term
305 : : integer :: ia,iatom,iband,ibg,ibgq,ibg1,icg,icgq,icg1,ider,idir0,idir1,idir_cprj,ipw
306 : : integer :: ierr,ii,ikg,ikg1,ikpt,ikpt_me,ilmn,iorder_cprj,ipert1
307 : : integer :: ispden,isppol,istwf_k,istr,istr1,itypat,jband,jj,kdir1,kpert1,master,mcgq,mcprjq
308 : : integer :: mdir1,me,mpert1,my_natom,my_comm_atom,my_nsppol,nband_k,nband_kocc,need_ylmgr1,my_nproc_band
309 : : integer :: nfftot,nkpg,nkpg1,nkpt_me,npw_,npw_k,npw1_k,nspden_rhoij
310 : : integer :: nvh1,nvxc1,nzlmopt_ipert,nzlmopt_ipert1,optlocal,optnl
311 : : integer :: option,opt_gvnlx1,qphase_rhoij,sij_opt,spaceworld,usevnl,wfcorr,ik_ddk
312 : : integer :: nband_me, iband_me, jband_me, iband_
313 : : integer :: do_scprod, do_bcast
314 : : integer :: startband, endband
315 : : integer :: ndat,idat,ispinor
316 : : integer :: gpu_option
317 : : real(dp) :: arg,doti,dotr,dot1i,dot1r,dot2i,dot2r,dot3i,dot3r,elfd_fact,invocc,lambda,wtk_k
318 : : logical :: force_recompute,has_dcwf,has_dcwf2,has_drho,has_ddk_file,has_vectornd
319 : : logical :: is_metal,is_metal_or_qne0,need_ddk_file,need_pawij10
320 : : logical :: need_wfk,need_wf1,nmxc,paral_atom,qne0,t_exist,use_ompgpu,with_vxctau
321 : : character(len=500) :: msg
322 : : character(len=fnlen) :: fiwfddk(3)
323 : : complex(dp), parameter :: cminusone = (-1._dp,0._dp)
324 : 361 : type(gs_hamiltonian_type) :: gs_hamkq
325 : 361 : type(rf_hamiltonian_type) :: rf_hamkq
326 : 361 : type(MPI_type) :: mpi_enreg_seq
327 : : !arrays
328 : : integer :: ddkfil(3),my_spintab(2),nband_tmp(1),npwar1_tmp(1)
329 : 361 : integer,allocatable,target :: bands_treated_now_ndat(:,:),band_procs(:)
330 : 361 : integer,pointer :: bands_treated_now(:)
331 : 361 : integer,allocatable :: jpert1(:),jdir1(:),kg1_k(:,:),kg_k(:,:)
332 : 361 : integer,pointer :: my_atmtab(:)
333 : : real(dp) :: dum1(1,1),dum2(1,1),dum3(1,1),epawnst(2),kpoint(3),kpq(3)
334 : 361 : real(dp),allocatable :: vdotr(:),vdoti(:),vdot1r(:),vdot1i(:),vdot2r(:),vdot2i(:)
335 : : real(dp) :: sumelfd(2),symfact(3),tsec(2),ylmgr_dum(1,3,1)
336 : 361 : real(dp),allocatable :: buffer(:),ch1c(:,:,:,:),cs1c(:,:,:,:)
337 : 361 : real(dp),allocatable,target :: ch1c_tmp(:,:)
338 : 361 : real(dp),allocatable :: cs1c_tmp(:,:)
339 : 361 : real(dp),allocatable :: dcwavef(:,:)
340 : 361 : real(dp), contiguous, pointer :: cwave0(:,:),cwavef(:,:)!,dcwavef(:,:)
341 : 361 : real(dp),allocatable :: cg_ddk(:,:,:)
342 : 361 : real(dp),allocatable :: doccde_k(:),doccde_kq(:)
343 : 361 : real(dp),allocatable :: dnhat1(:,:),drhoaug1(:,:,:,:)
344 : 361 : real(dp),allocatable :: drhor1(:,:),drho1wfg(:,:),drho1wfr(:,:,:)
345 : 722 : real(dp),allocatable :: d2nl_elfd(:,:),dkinpw(:)
346 : 361 : real(dp),allocatable :: d2nl_k(:,:),d2ovl_drho(:,:,:,:,:),d2ovl_k(:,:)
347 : 361 : real(dp),allocatable :: eig_k(:),eig_kq(:),eig1_k(:)
348 : 361 : real(dp),allocatable,target :: e1kbfr_spin(:,:,:,:,:,:),ffnlk(:,:,:,:),ffnl1(:,:,:,:)
349 : 361 : real(dp),allocatable,target :: gh1(:,:),gs1(:,:),gvnlx1(:,:),gvnlx2(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:)
350 : 361 : real(dp),allocatable,target :: gvnlx1_tmp(:,:)
351 : 361 : real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:),rhotmp(:,:),rocceig(:,:)
352 : 361 : real(dp),allocatable :: vectornd_pac(:,:,:,:,:),vectornd_pac_idir(:,:,:,:),vlocal(:,:,:,:),vtrial_(:,:)
353 : 361 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
354 : 361 : real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:),vtmp1(:,:),vxc10(:,:)
355 : 361 : real(dp),allocatable,target :: work(:,:,:),e1kb_work(:,:,:,:)
356 : 361 : real(dp),pointer :: e1kbfr(:,:,:,:,:),e1kb_ptr(:,:,:,:)
357 : 361 : real(dp), contiguous, pointer :: ffnl1_idir1(:,:,:,:)
358 : 361 : real(dp),pointer :: vhartr01(:),vpsp1_idir1(:),xccc3d1_idir1(:)
359 : 361 : type(pawcprj_type),allocatable :: dcwaveprj(:,:)
360 : 361 : type(pawcprj_type),allocatable,target :: cwaveprj0(:,:)
361 : 361 : type(pawcprj_type),pointer :: cwaveprj0_idir1(:,:)
362 : 361 : type(paw_ij_type),allocatable :: paw_ij10(:,:)
363 : 361 : type(pawrhoij_type),target,allocatable :: pawdrhoij1(:,:)
364 : 361 : type(pawrhoij_type),pointer :: pawdrhoij1_unsym(:,:)
365 [ + + + + : 35378 : type(wfk_t) :: ddks(3)
+ + + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + - ]
366 : :
367 : : ! *********************************************************************
368 : :
369 : : DBG_ENTER("COLL")
370 : :
371 : : ABI_NVTX_START_RANGE(NVTX_DFPT_NSTPAW)
372 : :
373 : : !Keep track of total time spent in dfpt_nstpaw
374 : 361 : call timab(566,1,tsec)
375 : :
376 : : !Not valid for PrintBandByBand
377 [ - + ]: 361 : if (dtset%prtbbb/=0) then
378 : 0 : ABI_BUG('not yet valid for prtbbb/=0!')
379 : : end if
380 : :
381 : : !NCPP restrictions
382 [ - + ]: 361 : if (usepaw==0) then
383 : : ! cprj cannot be used
384 [ # # ]: 0 : if (usecprj/=0) then
385 : 0 : ABI_BUG('NCPP: usecprj should be 0!')
386 : : end if
387 : : ! d2ovl cannot be used
388 [ # # # # ]: 0 : if (size(d2ovl)/=0) then
389 : 0 : ABI_BUG('NCPP: d2ovl should not be allocated!')
390 : : end if
391 : : end if
392 : :
393 : : !PAW restrictions
394 [ + - ]: 361 : if (usepaw==1) then
395 : : ! Test on FFT grid sizes
396 [ - + ]: 361 : if (pawfgr%nfft/=nfftf) then
397 : 0 : ABI_BUG('PAW: wrong values for nfft, nfftf!')
398 : : end if
399 : : ! Test gradients of cprj
400 [ + + - + ]: 361 : if (ipert<=dtset%natom.and.ncpgr/=3) then
401 : 0 : ABI_BUG('PAW: wrong value of ncpgr for ipert<=natom!')
402 : : end if
403 [ + + + + : 361 : if (ipert==dtset%natom+1.and.ncpgr/=1.and.dtset%orbmag==0) then
- + ]
404 : 0 : ABI_BUG('PAW: wrong value of ncpgr for ipert=natom+1!')
405 : : end if
406 [ + + - + ]: 361 : if (ipert==dtset%natom+2.and.ncpgr/=3) then
407 : 0 : ABI_BUG('PAW: wrong value of ncpgr for ipert=natom+2!')
408 : : end if
409 [ + + + + : 361 : if ((ipert==dtset%natom+3.or.ipert==dtset%natom+4).and.ncpgr/=1) then
- + ]
410 : 0 : ABI_BUG('PAW: wrong value of ncpgr for ipert=natom+3 or 4!')
411 : : end if
412 : : ! Test on availability of DijHartree and XC on-site potentials
413 [ + + + + ]: 361 : if (mpi_enreg%my_natom>0.and.ipert/=dtset%natom+1) then
414 [ + - - + ]: 218 : if (paw_ij1(1)%has_dijhartree==0.or.paw_an1(1)%has_vxc==0) then
415 : 0 : msg='PAW: paw_ij1%dijhartree and paw=_an1%vxc1 should be allocated !'
416 : : end if
417 : : end if
418 : : end if
419 : :
420 : : !Set up parallelism
421 : 361 : gpu_option=dtset%gpu_option
422 : 361 : ndat=dtset%bandpp
423 : 361 : master=0;me=mpi_enreg%me_kpt
424 : 361 : spaceworld=mpi_enreg%comm_cell
425 : 361 : paral_atom=(mpi_enreg%my_natom/=dtset%natom)
426 : 361 : my_comm_atom=mpi_enreg%comm_atom
427 : 361 : my_natom=mpi_enreg%my_natom
428 : 361 : my_nproc_band=mpi_enreg%nproc_band
429 : 361 : my_atmtab=>mpi_enreg%my_atmtab
430 [ + + ]: 1083 : my_spintab=mpi_enreg%my_isppoltab
431 [ + + + + ]: 1083 : my_nsppol=count(my_spintab==1)
432 : 361 : use_ompgpu=(gpu_option==ABI_GPU_OPENMP)
433 : :
434 : : !Fake MPI data to be used in sequential calls to parallel routines
435 : 361 : call initmpi_seq(mpi_enreg_seq)
436 : 361 : mpi_enreg_seq%my_natom=dtset%natom
437 : :
438 : : !Compute effective number of k-points
439 : : !AZ: two changes for the PAW+collinear ASR bug:
440 : : ! nkpt_me=nkpt_rbz*nsppol instead of nkpt_me=nkpt_rbz,
441 : : ! and ikpt_me is moved outside the LOOP OVER SPINS (see below)
442 : : nkpt_me=nkpt_rbz*nsppol !AZ:multiplied by nsppol
443 : : if(xmpi_paral==1)then
444 : 361 : nkpt_me=0
445 [ + + ]: 757 : do isppol=1,nsppol
446 [ + + ]: 15552 : do ikpt=1,nkpt_rbz
447 : 14795 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
448 [ + + ]: 15191 : if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) nkpt_me=nkpt_me+1
449 : : end do
450 : : end do
451 : : end if
452 : :
453 : : !Sizes for WF at k+q
454 : 361 : mcgq=mpw1*nspinor*mband_mem_rbz*mkqmem*nsppol
455 : 361 : mcprjq=nspinor*mband_mem_rbz*mkqmem*nsppol*usecprj
456 : :
457 [ + + - + : 16239 : ABI_MALLOC(bands_treated_now_ndat, (maxval(nband_rbz),ndat))
+ - - + +
- + - + -
- + - + ]
458 : :
459 : : #ifdef HAVE_OPENMP_OFFLOAD
460 : : !$OMP TARGET ENTER DATA MAP(to:cg,cg1) IF(gpu_option==ABI_GPU_OPENMP)
461 : : #endif
462 : :
463 : : !Check ddk files (needed to compute electric field perturbations)
464 : 361 : ddkfil(:)=0
465 [ + + ]: 1444 : do idir1=1,3
466 : 1083 : ddkcase=idir1+dtset%natom*3
467 : 1083 : call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(idir1))
468 : 1083 : t_exist = file_exists(fiwfddk(idir1))
469 : :
470 [ + - ]: 1083 : if (.not. t_exist) then
471 : : ! Try netcdf file.
472 : 1083 : t_exist = file_exists(nctk_ncify(fiwfddk(idir1)))
473 [ + + ]: 1083 : if (t_exist) then
474 : 402 : fiwfddk(idir1) = nctk_ncify(fiwfddk(idir1))
475 [ + - ]: 402 : write(msg,"(3a)")"- File: ",trim(fiwfddk(idir1))," does not exist but found netcdf file with similar name."
476 : 402 : call wrtout(std_out,msg,'COLL')
477 : : end if
478 : : end if
479 : :
480 : 361 : if (t_exist) then
481 : 402 : ddkfil(idir1)=20+idir1 ! Note the use of unit numbers 21, 22 and 23
482 : : end if
483 : : end do
484 [ + + + + ]: 1042 : has_ddk_file=(any(ddkfil(:)>0))
485 : :
486 : : !Define the set of perturbations (j1)=(ipert1,idir1)
487 : : !The first perturbation must be (j1)=(j2)=(ipert,idir)
488 : : !because we need to compute <g|H^(j2)-Eps.S^(j2)|u0> first.
489 [ + + ]: 361 : if (ipert/=dtset%natom+1) then
490 : 220 : mpert1=0
491 [ - + + - : 660 : ABI_MALLOC(jpert1,(mpert))
+ - + - -
+ - + ]
492 [ + + ]: 2502 : jpert1 = 0
493 [ + + + - ]: 220 : if (ipert/=dtset%natom+2.or.has_ddk_file) then
494 : 220 : mpert1=mpert1+1;jpert1(mpert1)=ipert
495 : : end if
496 [ + + ]: 2502 : do ipert1=1,mpert
497 [ + + ]: 2502 : if (ipert1/=ipert) then
498 [ + - ]: 2062 : if(dtset%usepaw==1) then
499 [ + + + + : 2062 : if((ipert1<=dtset%natom.or.(ipert1==dtset%natom+2.and.has_ddk_file).or.&
+ + + + +
- + - + +
+ + + + -
+ ]
500 : : & ((ipert>dtset%natom.and.ipert/=dtset%natom+5.and.ipert/=dtset%natom+6).and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)).or. &
501 : : & ((ipert1==dtset%natom+2).and.has_ddk_file))) then
502 : 864 : mpert1=mpert1+1;jpert1(mpert1)=ipert1
503 : : end if
504 : : else ! dtset%usepaw==0
505 [ # # # # : 0 : if ((ipert1<=dtset%natom.or.(ipert1==dtset%natom+2.and.has_ddk_file)).or.&
# # # # #
# # # #
# ]
506 : : & ((ipert==dtset%natom+3.or.ipert==dtset%natom+4).and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4))) then
507 : 0 : mpert1=mpert1+1;jpert1(mpert1)=ipert1
508 : : end if
509 : : end if
510 : : end if
511 : : end do
512 : : else
513 : 141 : mpert1=1
514 [ - + ]: 141 : ABI_MALLOC(jpert1,(mpert1))
515 : 141 : jpert1(1)=dtset%natom+1
516 : : end if
517 : 361 : mdir1=3
518 [ - + ]: 361 : ABI_MALLOC(jdir1,(mdir1))
519 [ + + ]: 1444 : jdir1(1:3)= (/ (idir1,idir1=1,3) /)
520 : 361 : jdir1(1)=idir;jdir1(idir)=1
521 : :
522 : : !Index of strain perturbation, if any
523 [ + + ]: 361 : istr=idir;if (ipert==dtset%natom+4) istr=idir+3
524 : :
525 : : !Open ddk WF file(s) in sequential mode
526 [ + + ]: 361 : if (has_ddk_file) then
527 [ + + ]: 536 : do kdir1=1,mdir1
528 : 402 : idir1=jdir1(kdir1)
529 [ + - ]: 536 : if (ddkfil(idir1)/=0) then
530 [ + - ]: 402 : write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(idir1))
531 [ + + ]: 1206 : call wrtout([std_out, ab_out],msg)
532 : 402 : call ddks(idir1)%open_read(fiwfddk(idir1), formeig1, dtset%iomode, ddkfil(idir1), xmpi_comm_self)
533 : : end if
534 : : end do
535 : :
536 [ - + + - : 536 : ABI_MALLOC(cg_ddk,(2,mpw1*nspinor*mband_mem_rbz,3))
+ - + - +
- - + -
+ ]
537 [ + + + + : 4045928 : cg_ddk = zero ! not all may be initialized below if only certain ddk directions are provided
+ + ]
538 : : end if
539 : :
540 : : !Zero only portion of matrix to be computed here
541 [ + + + + : 25251 : d2nl(:,:,1:dtset%natom+4,idir,ipert)=zero
+ + ]
542 [ + - + + : 25251 : if (usepaw==1) d2ovl(:,:,1:dtset%natom+4,idir,ipert)=zero
+ + + + ]
543 : :
544 : : !Update list of computed matrix elements
545 [ + + ]: 1586 : do kpert1=1,mpert1
546 : 1225 : ipert1=jpert1(kpert1)
547 [ + + ]: 5261 : do kdir1=1,mdir1
548 : 3675 : idir1=jdir1(kdir1)
549 : : if ((ipert1<=dtset%natom).or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4.or.&
550 : : & (ipert1==dtset%natom+1.and.((ddkfil(idir1)/=0).or.(dtset%rfdir(idir1)/=0.and.idir1<=idir))).or.&
551 [ + + + + : 3675 : & (ipert1==dtset%natom+2.and.ddkfil(idir1)/=0).or.&
+ + + + +
- + + + +
+ + + - ]
552 : 1225 : & ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.ddkfil(idir1)/=0)) then
553 : 3522 : blkflg(idir1,ipert1,idir,ipert)=1
554 : : end if
555 : : end do
556 : : end do
557 : :
558 : : !Initialize most of the (1st-order) Hamiltonian
559 : : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
560 : : !2) Perform the setup needed for the non-local factors:
561 : : call gs_hamkq%init(psps,pawtab,nspinor,nsppol,nspden,dtset%natom,&
562 : : & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
563 : : & paw_ij=paw_ij,mpi_atmtab=my_atmtab,comm_atom=my_comm_atom,mpi_spintab=mpi_enreg%my_isppoltab,&
564 : 361 : & usecprj=usecprj,nucdipmom=dtset%nucdipmom,gpu_option=gpu_option)
565 : 361 : has_vectornd = (with_vectornd .EQ. 1)
566 [ + + ]: 361 : if(has_vectornd) then
567 : : ! vlocal is needed when vectornd is present, for zora
568 [ - + + - : 342 : ABI_MALLOC(vlocal,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
569 [ - + + - : 228 : ABI_MALLOC(vtrial_,(nfftf,nspden))
- + + - +
- + - - +
- + ]
570 : : ! need a mutable version of vtrial
571 [ - + + + : 16185390 : vtrial_=vtrial
+ + ]
572 [ - + + - : 399 : ABI_MALLOC(vectornd_pac,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc,3))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + - + ]
573 [ - + + - : 285 : ABI_MALLOC(vectornd_pac_idir,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - - + -
+ - + ]
574 : : end if
575 : 361 : with_vxctau = ( usevxctau > 0 )
576 [ + + ]: 361 : if(with_vxctau) then
577 [ - + + - : 42 : ABI_MALLOC(vxctaulocal,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc,4))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + - + ]
578 : : end if
579 : :
580 : : !Variables common to all perturbations
581 [ + - - + : 173419 : arg=maxval(occ_rbz)-minval(occ_rbz)
- - + + +
- - + - -
+ + ]
582 : 361 : qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14)
583 [ + + + + ]: 361 : is_metal=((dtset%occopt>=3.and.dtset%occopt<=8).or.(abs(arg)>tol8))
584 : 361 : is_metal_or_qne0=((is_metal).or.(qne0))
585 [ + - + - ]: 361 : nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
586 [ - + + - : 1805 : ABI_MALLOC(ch1c,(2,mband_mem_rbz,dtset%mband,nkpt_me))
- + + - -
+ + - + -
+ - - + -
+ ]
587 [ - + + - : 1083 : ABI_MALLOC(ch1c_tmp,(2,mband_mem_rbz))
+ - + - -
+ - + ]
588 [ + + + + : 2241312 : ch1c(:,:,:,:)=zero
+ + + + ]
589 : : #ifdef HAVE_OPENMP_OFFLOAD
590 : : !$OMP TARGET ENTER DATA MAP(alloc:ch1c,ch1c_tmp) IF(gpu_option==ABI_GPU_OPENMP)
591 : : !$OMP TARGET UPDATE TO(ch1c) IF(gpu_option==ABI_GPU_OPENMP)
592 : : #endif
593 : 361 : nzlmopt_ipert=0;nzlmopt_ipert1=0
594 [ + - ]: 361 : if (usepaw==1) then
595 [ - + + - : 1805 : ABI_MALLOC(d2ovl_drho,(2,3,mpert,3,mpert))
+ - - + +
- + - + -
- + - + ]
596 [ + + + + : 1102329 : d2ovl_drho=zero
+ + + + +
+ ]
597 [ + - ]: 361 : if (dtset%pawnzlm/=0) then
598 [ + + ]: 361 : nzlmopt_ipert=1;if (dtset%nstep<2) nzlmopt_ipert=-1
599 : 361 : nzlmopt_ipert1=-1
600 : : end if
601 [ + + ]: 361 : if (is_metal_or_qne0) then
602 [ - + + - : 1260 : ABI_MALLOC(cs1c,(2,dtset%mband,mband_mem_rbz,nkpt_me))
- + + - -
+ + - + -
+ - - + -
+ ]
603 [ - + + - : 756 : ABI_MALLOC(cs1c_tmp,(2,dtset%mband))
+ - + - -
+ - + ]
604 [ + + + + : 1961378 : cs1c(:,:,:,:)=zero
+ + + + ]
605 : : end if
606 : : end if
607 : :
608 : : !Force the recomputation of on-site potentials and DijHartree
609 [ + - ]: 361 : if (usepaw==1) then
610 : 361 : call paw_an_reset_flags(paw_an1)
611 : 361 : call paw_ij_reset_flags(paw_ij1,dijhartree=.true.)
612 : : end if
613 : :
614 : : !LOOP OVER PERTURBATION TYPES (j1)
615 [ + + ]: 1586 : do kpert1=1,mpert1
616 : 1225 : ipert1=jpert1(kpert1)
617 : :
618 : : ! Flag for use of DDK file
619 [ + + + - : 1225 : need_ddk_file=(has_ddk_file.and.(ipert1==dtset%natom+1.or.ipert1==dtset%natom+2))
+ + ]
620 : :
621 : : ! Factor to be applied for electric Field (Eff. charges and piezo. tensor are "minus" d2E)
622 : 1225 : elfd_fact=one
623 [ + + + + : 1225 : if ((ipert <=dtset%natom.or.ipert ==dtset%natom+3.or.ipert ==dtset%natom+4).and. &
+ + + + ]
624 : 1225 : (ipert1==dtset%natom+2)) elfd_fact=-one
625 [ + + + + : 1225 : if ((ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and. &
+ + + + ]
626 : 1225 : (ipert ==dtset%natom+2)) elfd_fact=-one
627 : :
628 : : ! We want to compute delta_u^(j1))=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
629 : : ! see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)
630 : 1225 : has_dcwf=.false.;has_dcwf2=.false.;has_drho=.false.
631 [ + - ]: 1225 : if (usepaw==1) then
632 : 1225 : has_dcwf =(ipert1/=dtset%natom+2)
633 : 1225 : has_dcwf2=(ipert /=dtset%natom+2)
634 [ + + + + ]: 1225 : has_drho =(has_dcwf.and.ipert1/=dtset%natom+1)
635 : : end if
636 : :
637 : : ! Select which WF are needed
638 : 1225 : need_wfk=.true.
639 : 1225 : need_wf1=.true.
640 : :
641 : : ! Initialize data for NL 1st-order (j1) hamiltonian
642 : 1225 : call rf_hamkq%init(cplex,gs_hamkq,ipert1,mpi_spintab=[0,0])
643 : :
644 : : ! The following contributions are needed only for non-DDK perturbation:
645 : : ! - Frozen part of 1st-order Dij
646 : : ! - Contribution from local potential to dynamical matrix (due to Vxc^(j1)(tild_nc)+VH^(j1)(tild_nZc))
647 [ + + + - ]: 1225 : if (ipert/=dtset%natom+1.and.ipert1/=dtset%natom+1) then
648 : :
649 : : ! Allocations
650 : 1084 : force_recompute=(usepaw==0) ! This is dangerous...
651 : 1084 : nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
652 [ + + + + : 1084 : nvxc1=0;if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvxc1=nspden
+ + ]
653 [ + + + + ]: 1084 : nvh1=0;if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvh1=1
654 [ - + + - : 4336 : ABI_MALLOC(vxc10,(cplex*nfftf,nvxc1))
+ + + - +
- + + - +
- + ]
655 [ + + + - : 3252 : ABI_MALLOC(vhartr01,(cplex*nfftf*nvh1))
+ - + + -
+ - + ]
656 : 1084 : need_pawij10=(usepaw==1)
657 [ + - ]: 1084 : if (need_pawij10) then
658 [ + + + - : 13413 : ABI_MALLOC(paw_ij10,(my_natom,mdir1))
+ - + + -
+ - + + +
+ + ]
659 [ - + + - : 8672 : ABI_MALLOC(e1kbfr_spin,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,cplex,mdir1,my_nsppol))
- + + - -
+ + - - +
+ - + - -
+ + - + -
+ - - + -
+ - + ]
660 : : else
661 [ # # ]: 0 : ABI_MALLOC(paw_ij10,(0,0))
662 : : end if
663 : :
664 : : ! LOOP OVER PERTURBATION DIRECTIONS
665 [ + + ]: 4336 : do kdir1=1,mdir1
666 : 3252 : idir1=jdir1(kdir1)
667 [ + + ]: 3252 : istr1=idir1;if (ipert1==dtset%natom+4) istr1=idir1+3
668 : :
669 : : ! Get first-order local potential and first-order pseudo core density
670 [ + + + + : 3252 : if (ipert==ipert1.and.idir==idir1.and.(.not.force_recompute)) then
+ - ]
671 : 220 : vpsp1_idir1 => vpsp1
672 : 220 : xccc3d1_idir1 => xccc3d1
673 : : else
674 [ - + + - : 9096 : ABI_MALLOC(vpsp1_idir1,(cplex*nfftf))
+ - + - -
+ - + ]
675 [ + + + - : 9096 : ABI_MALLOC(xccc3d1_idir1,(cplex*n3xccc))
+ - + + -
+ - + ]
676 [ + - ]: 3032 : if (usepaw==1) then
677 : : call dfpt_atm2fft(gs_hamkq%atindx,cplex,gmet,gprimd,gsqcut,istr1,ipert1,&
678 : : & mgfftf,psps%mqgrid_vl,dtset%natom,1,nfftf,ngfftf,dtset%ntypat,&
679 : : & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
680 : : & atmrhor1=xccc3d1_idir1,atmvlocr1=vpsp1_idir1,optv_in=1,optn_in=n3xccc/nfftf,optn2_in=1,&
681 [ - + - + ]: 3032 : & vspl=psps%vlspl)
682 : : else
683 [ # # # # ]: 0 : if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
684 : :
685 : : !To compute Absolute Deformation Potentials together with FxE tensor
686 : : !the reference has to be the same as in the FxE routines
687 [ # # ]: 0 : g0term=0; if (dtset%rfstrs_ref==1) g0term=1
688 : :
689 : : call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
690 : : & nattyp,nfftf,ngfftf,dtset%ntypat,ph1df,psps%qgrid_vl,ucvol,&
691 [ # # ]: 0 : & psps%vlspl,vpsp1_idir1,g0term=g0term)
692 : : else
693 : : call dfpt_vlocal(gs_hamkq%atindx,cplex,gmet,gsqcut,dtset%icutcoul,idir1,ipert1,mpi_enreg,psps%mqgrid_vl,&
694 : : & dtset%natom,nattyp,nfftf,ngfftf,dtset%nkpt,dtset%ntypat,ngfftf(1),ngfftf(2),ngfftf(3),&
695 : : & ph1df,psps%qgrid_vl,dtset%qptn,dtset%rcut,rprimd,ucvol,dtset%vcutgeo,psps%vlspl,vpsp1_idir1,xred,&
696 [ # # ]: 0 : & zion=dtset%ziontypat)
697 : : end if
698 [ # # ]: 0 : if(psps%n1xccc/=0)then
699 : : call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,ngfftf(1),psps%n1xccc,&
700 : : & ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,&
701 [ # # ]: 0 : & psps%xcccrc,psps%xccc1d,xccc3d1_idir1,xred)
702 : : end if
703 : : end if
704 : : end if
705 : :
706 : : ! Compute 1st-order non-local factors (Dij^(j1)_fr)
707 [ + - ]: 3252 : if (need_pawij10) then
708 : 3252 : call paw_ij_nullify(paw_ij10(:,idir1))
709 : : call paw_ij_init(paw_ij10(:,idir1),cplex,nspinor,dtset%nsppol,dtset%nspden,&
710 : : & 0,dtset%natom,dtset%ntypat,dtset%typat,pawtab,has_dijfr=1,&
711 : 3252 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
712 [ + + + + ]: 3252 : if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then
713 : 3032 : option=0
714 : : call pawdijfr(gprimd,idir1,ipert1,my_natom,dtset%natom,nfftf,ngfftf,&
715 : : & nspden,nsppol,dtset%ntypat,option,paw_ij10(:,idir1),pawang,pawfgrtab,pawrad,&
716 : : & pawtab,cplex,dtset%qptn,rprimd,ucvol,vpsp1_idir1,vtrial,vxc,xred,&
717 [ - + - - : 3032 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
- - - + -
- - - - -
- - - - -
- - - - -
- - - + ]
718 : : else
719 [ + + ]: 663 : do iatom=1,my_natom
720 : 443 : paw_ij10(iatom,idir1)%has_dijfr=paw_ij1(iatom)%has_dijfr
721 [ + - + - : 24731 : if (paw_ij1(iatom)%has_dijfr==2) paw_ij10(iatom,idir1)%dijfr=paw_ij1(iatom)%dijfr
+ - - + -
- - - - -
+ + + + ]
722 : : end do
723 : : end if
724 : : end if
725 : :
726 : : ! Get first-order exchange-correlation potential (core-correction contribution only)
727 [ + + ]: 3252 : if (nvxc1>0) then
728 [ + + + + ]: 2850 : if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
729 : 624 : option=0
730 : : call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,&
731 : : & nhat,nhat1,nkxc,nmxc,nspden,n3xccc,option,dtset%qptn,rhor,rhor1,rprimd,&
732 [ - + ]: 624 : & usepaw,usexcnhat,vxc10,xccc3d1_idir1)
733 : : else
734 : : ! Non-collinear magnetism (should the second nkxc be nkxc_cur ?)
735 [ - + ]: 2226 : if (nspden==4) then
736 : 0 : option=0
737 : : call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum1,0,dum2,0,dum3,0,nkxc,&
738 : : & nmxc,nspden,n3xccc,1,option,dtset%qptn,dum1,dum1,rprimd,0,vxc,&
739 [ # # ]: 0 : & vxc10,xccc3d1_idir1)
740 : : else
741 : : call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum2,0,dum3,0,nkxc,&
742 [ - + ]: 2226 : & nmxc,nspden,n3xccc,0,dtset%qptn,dum1,rprimd,0,vxc10,xccc3d1_idir1)
743 : : end if
744 : : end if
745 : : end if
746 : :
747 : : ! Get first-order Hartree potential (metric tensor contribution only)
748 [ + + ]: 3252 : if (nvh1>0) then
749 : : call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,dtset%natom,&
750 [ - + ]: 624 : & nfftf,ngfftf,rhog,rprimd,vhartr01)
751 : : end if
752 : :
753 : : ! Get Hartree + xc + local contributions to dynamical matrix or elastic tensor
754 [ + + + + : 3252 : if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
+ + ]
755 [ - + ]: 2850 : if (usepaw==0) then
756 : : ! vxc1 is integrated with the total 1st-order density (rhor1 )
757 : : ! vpsp1 is integrated with the 1st-order pseudo density (rhor1)
758 : 0 : call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol)
759 [ # # ]: 0 : call dotprod_vn(cplex,rhor1,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol)
760 [ + + ]: 2850 : else if (usexcnhat/=0) then
761 : : ! vxc1 is integrated with the total 1st-order density (rhor1 including nhat1)
762 : : ! vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1)
763 [ - + + - : 6129 : ABI_MALLOC(rhotmp,(cplex*nfftf,1))
+ - + - -
+ - + ]
764 [ + + ]: 25608963 : rhotmp(:,1)=rhor1(:,1)-nhat1(:,1)
765 : 2043 : call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol)
766 [ - + ]: 2043 : call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol)
767 [ - + ]: 2043 : ABI_FREE(rhotmp)
768 : : else
769 : : ! vxc1 is integrated with the 1st-order pseudo density (rhor1 without nhat1)
770 : : ! vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1)
771 [ - + + - : 3228 : ABI_MALLOC(rhotmp,(cplex*nfftf,nspden))
- + + - +
- + - - +
- + ]
772 [ + + + + ]: 38357748 : rhotmp(:,:)=rhor1(:,:)-nhat1(:,:)
773 : 807 : call dotprod_vn(cplex,rhotmp,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol)
774 [ - + ]: 807 : call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol)
775 [ - + ]: 807 : ABI_FREE(rhotmp)
776 : : end if
777 [ + + ]: 2850 : if (nvh1>0) then
778 [ - + ]: 624 : call dotprod_vn(cplex,rhor1,dot3r,dot3i,nfftf,nfftot,1,2,vhartr01,ucvol)
779 : : else
780 : 2226 : dot3r=zero ; dot3i=zero
781 : : end if
782 : : ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2
783 : 2850 : dotr=dot1r+dot2r+dot3r;doti=dot1i+dot2i+dot3i
784 : : ! In case ipert = natom+2, these lines compute the local part
785 : : ! of the Born effective charges from phonon and electric
786 : : ! field type perturbations, see eq. 43 of X. Gonze and C. Lee, PRB 55, 10355 (1997) [[cite:Gonze1997a]]
787 : : ! The minus sign is due to the fact that the effective charges
788 : : ! are minus the second derivatives of the energy
789 : : if (ipert/=dtset%natom+1) then
790 : 2850 : d2lo(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr
791 : 2850 : d2lo(2,idir1,ipert1,idir,ipert)=elfd_fact*doti
792 : : end if
793 : : end if ! ipert1<=natom
794 : :
795 [ + + + + : 4336 : if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then
- + ]
796 [ - + ]: 3032 : ABI_FREE(vpsp1_idir1)
797 [ - + ]: 3032 : ABI_FREE(xccc3d1_idir1)
798 : : end if
799 : :
800 : : ! End loop on directions
801 : : end do
802 : :
803 : : ! Free memory
804 [ - + ]: 1084 : ABI_FREE(vxc10)
805 [ - + ]: 1084 : ABI_FREE(vhartr01)
806 : :
807 : : else ! ddk perturbation
808 [ + + + + ]: 1410 : d2lo(1:2,1:mdir1,ipert1,idir,ipert)=zero
809 : 141 : need_pawij10=.false.
810 [ - + ]: 141 : ABI_MALLOC(paw_ij10,(0,0))
811 : : end if
812 : :
813 : : ! Prepare RF PAW files for reading and writing if mkmem, mkqmem or mk1mem==0
814 : 1225 : iorder_cprj=0
815 : :
816 : : ! Allocate arrays used to accumulate density change due to overlap
817 [ + + ]: 1225 : if (has_drho) then
818 [ - + + - : 5700 : ABI_MALLOC(drhoaug1,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),mdir1))
- + + - -
+ + - + -
+ - + - -
+ - + -
+ ]
819 [ - + + - : 4750 : ABI_MALLOC(drho1wfr,(cplex*dtset%nfft,dtset%nspden,mdir1))
- + + - +
- + - + -
- + - + -
+ ]
820 [ + + + - : 11835 : ABI_MALLOC(pawdrhoij1,(my_natom,mdir1))
+ - + + -
+ - + - +
+ + + + ]
821 [ + + ]: 950 : if (paral_atom) then
822 [ - + + - : 11148 : ABI_MALLOC(pawdrhoij1_unsym,(dtset%natom,mdir1))
+ - + - -
+ - + + +
+ + ]
823 : : else
824 : 402 : pawdrhoij1_unsym => pawdrhoij1
825 : : end if
826 [ + + ]: 3800 : do kdir1=1,mdir1
827 : 2850 : idir1=jdir1(kdir1)
828 [ + + + + ]: 32396880 : drho1wfr(:,:,idir1)=zero
829 : : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
830 : 2850 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
831 : : call pawrhoij_alloc(pawdrhoij1(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,&
832 : : & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,use_rhoijp=1,use_rhoij_=0,&
833 : 2850 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=my_atmtab)
834 [ + + ]: 6650 : if (paral_atom) then
835 : : call pawrhoij_alloc(pawdrhoij1_unsym(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,&
836 : : & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,&
837 : 1644 : & use_rhoijp=0,use_rhoij_=1)
838 : : else
839 : 1206 : call pawrhoij_init_unpacked(pawdrhoij1_unsym(:,idir1))
840 : : end if
841 : : end do
842 : : end if
843 : :
844 : : ! Initialize shifts for global arrays
845 : 1225 : bdtot_index=0
846 : 1225 : bd2tot_index=0
847 : 1225 : ibg=0;icg=0
848 : 1225 : ibg1=0;icg1=0
849 : 1225 : ibgq=0;icgq=0
850 : :
851 : :
852 : : ! Has to get 1st-order non-local factors before the loop over spins
853 : : ! because this needs a communication over comm_atom (=comm_spinkpt)
854 [ + + ]: 1225 : if (need_pawij10) then
855 [ - + ]: 1084 : if (my_nsppol<nsppol) then
856 [ # # # # : 0 : ABI_MALLOC(e1kb_work,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,cplex))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
857 : : end if
858 : 1084 : ii=0
859 [ + + ]: 2234 : do isppol=1,nsppol
860 [ + - ]: 1150 : if (my_spintab(isppol)==1) ii=ii+1
861 [ - + ]: 1150 : if (my_spintab(isppol)/=1) e1kb_ptr => e1kb_work
862 [ + + ]: 5684 : do kdir1=1,mdir1
863 : 3450 : idir1=jdir1(kdir1)
864 [ + - ]: 3450 : if (my_spintab(isppol)==1) e1kb_ptr => e1kbfr_spin(:,:,:,:,idir1,ii)
865 : 4600 : call pawdij2e1kb(paw_ij10(:,idir1),isppol,my_comm_atom,e1kbfr=e1kb_ptr,mpi_atmtab=my_atmtab)
866 : : end do
867 : : end do
868 [ - + ]: 1084 : if (my_nsppol<nsppol) then
869 [ # # ]: 0 : ABI_FREE(e1kb_work)
870 : : end if
871 : : end if
872 : :
873 : : !AZ: do not reset ikpt_me inside the spin loop (PAW+collinear ASR bug)
874 : : ! otherwise the spin down overwrites the spin up component
875 : 1225 : ikpt_me=0
876 : :
877 : : ! LOOP OVER SPINS
878 [ + + ]: 2531 : do isppol=1,nsppol
879 : :
880 : : !ikpt_me=0 !AZ:commented, PAW+collinear ASR bug
881 : :
882 : : ! Rewind (k+G) data if needed
883 : 1306 : ikg=0;ikg1=0
884 : :
885 : : ! Continue to initialize the GS/RF Hamiltonian
886 : 1306 : call gs_hamkq%load_spin(isppol,with_nonlocal=.true.)
887 [ + + ]: 1306 : if (need_pawij10) then
888 : 1150 : ii=min(isppol,size(e1kbfr_spin,6))
889 [ + - ]: 1150 : if (ii>0) e1kbfr => e1kbfr_spin(:,:,:,:,:,ii)
890 : : end if
891 : :
892 : : ! if vectornd is present, set it up for addition to gs_hamkq and rf_hamkq.
893 : : ! Note that it must be done for the three Cartesian directions. Also, the following
894 : : ! code assumes explicitly and implicitly that nvloc = 1. This should eventually be generalized.
895 [ + + ]: 1306 : if(has_vectornd) then
896 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
897 : 60 : & dtset%nspden, gs_hamkq%nvloc, 0, pawfgr, mpi_enreg, vtrial_,vlocal)
898 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
899 : 60 : & dtset%nspden, gs_hamkq%nvloc, 3, pawfgr, mpi_enreg, vectornd,vectornd_pac)
900 : 60 : call gs_hamkq%load_spin(isppol, vlocal=vlocal,vectornd=vectornd_pac)
901 [ + + + + : 6586098 : vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,idir)
+ + + + ]
902 : 60 : call rf_hamkq%load_spin(isppol, vectornd=vectornd_pac_idir)
903 : : end if
904 : : !! add vxctau for mGGA to GS hamiltonian and RF hamiltonian
905 [ + + ]: 1306 : if (with_vxctau) then
906 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
907 : 6 : dtset%nspden, gs_hamkq%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
908 : 6 : call gs_hamkq%load_spin(isppol, vxctaulocal=vxctaulocal)
909 : 6 : call rf_hamkq%load_spin(isppol, vxctaulocal=vxctaulocal)
910 : : end if
911 : :
912 : : ! Initialize accumulation of density
913 [ + + + + : 37004704 : if (has_drho) drhoaug1(:,:,:,:)=zero
+ + + + +
+ ]
914 : :
915 : : ! LOOP OVER K-POINTS
916 [ + + ]: 32947 : do ikpt=1,nkpt_rbz
917 : :
918 : : ! Load dimensions for this k-point
919 : 31641 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
920 : 31641 : istwf_k=istwfk_rbz(ikpt)
921 : 31641 : npw_k=npwarr(ikpt)
922 : 31641 : npw1_k=npwar1(ikpt)
923 : 31641 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
924 : :
925 : : ! Skip loop if this k-point is not to be treated by this proc
926 [ + + ]: 31641 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
927 : 8416 : bdtot_index=bdtot_index+nband_k
928 : 8416 : bd2tot_index=bd2tot_index+2*nband_k**2
929 : 8416 : cycle ! Skip the rest of the k-point loop
930 : : end if
931 : :
932 : : ! Allocate/initialize local arrays and scalars for this k-point
933 [ - + ]: 23225 : ABI_MALLOC(d2nl_k,(2,3))
934 [ - + ]: 23225 : ABI_MALLOC(d2ovl_k,(2,3))
935 [ - + + - : 69675 : ABI_MALLOC(eig_k,(nband_k))
+ - + - -
+ - + ]
936 [ - + + - : 46450 : ABI_MALLOC(eig_kq,(nband_k))
- + - + ]
937 [ + - - + : 69675 : ABI_MALLOC(eig1_k,(2*nband_k**2))
+ - + - +
- - + -
+ ]
938 [ - + + - : 46450 : ABI_MALLOC(occ_k,(nband_k))
- + - + ]
939 [ - + + - : 46450 : ABI_MALLOC(vdotr,(nband_k))
- + - + ]
940 [ - + + - : 46450 : ABI_MALLOC(vdoti,(nband_k))
- + - + ]
941 [ - + + - : 46450 : ABI_MALLOC(vdot1r,(nband_k))
- + - + ]
942 [ - + + - : 46450 : ABI_MALLOC(vdot1i,(nband_k))
- + - + ]
943 [ - + + - : 46450 : ABI_MALLOC(vdot2r,(nband_k))
- + - + ]
944 [ - + + - : 46450 : ABI_MALLOC(vdot2i,(nband_k))
- + - + ]
945 : : #ifdef HAVE_OPENMP_OFFLOAD
946 : : !$OMP TARGET ENTER DATA MAP(alloc:vdotr,vdoti,vdot1r,vdot1i,vdot2r,vdot2i) IF(gpu_option==ABI_GPU_OPENMP)
947 : : #endif
948 [ + + + + ]: 232250 : d2nl_k(:,:)=zero
949 [ + + + + ]: 232250 : d2ovl_k(:,:)=zero
950 [ + + ]: 163505 : eig_k (:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
951 [ + + ]: 163505 : eig_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
952 [ + + ]: 2616305 : eig1_k(:)=eigen1(1+bd2tot_index:2*nband_k**2+bd2tot_index)
953 [ + + ]: 163505 : occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
954 [ + + + + ]: 163505 : nband_kocc=count(abs(occ_k(:))>tol8)
955 [ + + ]: 92900 : kpoint(:)=kpt_rbz(:,ikpt)
956 [ + + + + ]: 87344 : kpq(:)=kpoint(:);if (ipert1<dtset%natom+3) kpq(:)=kpq(:)+dtset%qptn(1:3)
957 : 23225 : wtk_k=wtk_rbz(ikpt)
958 : 23225 : need_ylmgr1=0;dimylmgr1=0
959 : 23225 : nkpg=0;nkpg1=0
960 : 23225 : ikpt_me=ikpt_me+1
961 [ + + ]: 23225 : if (is_metal) then
962 : : ! For each pair of active bands (m,n), generates the ratios
963 : : ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
964 [ - + + - : 11976 : ABI_MALLOC(doccde_k,(nband_k))
- + - + -
+ ]
965 [ - + + - : 11976 : ABI_MALLOC(doccde_kq,(nband_k))
- + - + -
+ ]
966 [ - + + - : 11976 : ABI_MALLOC(occ_kq,(nband_k))
- + - + -
+ ]
967 [ - + + - : 23952 : ABI_MALLOC(rocceig,(nband_k,nband_k))
- + + - +
- + - - +
- + - + ]
968 [ + + ]: 72148 : doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index)
969 [ + + ]: 72148 : doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index)
970 [ + + ]: 72148 : occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index)
971 : 5988 : call occeig(doccde_k,doccde_kq,eig_k,eig_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig)
972 : : end if
973 : :
974 : : ! Take care of the npw and kg records in WF and DDK files
975 [ + + ]: 23225 : if (need_ddk_file) then
976 [ + + ]: 6208 : do kdir1=1,mdir1
977 : 4656 : idir1=jdir1(kdir1)
978 [ + - ]: 6208 : if (ddkfil(idir1)/=0)then
979 : : !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt)
980 : 4656 : ik_ddk = indkpt1(ikpt)
981 : 4656 : npw_ = ddks(idir1)%hdr%npwarr(ik_ddk)
982 [ - + ]: 4656 : if (npw_/=npw_k) then
983 : : write(msg, '(a,i0,a,i0,a,i0,a,a,i0,a,a,i0)')&
984 : 0 : 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,&
985 : 0 : 'the number of plane waves in the ddk file is equal to', npw_,ch10,&
986 : 0 : 'while it should be ',npw_k
987 : 0 : ABI_ERROR(msg)
988 : : end if
989 : :
990 : : ! NB: this will fail if the bands are not contiguous.
991 : 4656 : startband = nband_k
992 : 4656 : endband = 1
993 [ + + ]: 34944 : do iband=1,nband_k
994 [ + - ]: 34944 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then
995 : 30288 : if (iband < startband) startband = iband
996 : 30288 : if (iband > endband) endband = iband
997 : : end if
998 : : end do
999 : : ! NB: eig_k is band distributed in call to read_band_block, though array has full size,
1000 : : ! only certain columns for my iband are filled, then used below
1001 : : call ddks(idir1)%read_band_block((/startband,endband/),ik_ddk,isppol,xmpio_collective, &
1002 [ + + ]: 13968 : & cg_k=cg_ddk(:,:,idir1))
1003 : : end if ! ddk file is already present
1004 : : end do
1005 : : end if
1006 : :
1007 : : ! Allocate arrays used for NL form factors
1008 [ - + + - : 69675 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + -
+ ]
1009 [ - + + - : 69675 : ABI_MALLOC(kg1_k,(3,npw1_k))
+ - + - -
+ - + -
+ ]
1010 [ - + + - : 92900 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
- + + - +
- + - - +
- + - + ]
1011 [ - + + - : 92900 : ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
- + + - +
- + - - +
- + - + ]
1012 [ - + + + : 23225 : if (psps%useylm==1.and.(need_ddk_file.or.ipert1==dtset%natom+1)) need_ylmgr1=1
+ + ]
1013 [ - + + + : 23225 : if (psps%useylm==1.and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)) need_ylmgr1=1
+ + ]
1014 : 23225 : dimylmgr1=max(useylmgr1,need_ylmgr1)
1015 [ - + + - : 92900 : ABI_MALLOC(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*dimylmgr1))
+ + + - +
- + + - +
- + - + ]
1016 : :
1017 : : ! Get plane-wave vectors and related data at k
1018 [ + + + + ]: 13883093 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
1019 [ + - ]: 23225 : if (psps%useylm==1) then
1020 [ + + ]: 122115 : do jj=1,psps%mpsang*psps%mpsang
1021 [ + + ]: 14970833 : ylm_k(1:npw_k,jj)=ylm(1+ikg:npw_k+ikg,jj)
1022 : : end do
1023 : : end if
1024 : :
1025 : : ! Get plane-wave vectors and related data at k+q
1026 [ + + + + ]: 13883093 : kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
1027 [ + - ]: 23225 : if (psps%useylm==1) then
1028 [ + + ]: 122115 : do jj=1,psps%mpsang*psps%mpsang
1029 [ + + ]: 14970833 : ylm1_k(1:npw1_k,jj)=ylm1(1+ikg1:npw1_k+ikg1,jj)
1030 : : end do
1031 [ + + + + ]: 23225 : if (need_ylmgr1==1.and.useylmgr1/=0) then
1032 [ + + ]: 35985 : do jj=1,psps%mpsang*psps%mpsang
1033 [ + + ]: 124023 : do ia=1,3
1034 [ + + ]: 17815470 : ylmgr1_k(1:npw1_k,ia,jj)=ylmgr1(1+ikg1:npw1_k+ikg1,ia,jj)
1035 : : end do
1036 : : end do
1037 : : end if
1038 : : end if
1039 : :
1040 : : ! If Ylm gradients at k+q are needed and not in memory, compute them
1041 [ + - + + ]: 7391 : if (need_ylmgr1==1.and.useylmgr1==0) then
1042 : 752 : option=-1;npwar1_tmp(1)=npw1_k;nband_tmp(1)=nband_k
1043 : : !Subtlety: initylmg is called in sequential mode
1044 : : call initylmg(gprimd,kg1_k,kpq,1,mpi_enreg_seq,psps%mpsang,&
1045 : : & npw1_k,nband_tmp,1,npwar1_tmp,nsppol,option,rprimd,&
1046 : 752 : & ylm1_k,ylmgr1_k)
1047 : : end if
1048 : :
1049 : : ! Compute (k+G) vectors
1050 [ + + ]: 23225 : nkpg=0;if(ipert1<=dtset%natom) nkpg=3*dtset%nloalg(3)
1051 [ - + + - : 92900 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ + + - +
- + + - +
- + - + ]
1052 [ + + ]: 23225 : if (nkpg>0) then
1053 : 15834 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1054 : : end if
1055 : :
1056 : : ! Compute (k+q+G) vectors
1057 [ + + + - ]: 23225 : nkpg1=0;if(ipert1<=dtset%natom.or.need_ylmgr1==1) nkpg1=3*dtset%nloalg(3)
1058 [ - + + - : 92900 : ABI_MALLOC(kpg1_k,(npw1_k,nkpg1))
- + + - +
- + - - +
- + - + ]
1059 [ + - ]: 23225 : if (nkpg1>0) then
1060 : 23225 : call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k)
1061 : : end if
1062 : :
1063 : : ! Allocate kinetic contributions
1064 [ - + + - : 69675 : ABI_MALLOC(dkinpw,(npw_k))
+ - + - -
+ - + -
+ ]
1065 [ - + + - : 69675 : ABI_MALLOC(kinpw1,(npw1_k))
+ - + - -
+ - + -
+ ]
1066 [ + + ]: 3488192 : dkinpw=zero
1067 : : ! Compute (1/2) (2 Pi)**2 (k+q+G)**2:
1068 : : ! call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k)
1069 : 23225 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
1070 : :
1071 : : ! Compute nonlocal form factors ffnl at (k+G), for all atoms
1072 : 23225 : ider=0;idir0=0
1073 [ + + ]: 23225 : dimffnl=0;if (ipert1<=dtset%natom) dimffnl=1
1074 [ - + + - : 139350 : ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
+ + + - -
+ + - - +
+ - + - +
+ - + - +
- + ]
1075 [ + + ]: 23225 : if (ipert1<=dtset%natom) then
1076 : : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,&
1077 : : & ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,&
1078 : : & psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,&
1079 : 15834 : & psps%useylm,ylm_k,ylmgr_dum)
1080 : : #ifdef HAVE_OPENMP_OFFLOAD
1081 : : !$OMP TARGET ENTER DATA MAP(to:ffnlk) IF(gpu_option==ABI_GPU_OPENMP)
1082 : : #endif
1083 : : end if
1084 : :
1085 : : ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
1086 [ + + ]: 23225 : ider=1;if (ipert1<=dtset%natom) ider=0
1087 [ + + + + ]: 23225 : if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)then
1088 [ + - ]: 1852 : dimffnl1=1;if (ider>=1) dimffnl1=2+5*psps%useylm
1089 [ + - + - ]: 1852 : idir0=0;if (ider>0.and.psps%useylm==1) idir0=-7
1090 : : else
1091 [ + + ]: 21373 : dimffnl1=1;if (ider>=1) dimffnl1=2+2*psps%useylm
1092 [ + + + - ]: 21373 : idir0=0;if (ider>0.and.psps%useylm==1) idir0=4
1093 : : end if
1094 [ - + + - : 139350 : ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1095 : : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,&
1096 : : & ider,idir0,psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,&
1097 : : & psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,&
1098 : 23225 : & psps%useylm,ylm1_k,ylmgr1_k)
1099 : : #ifdef HAVE_OPENMP_OFFLOAD
1100 : : !$OMP TARGET ENTER DATA MAP(to:ffnl1) IF(gpu_option==ABI_GPU_OPENMP)
1101 : : #endif
1102 : :
1103 : : ! Extract non-local form factors for H^(j1)
1104 [ + + ]: 23225 : if (ipert1<=dtset%natom) then
1105 : 15834 : ffnl1_idir1 => ffnl1(:,:,:,:)
1106 : 15834 : dimffnl1_idir1=dimffnl1
1107 : : else
1108 : 7391 : dimffnl1_idir1=1+ider
1109 [ - + + - : 44346 : ABI_MALLOC(ffnl1_idir1,(npw1_k,dimffnl1_idir1,psps%lmnmax,psps%ntypat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1110 [ - + ]: 7391 : ii=1;if (psps%useylm==0) ii=1+ider
1111 [ + + ]: 20842 : do itypat=1,psps%ntypat
1112 [ + + ]: 137436 : do ilmn=1,psps%lmnmax
1113 [ + - + + : 42771229 : ffnl1_idir1(1:npw1_k,1:ii,ilmn,itypat)=ffnl1(1:npw1_k,1:ii,ilmn,itypat)
+ + + + +
+ ]
1114 : : end do
1115 : : end do
1116 : : #ifdef HAVE_OPENMP_OFFLOAD
1117 : : !$OMP TARGET ENTER DATA MAP(to:ffnl1_idir1) IF(gpu_option==ABI_GPU_OPENMP)
1118 : : #endif
1119 : : end if
1120 : :
1121 : : ! Load k-dependent part in the Hamiltonian datastructure
1122 [ - + + - : 92900 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk))
- + + - +
- + - - +
- + - + ]
1123 : : call gs_hamkq%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,&
1124 : 23225 : & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
1125 : : #ifdef HAVE_OPENMP_OFFLOAD
1126 : : !$OMP TARGET ENTER DATA MAP(to:ph3d) IF(gpu_option==ABI_GPU_OPENMP)
1127 : : #endif
1128 [ + + + + ]: 116125 : if (size(ffnlk)>0) then
1129 : 15834 : call gs_hamkq%load_k(ffnl_k=ffnlk)
1130 : : else
1131 : 7391 : call gs_hamkq%load_k(ffnl_k=ffnl1)
1132 : : end if
1133 : :
1134 : : ! Load k+q-dependent part in the Hamiltonian datastructure
1135 : : ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
1136 : : call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
1137 : : & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1_idir1,&
1138 : 23225 : & compute_gbound=.true.)
1139 [ + + ]: 23225 : if (qne0) then
1140 [ - + + - : 45440 : ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk))
- + + - +
- + - - +
- + - + ]
1141 : 11360 : call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.)
1142 : : end if
1143 : :
1144 : : ! Load k-dependent part in the 1st-order Hamiltonian datastructure
1145 : 23225 : call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw)
1146 : :
1147 : : ! Allocate memory space for one band
1148 [ - + + - : 69675 : ABI_MALLOC(gh1,(2,npw1_k*nspinor*ndat))
+ - + - -
+ - + -
+ ]
1149 : : #ifdef HAVE_OPENMP_OFFLOAD
1150 : : !$OMP TARGET ENTER DATA MAP(alloc:gh1) IF(gpu_option==ABI_GPU_OPENMP)
1151 : : #endif
1152 : 23225 : nullify(cwaveprj0_idir1)
1153 [ + - ]: 23225 : if (usecprj==1) then
1154 [ - + + - : 168462 : ABI_MALLOC(cwaveprj0,(dtset%natom,nspinor*ndat))
- + + - +
- + - - +
- + - + +
+ + + ]
1155 : 23225 : call pawcprj_alloc(cwaveprj0,ncpgr,gs_hamkq%dimcprj)
1156 : : end if
1157 [ + + ]: 23225 : if (has_dcwf) then
1158 [ - + + - : 43346 : ABI_MALLOC(gs1,(2,npw1_k*nspinor*ndat))
- + - + ]
1159 : : #ifdef HAVE_OPENMP_OFFLOAD
1160 : : !$OMP TARGET ENTER DATA MAP(alloc:gs1) IF(gpu_option==ABI_GPU_OPENMP)
1161 : : #endif
1162 : : else
1163 [ - + ]: 1552 : ABI_MALLOC(gs1,(1,1))
1164 : : end if
1165 : :
1166 [ - + + - : 69675 : ABI_MALLOC(band_procs, (nband_k))
+ - + - -
+ - + ]
1167 : : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,nband_k,&
1168 : 23225 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
1169 : :
1170 [ - + + - : 46450 : ABI_MALLOC(gvnlx1, (2,npw1_k*nspinor*ndat))
- + - + -
+ ]
1171 [ - + + - : 69675 : ABI_MALLOC(gvnlx1_tmp,(2,npw1_k*nspinor))
+ - + - -
+ - + -
+ ]
1172 : : #ifdef HAVE_OPENMP_OFFLOAD
1173 : : !$OMP TARGET ENTER DATA MAP(alloc:gvnlx1(1:2,1:npw1_k*nspinor*ndat)) IF(gpu_option==ABI_GPU_OPENMP)
1174 : : !$OMP TARGET ENTER DATA MAP(alloc:gvnlx1_tmp(1:2,1:npw1_k*nspinor)) IF(gpu_option==ABI_GPU_OPENMP)
1175 : : #endif
1176 [ + + ]: 23225 : if (has_dcwf.and.is_metal_or_qne0) then
1177 [ - + + - : 48912 : ABI_MALLOC(gvnlx2, (2,npw1_k*nspinor*nband_k))
+ - + - -
+ - + -
+ ]
1178 : : #ifdef HAVE_OPENMP_OFFLOAD
1179 : : !$OMP TARGET ENTER DATA MAP(alloc:gvnlx2(1:2,1:npw1_k*nspinor*nband_k)) IF(gpu_option==ABI_GPU_OPENMP)
1180 : : #endif
1181 : : end if
1182 : :
1183 : : ! LOOP OVER BANDS
1184 : 23225 : iband_me = 0
1185 [ - + - - : 163505 : do iband=1,nband_k,ndat
- + + + ]
1186 : :
1187 : : ! Skip band if not to be treated by this proc
1188 : : if (xmpi_paral==1) then
1189 [ + + ]: 140280 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle
1190 : : end if
1191 : 138300 : iband_me = iband_me + ndat
1192 : :
1193 [ + - ]: 138300 : if(gpu_option==ABI_GPU_DISABLED) then
1194 [ + + + + ]: 3852420 : ch1c_tmp(:,:) = zero
1195 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1196 : 0 : call gpu_set_to_zero(ch1c_tmp,int(2,c_size_t)*mband_mem_rbz)
1197 : : end if
1198 : :
1199 : : ! Extract GS wavefunctions
1200 : : if (need_wfk) then
1201 : : cwave0(1:2,1:npw_k*nspinor*ndat) => &
1202 : 138300 : & cg(:,1+(iband_me-ndat)*npw_k*nspinor+icg:iband_me*npw_k*nspinor+icg)
1203 [ + - ]: 138300 : if (usecprj==1) then
1204 [ + + ]: 276600 : do idat=1,ndat
1205 : : call pawcprj_get(gs_hamkq%atindx1,cwaveprj0(:,1+(idat-1)*nspinor:idat*nspinor),cprj,dtset%natom,iband_me-(ndat-idat),ibg,ikpt,iorder_cprj,&
1206 : 276600 : & isppol,mband_mem_rbz,mkmem,dtset%natom,1,nband_me,nspinor,nsppol,dtfil%unpaw)
1207 : : end do
1208 : : ! in distributed cprj memory, no need for these? cg and cprj have same distribution
1209 : : !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1210 : : end if
1211 : : end if
1212 : :
1213 : : ! Extract 1st-order wavefunctions
1214 : : if (need_wf1) then
1215 : : cwavef(1:2,1:npw1_k*nspinor*ndat) => &
1216 : 138300 : & cg1(:,1+(iband_me-ndat)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)
1217 : : end if
1218 : :
1219 : : ! LOOP OVER PERTURBATION DIRECTIONS
1220 [ + + ]: 576425 : do kdir1=1,mdir1
1221 : 414900 : idir1=jdir1(kdir1)
1222 [ + + ]: 414900 : istr1=idir1;if(ipert1==dtset%natom+4) istr1=idir1+3
1223 : :
1224 : : ! Not able to compute if ipert1=(Elect. field) and no ddk WF file
1225 [ + + + - ]: 414900 : if (ipert1==dtset%natom+2.and.ddkfil(idir1)==0) cycle
1226 : :
1227 [ + - ]: 414900 : if(gpu_option==ABI_GPU_DISABLED) then
1228 [ + + + + ]: 228217878 : gvnlx1(:,:)=zero
1229 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1230 : 0 : call gpu_set_to_zero(gvnlx1,int(2,c_size_t)*npw1_k*nspinor*ndat)
1231 : : end if
1232 : :
1233 : : ! Extract 1st-order NL form factors derivatives for this idir1
1234 [ + + + - : 414900 : if (dimffnl1_idir1>=2.and.psps%useylm==1.and.ipert1>dtset%natom) then
+ - ]
1235 [ + + ]: 447684 : do itypat=1,psps%ntypat
1236 [ + + ]: 3306000 : do ilmn=1,psps%lmnmax
1237 [ + - + + : 1100353956 : ffnl1_idir1(1:npw1_k,2,ilmn,itypat)=ffnl1(1:npw1_k,1+istr1,ilmn,itypat)
+ + ]
1238 : : end do
1239 : : end do
1240 : : #ifdef HAVE_OPENMP_OFFLOAD
1241 : : !$OMP TARGET UPDATE TO(ffnl1_idir1) IF(gpu_option==ABI_GPU_OPENMP)
1242 : : #endif
1243 : : end if
1244 : :
1245 : : ! Setup gemm_nonlop
1246 [ - + ]: 414900 : if (gemm_nonlop_use_gemm) then
1247 : : call set_gemm_nonlop_ikpt(ikpt,gs_hamkq%npw_fft_k,gs_hamkq%istwf_k,gs_hamkq%indlmn,&
1248 : 0 : & gs_hamkq%ntypat,gs_hamkq%nattyp,gs_hamkq%gpu_option)
1249 : : end if ! gemm_nonlop_use_gemm
1250 : :
1251 : : ! Extract ground state projected WF and derivatives in idir1 direction
1252 [ + - ]: 414900 : if (usecprj==1) then
1253 : 414900 : cpopt=1
1254 : :
1255 : : ! === Atomic displ. perturbation
1256 [ + + ]: 414900 : if (ipert<=dtset%natom) then
1257 [ - + + - : 1536624 : ABI_MALLOC(cwaveprj0_idir1,(dtset%natom,nspinor*ndat))
- + + - +
- + - - +
- + + + +
+ ]
1258 : 211536 : call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj)
1259 [ + + ]: 211536 : if (ipert1<=dtset%natom) then
1260 : 197328 : call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1)
1261 : : else
1262 [ + - ]: 14208 : if (ipert1==dtset%natom+2) then
1263 : 14208 : idir_cprj=idir1;choice=5
1264 : : end if
1265 [ + - - + ]: 14208 : if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1266 : 0 : idir_cprj=istr1;choice=3
1267 : : end if
1268 : 14208 : call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1)
1269 : : call getcprj(choice,cpopt,cwave0,&
1270 : : & cwaveprj0_idir1,&
1271 : : & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,&
1272 : : & gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,&
1273 : : & gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,&
1274 : : & gs_hamkq%nloalg,gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,&
1275 [ + - - - : 14208 : & gs_hamkq%ph1d,gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm,is_kprime=.true.,gpu_option=gpu_option)
- - - - +
- - - - -
- - - + -
+ - + ]
1276 : : end if
1277 : :
1278 : : ! === Wave-vector perturbation
1279 [ + + ]: 203364 : else if (ipert==dtset%natom+1) then
1280 : 94608 : cwaveprj0_idir1 => cwaveprj0
1281 : :
1282 : : ! == Electric field perturbation
1283 [ + + ]: 108756 : else if (ipert==dtset%natom+2) then
1284 [ - + + - : 465360 : ABI_MALLOC(cwaveprj0_idir1,(dtset%natom,nspinor*ndat))
- + + - +
- + - - +
- + + + +
+ ]
1285 : 62640 : call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj)
1286 [ + + ]: 62640 : if (ipert1==dtset%natom+2) then
1287 : 11760 : call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1)
1288 : : else
1289 [ + + ]: 50880 : if (ipert1<=dtset%natom) then
1290 : 27360 : idir_cprj=idir1;choice=2
1291 : : end if
1292 [ + + + + ]: 50880 : if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1293 : 23520 : idir_cprj=istr1;choice=3
1294 : : end if
1295 : 50880 : call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1)
1296 : : call getcprj(choice,cpopt,cwave0,&
1297 : : & cwaveprj0_idir1,&
1298 : : & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,&
1299 : : & gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,&
1300 : : & gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,&
1301 : : & gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,&
1302 [ + - - - : 50880 : & gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm,is_kprime=.true.,gpu_option=gpu_option)
- - - - +
- - - - -
- - - + -
+ - + ]
1303 : : end if
1304 : :
1305 : : ! === Strain perturbation
1306 [ + + + - ]: 46116 : else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
1307 [ + + + + : 46116 : if ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.(istr==istr1)) then
+ + ]
1308 : 2484 : cwaveprj0_idir1 => cwaveprj0
1309 : : else
1310 [ - + + - : 384912 : ABI_MALLOC(cwaveprj0_idir1,(dtset%natom,nspinor*ndat))
- + + - +
- + - - +
- + + + +
+ ]
1311 : 43632 : call pawcprj_alloc(cwaveprj0_idir1,ncpgr,gs_hamkq%dimcprj)
1312 [ + + ]: 43632 : if (ipert1<=dtset%natom) then
1313 : 26892 : idir_cprj=idir1;choice=2
1314 : : end if
1315 [ + + ]: 43632 : if (ipert1==dtset%natom+2) then
1316 : 4320 : idir_cprj=idir1;choice=5
1317 : : end if
1318 [ + + + + ]: 43632 : if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1319 : 12420 : idir_cprj=istr1;choice=3
1320 : : end if
1321 : 43632 : call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1)
1322 : : call getcprj(choice,cpopt,cwave0,&
1323 : : & cwaveprj0_idir1,&
1324 : : & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,gs_hamkq%kg_kp,&
1325 : : & gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,ndat,&
1326 : : & gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,gs_hamkq%npw_kp,&
1327 : : & gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,gs_hamkq%ph3d_kp,&
1328 [ + - - - : 43632 : & gs_hamkq%ucvol,gs_hamkq%useylm,is_kprime=.true.,gpu_option=gpu_option)
- - - - +
- - - - -
- - - + -
+ - + ]
1329 : : end if
1330 : : end if ! ipert
1331 : :
1332 : : else ! usecprj=0: cwaveprj0_idir1 is not used
1333 : 0 : cwaveprj0_idir1 => cwaveprj0
1334 : : end if
1335 : :
1336 : : ! Eventually compute 1st-order kinetic operator
1337 [ + + ]: 414900 : if (ipert1==dtset%natom+1) then
1338 : 94608 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0)
1339 [ + + + + ]: 320292 : else if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1340 : 38424 : call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr1,kg_k,kpoint,npw_k)
1341 : : end if
1342 : :
1343 : : ! Finalize initialization of 1st-order NL hamiltonian
1344 [ + + ]: 414900 : if (need_pawij10) rf_hamkq%e1kbfr => e1kbfr(:,:,:,:,idir1)
1345 : :
1346 : : ! Read DDK wave function (if ipert1=electric field)
1347 [ + + ]: 414900 : if (ipert1==dtset%natom+2) then
1348 : 30288 : usevnl=1
1349 [ + - ]: 30288 : if (need_ddk_file) then
1350 [ + - ]: 30288 : if (ddkfil(idir1)/=0) then
1351 : : !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt)
1352 : : !ik_ddk = indkpt1(ikpt)
1353 : : !call ddks(idir1)%read_bks(iband, ik_ddk, isppol, xmpio_single, cg_bks=gvnlx1)
1354 [ + - + - : 17863404 : gvnlx1 = cg_ddk(:,1+(iband_me-ndat)*npw_k*nspinor:iband_me*npw_k*nspinor,idir1)
- + - - +
+ + + ]
1355 : : else
1356 [ # # # # ]: 0 : gvnlx1=zero
1357 : : end if
1358 [ + - ]: 30288 : if (ipert1==dtset%natom+2) then
1359 [ + + ]: 5964564 : do ii=1,npw1_k*nspinor*ndat ! Multiply ddk by +i (to be consistent with getgh1c)
1360 : 5934276 : arg=gvnlx1(1,ii)
1361 : 5934276 : gvnlx1(1,ii)=-gvnlx1(2,ii)
1362 : 5964564 : gvnlx1(2,ii)=arg
1363 : : end do
1364 : : end if
1365 : : else
1366 [ # # # # ]: 0 : gvnlx1=zero
1367 : : end if
1368 : : #ifdef HAVE_OPENMP_OFFLOAD
1369 : : !$OMP TARGET UPDATE TO(gvnlx1) IF(gpu_option==ABI_GPU_OPENMP)
1370 : : #endif
1371 : : else
1372 : 384612 : usevnl=0
1373 : : end if
1374 : :
1375 : : ! Get |H^(j2)-Eps_k_i.S^(j2)|u0_k_i> (VHxc-dependent part not taken into account) and S^(j2)|u0>
1376 : 414900 : lambda=eig_k(iband);berryopt=0;optlocal=0
1377 [ + + + + ]: 414900 : optnl=0;if (ipert1/=dtset%natom+1.or.idir==idir1) optnl=1
1378 [ + + ]: 414900 : opt_gvnlx1=0;if (ipert1==dtset%natom+2) opt_gvnlx1=2
1379 [ + + ]: 414900 : sij_opt=-1;if (has_dcwf) sij_opt=1
1380 [ - + ]: 414900 : if (usepaw==0) sij_opt=0
1381 : : call getgh1c(berryopt,cwave0,&
1382 : : & cwaveprj0_idir1,&
1383 : : & gh1,dum1,&
1384 : : & gs1,gs_hamkq,&
1385 : : & gvnlx1,&
1386 : : & idir1,ipert1,&
1387 : 414900 : & eig_k(iband:iband+ndat-1),mpi_enreg,ndat,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
1388 [ + + + + ]: 414900 : if (sij_opt==1.and.optnl==1) then
1389 [ + - ]: 321540 : if(gpu_option==ABI_GPU_DISABLED) then
1390 [ + + ]: 643080 : do idat=1,ndat
1391 : : gh1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor) = &
1392 [ + + + + ]: 177529530 : & gh1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor)-eig_k(iband+idat-1) * gs1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor)
1393 : : end do
1394 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1395 : : #ifdef HAVE_OPENMP_OFFLOAD
1396 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) PRIVATE(ipw,idat) MAP(to:gs1,gh1,eig_k)
1397 : : #endif
1398 [ # # ]: 0 : do idat=1,ndat
1399 [ # # ]: 0 : do ispinor=1,nspinor
1400 [ # # ]: 0 : do ipw=1,npw1_k
1401 : : gh1(1,(idat-1)*nspinor*npw1_k+(ispinor-1)*npw1_k+ipw) = &
1402 : : & gh1(1,(idat-1)*nspinor*npw1_k+(ispinor-1)*npw1_k+ipw) &
1403 : 0 : & - eig_k(iband+idat-1) * gs1(1,(idat-1)*nspinor*npw1_k+(ispinor-1)*npw1_k+ipw)
1404 : : gh1(2,(idat-1)*nspinor*npw1_k+(ispinor-1)*npw1_k+ipw) = &
1405 : : & gh1(2,(idat-1)*nspinor*npw1_k+(ispinor-1)*npw1_k+ipw) &
1406 : 0 : & - eig_k(iband+idat-1) * gs1(2,(idat-1)*nspinor*npw1_k+(ispinor-1)*npw1_k+ipw)
1407 : : end do
1408 : : end do
1409 : : end do
1410 : : end if
1411 : : end if
1412 : :
1413 : : ! If needed, compute here <delta_u^(j1)_k_i|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>
1414 : : ! with delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
1415 : : ! (see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42))
1416 : : ! This can be rewritten as:
1417 : : ! -1/2.<u0_k_i|S^(j1)| Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j>
1418 : : ! The sum over j can be computed with a single call to projbd routine
1419 : : ! At first call (when j1=j2), ch1c=<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i> is stored
1420 : : ! For the next calls, it is reused.
1421 [ + + + + : 414900 : if (has_dcwf.or.(ipert==ipert1.and.idir==idir1.and.usepaw==1)) then
+ + + - ]
1422 : : ! note: gvnlx1 used as temporary space
1423 [ + + + + ]: 388532 : if (ipert==ipert1.and.idir==idir1) then
1424 : 70972 : option=0;
1425 [ + - ]: 70972 : if(gpu_option==ABI_GPU_DISABLED) then
1426 [ + - + - : 31896758 : gvnlx1=gh1
- + - - -
- - - + +
+ + ]
1427 : : else if(gpu_option==ABI_GPU_OPENMP) then
1428 : : #ifdef HAVE_OPENMP_OFFLOAD
1429 : : call gpu_copy(gvnlx1, gh1, int(2,c_size_t)*npw1_k*nspinor*ndat)
1430 : : #endif
1431 : : end if
1432 : : else
1433 : 317560 : option=1;
1434 [ + - ]: 317560 : if(gpu_option==ABI_GPU_DISABLED) then
1435 [ + + + + ]: 180376624 : gvnlx1=zero
1436 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1437 : 0 : call gpu_set_to_zero(gvnlx1,int(2,c_size_t)*npw1_k*nspinor*ndat)
1438 : : end if
1439 : : end if
1440 : :
1441 : :
1442 : :
1443 [ + + ]: 777064 : do idat=1,ndat
1444 : :
1445 : 388532 : bands_treated_now => bands_treated_now_ndat(:,idat)
1446 [ + + ]: 3923992 : bands_treated_now = 0
1447 : 388532 : bands_treated_now(iband+idat-1) = 1
1448 : 388532 : call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr)
1449 [ + + ]: 4312524 : do iband_=1, nband_k
1450 [ + + ]: 3535460 : if (bands_treated_now(iband_) == 0) cycle
1451 : :
1452 : : ! distribute gvnlx1 to my subcomm
1453 [ + - ]: 394472 : if(gpu_option==ABI_GPU_DISABLED) then
1454 [ + + + + ]: 245566418 : gvnlx1_tmp = zero
1455 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1456 : 0 : call gpu_set_to_zero(gvnlx1_tmp,int(2,c_size_t)*npw1_k*nspinor)
1457 : : end if
1458 [ + + ]: 394472 : if (iband_ == iband+idat-1) then
1459 [ + - ]: 388532 : if(gpu_option==ABI_GPU_DISABLED) then
1460 [ + - + - : 212590942 : gvnlx1_tmp = gvnlx1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor)
- + - - +
+ + + ]
1461 : : else if(gpu_option==ABI_GPU_OPENMP) then
1462 : : #ifdef HAVE_OPENMP_OFFLOAD
1463 : : call gpu_copy(gvnlx1_tmp,&
1464 : : & gvnlx1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor),&
1465 : : & int(2,c_size_t)*npw1_k*nspinor)
1466 : : #endif
1467 : : end if
1468 : : end if
1469 : :
1470 : : ! TODO CHECK IF IT IS BAND_PROCS(IBAND_)
1471 : : !call xmpi_bcast(gvnlx1_tmp, band_procs(iband_), mpi_enreg%comm_band, ierr)
1472 [ + + ]: 394472 : if(xmpi_comm_size(mpi_enreg%comm_band)>1) then
1473 : 7920 : call xmpi_sum(gvnlx1_tmp,mpi_enreg%comm_band,ierr,use_omp_map=use_ompgpu)
1474 : : end if
1475 : :
1476 [ + + ]: 394472 : if (option == 1) then
1477 : : ! in case I need to reuse the ch1c (option 1) then load them here
1478 [ + - ]: 323248 : if(gpu_option==ABI_GPU_DISABLED) then
1479 [ + + + + ]: 8771560 : ch1c_tmp(:,1:nband_me) = ch1c(:,1:nband_me,iband_,ikpt_me)
1480 : : else
1481 : : #ifdef HAVE_OPENMP_OFFLOAD
1482 : : call gpu_copy(ch1c_tmp, ch1c(:,1:nband_me,iband_,ikpt_me),&
1483 : : & int(2,c_size_t)*nband_me)
1484 : : #endif
1485 : : end if
1486 : : end if
1487 : :
1488 : :
1489 : : ! Compute -Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j>
1490 : : call projbd(cgq,gvnlx1_tmp,-1,icgq,0,istwf_k,mcgq,0,nband_me,npw1_k,nspinor,&
1491 : 394472 : & dum1,ch1c_tmp,option,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft,gpu_option=gpu_option)
1492 : :
1493 : : !sum over all jband by combining the projbd
1494 [ + + ]: 394472 : if(xmpi_comm_size(mpi_enreg%comm_band)>1) then
1495 : 7920 : call xmpi_sum(gvnlx1_tmp,mpi_enreg%comm_band,ierr,use_omp_map=use_ompgpu)
1496 : : end if
1497 : :
1498 : : ! keep my own gvnlx
1499 [ + + ]: 394472 : if (iband_ == iband+idat-1) then
1500 : : ! if bands are parallelized, I have only projected against bands on my cpu
1501 : : ! Pc|work> = |work> - Sum_l <psi_{k+q, l}|work> |psi_{k+q, l}>
1502 : : ! = Sum_nproc_band (|work> - Sum_{my l} <psi_{k+q, l}|work> |psi_{k+q, l}>) - (nproc_band-1) |work>
1503 : : !TODO: make this a blas call? zaxpy
1504 [ + - ]: 388532 : if(gpu_option==ABI_GPU_DISABLED) then
1505 : : gvnlx1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor) = &
1506 [ + + + + ]: 212202410 : gvnlx1_tmp - (my_nproc_band-1)*gvnlx1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor)
1507 : : else if(gpu_option==ABI_GPU_OPENMP) then
1508 : : #ifdef HAVE_OPENMP_OFFLOAD
1509 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(ipw) MAP(to:gvnlx1,gvnlx1_tmp)
1510 : : do ipw=1,npw1_k*nspinor
1511 : : gvnlx1(1,(idat-1)*npw1_k*nspinor + ipw) = &
1512 : : gvnlx1_tmp(1,ipw) - (my_nproc_band-1)*gvnlx1(1,(idat-1)*npw1_k*nspinor + ipw)
1513 : : gvnlx1(2,(idat-1)*npw1_k*nspinor + ipw) = &
1514 : : gvnlx1_tmp(2,ipw) - (my_nproc_band-1)*gvnlx1(2,(idat-1)*npw1_k*nspinor + ipw)
1515 : : end do
1516 : : #endif
1517 : : end if
1518 : : end if
1519 : :
1520 [ + + ]: 783004 : if (option == 0) then
1521 : : ! save ch1c for all of the iband_ on each proc, for later use. First band index only for my nband_me which matches cgq
1522 [ + - ]: 71224 : if(gpu_option==ABI_GPU_DISABLED) then
1523 [ + + + + ]: 2229292 : ch1c(:,1:nband_me,iband_,ikpt_me) = ch1c_tmp(:,1:nband_me)
1524 : : else if(gpu_option==ABI_GPU_OPENMP) then
1525 : : #ifdef HAVE_OPENMP_OFFLOAD
1526 : : call gpu_copy(ch1c(:,1:nband_me,iband_,ikpt_me), ch1c_tmp, &
1527 : : & int(2,c_size_t)*nband_me)
1528 : : #endif
1529 : : end if
1530 : : end if
1531 : : end do ! iband_
1532 : : end do !idat
1533 : :
1534 [ + + ]: 388532 : if (has_dcwf) then
1535 [ + + + + ]: 384612 : if (ipert==ipert1.and.idir==idir1) then
1536 [ + - ]: 67052 : if(gpu_option==ABI_GPU_DISABLED) then
1537 [ - + + + : 30075190 : gvnlx1 = gvnlx1-gh1
+ + ]
1538 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1539 : : call abi_xaxpy(2*npw1_k*nspinor*ndat, cminusone, &
1540 : 0 : & gh1, 1, gvnlx1, 1, gpu_option=gpu_option)
1541 : : end if
1542 : : end if
1543 : 384612 : do_scprod=0
1544 [ + + ]: 769224 : do idat=1,ndat
1545 [ + + ]: 769224 : if (abs(occ_k(iband+idat-1))>tol8) do_scprod=1
1546 : : end do !idat
1547 [ + + ]: 384612 : if(do_scprod==1) then
1548 : : ! Compute: -<u0_k_i|S^(j1)| Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j>
1549 : : call dotprod_g_batch_full(vdotr,vdoti,istwf_k,npw1_k*nspinor,ndat,2,&
1550 : : & gs1,&
1551 : : & gvnlx1,&
1552 : 316362 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gpu_option=gpu_option)
1553 : : #ifdef HAVE_OPENMP_OFFLOAD
1554 : : !$OMP TARGET UPDATE FROM(vdotr,vdoti) IF(gpu_option==ABI_GPU_OPENMP)
1555 : : #endif
1556 : : ! Add contribution to DDB
1557 : : ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2
1558 : : ! (-1) factor already present
1559 [ + + ]: 632724 : do idat=1,ndat
1560 [ + - ]: 632724 : if (abs(occ_k(iband+idat-1))>tol8) then
1561 : 316362 : d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband+idat-1)*elfd_fact*vdotr(idat)
1562 : 316362 : d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband+idat-1)*elfd_fact*vdoti(idat)
1563 : : end if
1564 : : end do !idat
1565 : : end if
1566 : : end if
1567 : : end if
1568 : :
1569 : :
1570 [ + + ]: 410980 : if (has_dcwf.and.is_metal_or_qne0) then
1571 : :
1572 [ + - ]: 311364 : if(gpu_option==ABI_GPU_DISABLED) then
1573 [ + + + + ]: 1976397606 : gvnlx2 = zero
1574 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1575 : 0 : call gpu_set_to_zero(gvnlx2,int(2,c_size_t)*npw1_k*nspinor*nband_k)
1576 : : end if
1577 : :
1578 : 311364 : jband_me = 0
1579 [ + + ]: 3401280 : do jband=1,nband_k
1580 [ + - ]: 3401280 : if (mpi_enreg%proc_distrb(ikpt,jband,isppol)==me) then
1581 : 3089916 : jband_me = jband_me + 1
1582 : : ! gvnlx1 depends on j only, I have it, and everyone needs it
1583 [ + - ]: 3089916 : if(gpu_option==ABI_GPU_DISABLED) then
1584 [ + + + + ]: 1979176158 : gvnlx2(:,1+(jband-1)*npw1_k*nspinor:jband*npw1_k*nspinor)=cgq(:,1+npw1_k*nspinor*(jband_me-1)+icgq:npw1_k*nspinor*jband_me+icgq)
1585 : : else if(gpu_option==ABI_GPU_OPENMP) then
1586 : : #ifdef HAVE_OPENMP_OFFLOAD
1587 : : call gpu_copy(gvnlx2(:,1+(jband-1)*npw1_k*nspinor:jband*npw1_k*nspinor), &
1588 : : & cgq(:,1+npw1_k*nspinor*(jband_me-1)+icgq:npw1_k*nspinor*jband_me+icgq),&
1589 : : & int(2,c_size_t)*npw1_k*nspinor)
1590 : : #endif
1591 : : end if
1592 : : end if
1593 : : ! xmpi bcast the current jband to other procs in band pool
1594 : : end do
1595 : :
1596 [ - + ]: 311364 : if(xmpi_comm_size(mpi_enreg%comm_band)>1) then
1597 : 0 : call xmpi_sum(gvnlx2, mpi_enreg%comm_band, ierr, use_omp_map=use_ompgpu)
1598 : : end if
1599 : : end if
1600 : :
1601 [ + + ]: 829800 : do idat=1,ndat
1602 : :
1603 : 414900 : bands_treated_now => bands_treated_now_ndat(:,idat)
1604 : :
1605 : : ! If needed, compute here <delta_u^(j1)_k_i|H-Eps_k_i.S|u^(j2)_k_i>
1606 : : ! This is equal to <delta_u^(j1)_k_i|H-Eps_k_i.S|delta_u^(j2)_k_i> (I)
1607 : : ! +<delta_u^(j1)_k_i|H-Eps_k_i.S|u^paral^(j2)_k_i> (II)
1608 : : ! (u^paral^(j2)_k_i is the part of u^(j2)_k_i parallel to active space : metals)
1609 : : ! (I) can be rewritten as:
1610 : : ! Sum_j{ 1/4.<u0_k_i|S^(j1)|u0_k+q_j>.<u0_k+q_j|S^(j2)|u0_k_i>.(Eps_k+q_j-Eps_k_i) }
1611 : : ! (II) can be rewritten as:
1612 : : ! Sum_j{1/2.(occ_kq_j-occ_k_i).Eps1_k,q_ij.<u0_k_i|S^(j1)|u0_k+q_j> }
1613 : : ! where Eps1_k,q_ij=<u0_k+q_j|H^(j2)-1/2(Eps_k+q_j-Eps_k_i)S^(j2)|u0_k_i>
1614 : : ! At first call (when j1=j2), cs1c=<u0_k_i|S^(j1)|u0_k+q_j> is stored
1615 : : ! For the next calls, it is reused.
1616 [ + + ]: 829800 : if (has_dcwf.and.is_metal_or_qne0) then
1617 : : ! dotX is local to my proc, and should accumulate sum over all jband, for my iband_me
1618 : 311364 : dotr=zero;doti=zero
1619 : : ! flag to broadcast the j dependent vector in the band pool, to get full sum over j
1620 : 311364 : do_bcast = 0
1621 : : ! flag to do scalar product: only needed if we are saving cs1c or if the band is occupied
1622 : 311364 : do_scprod = 0
1623 [ + + + + ]: 311364 : if ((ipert==ipert1.and.idir==idir1)) then
1624 : 311364 : do_bcast = 1
1625 : 311364 : do_scprod = 1
1626 : : end if
1627 : 311364 : invocc=zero
1628 [ + + ]: 311364 : if (abs(occ_k(iband+idat-1))>tol8) then
1629 : 243114 : invocc=two/occ_k(iband+idat-1)
1630 : 243114 : do_scprod = 1
1631 : : end if
1632 : : ! does anyone else need the cgq(j) below?
1633 [ + + ]: 3401280 : do iband_ = 1, nband_k
1634 [ + + + + ]: 3401280 : if (bands_treated_now(iband_) > 0 .and. abs(occ_k(iband_))>tol8) then
1635 : 3089916 : do_bcast = 1
1636 : : end if
1637 : : end do
1638 : :
1639 : :
1640 [ + + ]: 311364 : if (do_bcast > 0) then
1641 : : ! Computation of cs1c=<u0_k_i|S^(j1)|u0_k+q_j>
1642 : : ! do _I_ need to calculate the dot1X?
1643 [ + - ]: 257046 : if (do_scprod > 0) then
1644 : : call dotprod_g_batch_half(vdot1r,vdot1i,istwf_k,npw1_k*nspinor,nband_k,2,&
1645 : : & gs1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor),&
1646 : 257046 : & gvnlx2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gpu_option=gpu_option)
1647 : : #ifdef HAVE_OPENMP_OFFLOAD
1648 : : !$OMP TARGET UPDATE FROM(vdot1r,vdot1i) IF(gpu_option==ABI_GPU_OPENMP)
1649 : : #endif
1650 [ + + + + : 257046 : if (ipert==ipert1.and.idir==idir1.and.has_dcwf2) then
+ - ]
1651 [ + + ]: 666168 : cs1c(1,:,iband_me-(ndat-idat),ikpt_me)=vdot1r(:)
1652 [ + + ]: 666168 : cs1c(2,:,iband_me-(ndat-idat),ikpt_me)=vdot1i(:)
1653 : : end if
1654 : : end if
1655 : : end if ! ipert==ipert1.and.idir==idir1 or some iband for some proc in pool is filled
1656 : :
1657 [ + + ]: 311364 : if (do_scprod > 0) then
1658 [ + + ]: 2516526 : do jband=1,nband_k
1659 : : ! Computation of term (I)
1660 [ + + ]: 2259480 : if (has_dcwf2) then
1661 : 2102808 : arg=eig_kq(jband)-eig_k(iband+idat-1)
1662 : 2102808 : dot2r=cs1c(1,jband,iband_me-(ndat-idat),ikpt_me)
1663 : 2102808 : dot2i=cs1c(2,jband,iband_me-(ndat-idat),ikpt_me)
1664 : 2102808 : dotr=dotr+(vdot1r(jband)*dot2r+vdot1i(jband)*dot2i)*arg
1665 : 2102808 : doti=doti+(vdot1i(jband)*dot2r-vdot1r(jband)*dot2i)*arg
1666 : : end if
1667 : : ! Computation of term (II) TODO: the next two ifs could be combined
1668 [ + + ]: 2516526 : if (is_metal) then
1669 [ + + ]: 1734168 : if (abs(rocceig(jband,iband+idat-1))>tol8) then
1670 : 538308 : ii=2*jband-1+(iband+idat-2)*2*nband_k
1671 : 538308 : arg=invocc*rocceig(jband,iband+idat-1)*(eig_k(iband+idat-1)-eig_kq(jband))
1672 : 538308 : dot2r=eig1_k(ii)
1673 : 538308 : dot2i=eig1_k(ii+1)
1674 : 538308 : dotr=dotr+arg*(vdot1r(jband)*dot2r-vdot1i(jband)*dot2i)
1675 : 538308 : doti=doti+arg*(vdot1r(jband)*dot2i+dot2r*vdot1i(jband))
1676 : : end if
1677 : : end if
1678 : : end do
1679 : : end if ! occ bands
1680 : :
1681 : 311364 : dotr=quarter*dotr
1682 : 311364 : doti=quarter*doti
1683 : :
1684 : : ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2))
1685 : : ! Note2: do not sum over bands here - comm_band is a sub communicator of spacecomm, and a full sum is done later
1686 : 311364 : d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband+idat-1)*two*elfd_fact*dotr
1687 : 311364 : d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband+idat-1)*two*elfd_fact*doti
1688 : : end if
1689 : : end do !idat
1690 : :
1691 : : ! Build the matrix element <u0_k_i|H^(j1)-Eps_k_i.S^(j1)|u^(j2)_k,q_i>
1692 : : ! and add contribution to DDB
1693 [ + + ]: 414900 : if (ipert1/=dtset%natom+1) then
1694 : : do_scprod=0
1695 [ + + ]: 640584 : do idat=1,ndat
1696 [ + + ]: 640584 : if (abs(occ_k(iband+idat-1))>tol8) do_scprod=1
1697 : : end do !idat
1698 [ + + ]: 320292 : if(do_scprod==1) then
1699 : : call dotprod_g_batch_full(vdotr,vdoti,istwf_k,npw1_k*nspinor,ndat,2,&
1700 : : & gh1,&
1701 : : & cwavef,&
1702 : 275868 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gpu_option=gpu_option)
1703 : : #ifdef HAVE_OPENMP_OFFLOAD
1704 : : !$OMP TARGET UPDATE FROM(vdotr,vdoti) IF(gpu_option==ABI_GPU_OPENMP)
1705 : : #endif
1706 : : end if
1707 [ + + ]: 640584 : do idat=1,ndat
1708 [ + + ]: 640584 : if (abs(occ_k(iband+idat-1))>tol8) then
1709 : : ! Case ipert1=natom+2 (electric field):
1710 : : ! gh1 contains H^(j1)|u0_k_i> (VHxc constant) which corresponds
1711 : : ! to i.d/dk in Eq. (38) of Gonze, PRB 55, 10355 (1997) [[cite:Gonze1997a]].
1712 : : ! * if ipert==natom+2, we apply directly Eq. (38)
1713 : : ! * if ipert/=natom+2, Born effective charges are minus D2E
1714 : 275868 : d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband+idat-1)*two*elfd_fact*vdotr(idat)
1715 : 275868 : d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband+idat-1)*two*elfd_fact*vdoti(idat)
1716 : : end if
1717 : : end do !idat
1718 : :
1719 : : ! Or compute localisation tensor (ddk)
1720 : : ! See M. Veithen thesis Eq(2.5)
1721 : : ! MT jan-2010: this is probably not correctly implemented for PAW !!!
1722 : : ! missing terms due to S^(1) and S^(2)
1723 : : else
1724 : : ! note: gh1 used as temporary space (to store idir ddk WF)
1725 [ + + ]: 94608 : if (idir==idir1) then
1726 [ + - ]: 31536 : if(gpu_option==ABI_GPU_DISABLED) then
1727 [ + - + + : 33145236 : gvnlx1(:,1:ndat*npw1_k*nspinor)=cwavef(:,1:ndat*npw1_k*nspinor)
+ + + + +
+ ]
1728 : : else if(gpu_option==ABI_GPU_OPENMP) then
1729 : : #ifdef HAVE_OPENMP_OFFLOAD
1730 : : call gpu_copy(gvnlx1, cwavef, int(2,c_size_t)*npw1_k*nspinor*ndat)
1731 : : #endif
1732 : : end if
1733 : : else
1734 [ - + - - ]: 63072 : if (need_ddk_file.and.ddkfil(idir1)/=0) then
1735 : : !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt)
1736 : : !ik_ddk = indkpt1(ikpt)
1737 : : !call ddks(idir1)%read_bks(iband, ik_ddk, isppol, xmpio_single, cg_bks=gvnlx1_tmp)
1738 [ # # # # ]: 0 : gvnlx1(:,1:ndat*npw1_k*nspinor) = cg_ddk(:,1+(iband_me-ndat)*npw1_k*nspinor:iband_me*npw1_k*nspinor,idir1)
1739 : : #ifdef HAVE_OPENMP_OFFLOAD
1740 : : !$OMP TARGET UPDATE TO(gvnlx1) IF(gpu_option==ABI_GPU_OPENMP)
1741 : : #endif
1742 : : else
1743 [ + - ]: 63072 : if(gpu_option==ABI_GPU_DISABLED) then
1744 [ + + + + ]: 33176772 : gvnlx1=zero
1745 [ # # ]: 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1746 : 0 : call gpu_set_to_zero(gvnlx1,int(2,c_size_t)*npw1_k*nspinor*ndat)
1747 : : end if
1748 : : end if
1749 : : end if
1750 : 94608 : do_scprod=0
1751 [ + + ]: 189216 : do idat=1,ndat
1752 [ + + ]: 189216 : if (abs(occ_k(iband+idat-1))>tol8) then
1753 : 63522 : do_scprod=1
1754 : : end if
1755 : : end do
1756 [ + + ]: 94608 : if (do_scprod==1) then
1757 : : call dotprod_g_batch_full(vdotr,vdoti,istwf_k,npw1_k*nspinor,ndat,2,&
1758 : : & gvnlx1,&
1759 : : & cwavef,&
1760 : 63522 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gpu_option=gpu_option)
1761 : : call dotprod_g_batch_full(vdot1r,vdot1i,istwf_k,npw1_k*nspinor,ndat,2,&
1762 : : & cwave0,&
1763 : : & gvnlx1,&
1764 : 63522 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gpu_option=gpu_option)
1765 : : call dotprod_g_batch_full(vdot2r,vdot2i,istwf_k,npw1_k*nspinor,ndat,2,&
1766 : : & cwavef,&
1767 : : & cwave0,&
1768 : 63522 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gpu_option=gpu_option)
1769 : : #ifdef HAVE_OPENMP_OFFLOAD
1770 : : !$OMP TARGET UPDATE FROM(vdotr,vdoti,vdot1r,vdot1i,vdot2r,vdot2i) IF(gpu_option==ABI_GPU_OPENMP)
1771 : : #endif
1772 [ + + ]: 127044 : do idat=1,ndat
1773 [ + - ]: 127044 : if (abs(occ_k(iband+idat-1))>tol8) then
1774 : 63522 : vdotr(idat)=vdotr(idat)-(vdot1r(idat)*vdot2r(idat)-vdot1i(idat)*vdot2i(idat))
1775 : 63522 : vdoti(idat)=vdoti(idat)-(vdot1r(idat)*vdot2i(idat)+vdot1i(idat)*vdot2r(idat))
1776 : 63522 : d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband+idat-1)*vdotr(idat)/(nband_kocc*two)
1777 : 63522 : d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband+idat-1)*vdoti(idat)/(nband_kocc*two)
1778 : : end if
1779 : : end do !idat
1780 : : end if
1781 : : end if
1782 : :
1783 : : ! Accumulate here 1st-order density change due to overlap operator changes (if any)
1784 [ + + ]: 414900 : if (has_drho) then
1785 : : ! Compute here delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
1786 : : ! (see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42))
1787 [ - + + - : 870012 : ABI_MALLOC(dcwavef,(2,npw1_k*nspinor*ndat))
+ - + - -
+ - + ]
1788 : : #ifdef HAVE_OPENMP_OFFLOAD
1789 : : !$OMP TARGET ENTER DATA MAP(alloc:dcwavef) IF(gpu_option==ABI_GPU_OPENMP)
1790 : : #endif
1791 [ - + + - : 2178024 : ABI_MALLOC(dcwaveprj,(dtset%natom,nspinor*ndat))
- + + - +
- + - - +
- + - + +
+ + + ]
1792 : 290004 : call pawcprj_alloc(dcwaveprj,0,gs_hamkq%dimcprj)
1793 : : ! NB: have to call getdc with all band processors to distribute cgq cprjq correctly
1794 : : call getdc1(iband,band_procs,bands_treated_now_ndat,cgq,cprjq,&
1795 : : & dcwavef,&
1796 : : & dcwaveprj,&
1797 : : & ibgq,icgq,istwf_k,mcgq,&
1798 : : & mcprjq,mpi_enreg,ndat,dtset%natom,nband_k,nband_me,npw1_k,nspinor,1,&
1799 : 290004 : & gs1,gpu_option=gpu_option)
1800 : :
1801 : 290004 : option=1;wfcorr=0
1802 [ + + ]: 290004 : if (abs(occ_k(iband))>tol8) then
1803 : : ! Accumulate 1st-order density due to delta_u^(j1)
1804 : : call dfpt_accrho(cplex,cwave0,&
1805 : : & dcwavef,&
1806 : : & dcwavef,&
1807 : : & cwaveprj0_idir1,&
1808 : : & dcwaveprj,&
1809 : : & eig_k(iband:iband+ndat-1),gs_hamkq,iband,idir1,ipert1,isppol,dtset%kptopt,&
1810 : : & mpi_enreg,ndat,dtset%natom,nband_k,1,npw_k,npw1_k,nspinor,occ_k,option,&
1811 [ - + ]: 252840 : & pawdrhoij1_unsym(:,idir1),drhoaug1(:,:,:,idir1),tim_fourwf,wfcorr,wtk_k)
1812 : : end if
1813 : :
1814 : 290004 : call pawcprj_free(dcwaveprj)
1815 [ + - + + : 1018008 : ABI_FREE(dcwaveprj)
- + - + -
+ ]
1816 : : #ifdef HAVE_OPENMP_OFFLOAD
1817 : : !$OMP TARGET EXIT DATA MAP(delete:dcwavef) IF(gpu_option==ABI_GPU_OPENMP)
1818 : : #endif
1819 [ - + ]: 290004 : ABI_FREE(dcwavef)
1820 : : end if ! has_drho
1821 : : !do idat=1,ndat
1822 : : ! !call pawcprj_output(cwaveprj0_idir1(:,1+(idat-1)*nspinor:idat*nspinor),prtgrads=1)
1823 : : !end do
1824 : :
1825 [ - + + + ]: 555180 : if((usecprj==1).and..not.(associated(cwaveprj0_idir1,cwaveprj0)))then
1826 : 317808 : call pawcprj_free(cwaveprj0_idir1)
1827 [ - + ]: 317808 : ABI_FREE(cwaveprj0_idir1)
1828 : : end if
1829 : :
1830 : : ! End of loops
1831 : : end do ! idir1
1832 : : end do ! iband
1833 : :
1834 [ - + ]: 23225 : ABI_FREE(band_procs)
1835 : :
1836 : : ! Accumulate contribution of this k-point
1837 [ + + + + ]: 232250 : d2nl (:,:,ipert1,idir,ipert)=d2nl (:,:,ipert1,idir,ipert)+d2nl_k (:,:)
1838 [ + - + + : 232250 : if (usepaw==1) d2ovl(:,:,ipert1,idir,ipert)=d2ovl(:,:,ipert1,idir,ipert)+d2ovl_k(:,:)
+ + ]
1839 : :
1840 : :
1841 : : ! Deallocations of arrays used for this k-point
1842 : : #ifdef HAVE_OPENMP_OFFLOAD
1843 : : !$OMP TARGET EXIT DATA MAP(delete:gvnlx1,gvnlx1_tmp) IF(gpu_option==ABI_GPU_OPENMP)
1844 : :
1845 : : !$OMP TARGET EXIT DATA MAP(delete:gh1) IF(gpu_option==ABI_GPU_OPENMP)
1846 : : !$OMP TARGET EXIT DATA MAP(delete:gs1) IF(has_dcwf .and. gpu_option==ABI_GPU_OPENMP)
1847 : : #endif
1848 [ + + ]: 23225 : if (has_dcwf.and.is_metal_or_qne0) then
1849 : : #ifdef HAVE_OPENMP_OFFLOAD
1850 : : !$OMP TARGET EXIT DATA MAP(delete:gvnlx2) IF(gpu_option==ABI_GPU_OPENMP)
1851 : : #endif
1852 [ - + ]: 16304 : ABI_FREE(gvnlx2)
1853 : : end if
1854 : :
1855 [ - + ]: 23225 : ABI_FREE(gvnlx1)
1856 [ - + ]: 23225 : ABI_FREE(gvnlx1_tmp)
1857 [ - + ]: 23225 : ABI_FREE(gh1)
1858 [ - + ]: 23225 : ABI_FREE(gs1)
1859 : : if (need_wfk) then
1860 : 23225 : nullify(cwave0)
1861 : : end if
1862 : : if (need_wf1) then
1863 : 23225 : nullify(cwavef)
1864 : : end if
1865 [ - + ]: 23225 : ABI_FREE(kg_k)
1866 [ - + ]: 23225 : ABI_FREE(kg1_k)
1867 [ - + ]: 23225 : ABI_FREE(ylm_k)
1868 [ - + ]: 23225 : ABI_FREE(ylm1_k)
1869 [ - + ]: 23225 : ABI_FREE(ylmgr1_k)
1870 [ - + ]: 23225 : ABI_FREE(kpg_k)
1871 [ - + ]: 23225 : ABI_FREE(kpg1_k)
1872 [ - + ]: 23225 : ABI_FREE(d2nl_k)
1873 [ - + ]: 23225 : ABI_FREE(d2ovl_k)
1874 [ - + ]: 23225 : ABI_FREE(eig_k)
1875 [ - + ]: 23225 : ABI_FREE(eig_kq)
1876 [ - + ]: 23225 : ABI_FREE(eig1_k)
1877 [ - + ]: 23225 : ABI_FREE(occ_k)
1878 : : #ifdef HAVE_OPENMP_OFFLOAD
1879 : : !$OMP TARGET EXIT DATA MAP(delete:vdotr,vdoti,vdot1r,vdot1i,vdot2r,vdot2i) IF(gpu_option==ABI_GPU_OPENMP)
1880 : : #endif
1881 [ - + ]: 23225 : ABI_FREE(vdotr)
1882 [ - + ]: 23225 : ABI_FREE(vdoti)
1883 [ - + ]: 23225 : ABI_FREE(vdot1r)
1884 [ - + ]: 23225 : ABI_FREE(vdot1i)
1885 [ - + ]: 23225 : ABI_FREE(vdot2r)
1886 [ - + ]: 23225 : ABI_FREE(vdot2i)
1887 [ + + ]: 23225 : if (is_metal) then
1888 [ - + ]: 5988 : ABI_FREE(doccde_k)
1889 [ - + ]: 5988 : ABI_FREE(doccde_kq)
1890 [ - + ]: 5988 : ABI_FREE(occ_kq)
1891 [ - + ]: 5988 : ABI_FREE(rocceig)
1892 : : end if
1893 [ - + ]: 23225 : ABI_FREE(dkinpw)
1894 [ - + ]: 23225 : ABI_FREE(kinpw1)
1895 : : #ifdef HAVE_OPENMP_OFFLOAD
1896 : : if (ipert1<=dtset%natom) then
1897 : : !$OMP TARGET EXIT DATA MAP(delete:ffnlk) IF(gpu_option==ABI_GPU_OPENMP)
1898 : : end if
1899 : : !$OMP TARGET EXIT DATA MAP(delete:ffnl1) IF(gpu_option==ABI_GPU_OPENMP)
1900 : : !$OMP TARGET EXIT DATA MAP(delete:ph3d) IF(gpu_option==ABI_GPU_OPENMP)
1901 : : #endif
1902 [ - + ]: 23225 : ABI_FREE(ph3d)
1903 [ + + ]: 23225 : if (allocated(ph3d1)) then
1904 : 11360 : ABI_FREE(ph3d1)
1905 : : end if
1906 [ - + ]: 23225 : ABI_FREE(ffnlk)
1907 [ - + ]: 23225 : ABI_FREE(ffnl1)
1908 [ + + ]: 23225 : if (ipert1>dtset%natom) then
1909 : : #ifdef HAVE_OPENMP_OFFLOAD
1910 : : !$OMP TARGET EXIT DATA MAP(delete:ffnl1_idir1) IF(gpu_option==ABI_GPU_OPENMP)
1911 : : #endif
1912 [ - + ]: 7391 : ABI_FREE(ffnl1_idir1)
1913 : : end if
1914 : 23225 : nullify(ffnl1_idir1)
1915 [ + - ]: 23225 : if (usecprj==1) then
1916 : 23225 : call pawcprj_free(cwaveprj0)
1917 [ + - + + : 75160 : ABI_FREE(cwaveprj0)
- + - + -
+ ]
1918 : : end if
1919 : 23225 : nullify(cwaveprj0_idir1)
1920 : : ! Shift arrays
1921 : 23225 : bdtot_index=bdtot_index+nband_k
1922 : 23225 : bd2tot_index=bd2tot_index+2*nband_k**2
1923 [ + - ]: 23225 : if (mkmem/=0) then
1924 : 23225 : ibg=ibg+nspinor*nband_me
1925 : 23225 : icg=icg+npw_k*nspinor*nband_me
1926 : 23225 : ikg=ikg+npw_k
1927 : : end if
1928 [ + - ]: 23225 : if (mkqmem/=0) then
1929 : 23225 : ibgq=ibgq+nspinor*nband_me
1930 : 23225 : icgq=icgq+npw1_k*nspinor*nband_me
1931 : : end if
1932 [ + - ]: 24531 : if (mk1mem/=0) then
1933 : 23225 : ibg1=ibg1+nspinor*nband_me
1934 : 23225 : icg1=icg1+npw1_k*nspinor*nband_me
1935 : 23225 : ikg1=ikg1+npw1_k
1936 : : end if
1937 : :
1938 : :
1939 : : end do ! End loop over K-POINTS
1940 : : !----------------------------------------------------------------
1941 : :
1942 : : ! Transfer 1st-order density change due to overlap; also take into account the spin.
1943 [ + + ]: 2531 : if(has_drho) then
1944 [ + + ]: 4008 : do kdir1=1,mdir1
1945 : 3006 : idir1=jdir1(kdir1)
1946 : : call fftpac(isppol,mpi_enreg,nspden,cplex*dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),&
1947 : : & cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),&
1948 : 4008 : & dtset%ngfft,drho1wfr(:,:,idir1),drhoaug1(:,:,:,idir1),1)
1949 : : end do
1950 : : end if
1951 : :
1952 : : end do ! End loop over SPINS
1953 : : !----------------------------------------------------------------
1954 : :
1955 : : ! Free memory used for this type of perturbation
1956 : 1225 : call rf_hamkq%free()
1957 [ + + ]: 1225 : if (has_vectornd) then
1958 [ - + ]: 57 : ABI_FREE(vectornd_pac_idir)
1959 : : end if
1960 [ + + ]: 1225 : if (allocated(vlocal)) then
1961 : 57 : ABI_FREE(vlocal)
1962 : : end if
1963 [ + + ]: 1225 : if (allocated(vtrial_)) then
1964 : 57 : ABI_FREE(vtrial_)
1965 : : end if
1966 [ + + ]: 1225 : if (has_drho) then
1967 [ - + ]: 950 : ABI_FREE(drhoaug1)
1968 : : end if
1969 [ + + ]: 1225 : if (need_pawij10) then
1970 [ + + ]: 4336 : do kdir1=1,mdir1
1971 : 3252 : idir1=jdir1(kdir1)
1972 : 4336 : call paw_ij_free(paw_ij10(:,idir1))
1973 : : end do
1974 [ - + ]: 1084 : ABI_FREE(e1kbfr_spin)
1975 : : end if
1976 [ + - + + : 8134 : ABI_FREE(paw_ij10)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
1977 : :
1978 : : ! In case of parallelism, sum 1st-order density and occupation matrix over processors
1979 [ + + ]: 1225 : if (has_drho.and.xmpi_paral==1) then
1980 : :
1981 : : ! Accumulate 1st-order density
1982 : 950 : call timab(48,1,tsec)
1983 : 950 : bufsz=cplex*dtset%nfft*nspden*mdir1
1984 [ - + + - : 2850 : ABI_MALLOC(buffer,(bufsz))
+ - + - -
+ - + -
+ ]
1985 [ + + ]: 1900 : buffer(1:bufsz)=reshape(drho1wfr,(/bufsz/))
1986 : 950 : call xmpi_sum(buffer,bufsz,spaceworld,ierr)
1987 [ + + ]: 3800 : drho1wfr(:,:,:)=reshape(buffer(1:bufsz),(/cplex*dtset%nfft,nspden,mdir1/))
1988 [ - + ]: 950 : ABI_FREE(buffer)
1989 : 950 : call timab(48,2,tsec)
1990 : :
1991 : : ! Accumulate 1st-order PAW occupancies
1992 [ + - ]: 1900 : if (usepaw==1) then
1993 : 950 : call pawrhoij_mpisum_unpacked(pawdrhoij1_unsym,spaceworld)
1994 : : end if
1995 : :
1996 : : end if
1997 : :
1998 : : ! Compute second part of overlap contribution (due to VHxc^(j2)(tild_n+hat_n))
1999 : 361 : if (has_drho) then
2000 : :
2001 [ - + + - : 3800 : ABI_MALLOC(drhor1,(cplex*nfftf,nspden))
- + + - +
- + - - +
- + ]
2002 [ - + + - : 2850 : ABI_MALLOC(dnhat1,(cplex*nfftf,nspden))
- + + - -
+ - + ]
2003 : :
2004 : : ! LOOP OVER PERTURBATION DIRECTIONS
2005 [ + + ]: 3800 : do kdir1=1,mdir1
2006 : 2850 : idir1=jdir1(kdir1)
2007 : :
2008 : : ! Build and symmetrize 1st-order density change due to change of overlap
2009 [ - + + - : 8550 : ABI_MALLOC(drho1wfg,(2,dtset%nfft))
+ - + - -
+ - + ]
2010 : : call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,&
2011 : : & nspden,nsppol,nsym1,phnons1,drho1wfg,drho1wfr(:,:,idir1),&
2012 : 2850 : & rprimd,symaf1,symrl1,tnons1)
2013 [ + - ]: 2850 : if (dtset%pawstgylm/=0) then
2014 : 2850 : option=0
2015 : : call pawnhatfr(option,idir1,ipert1,my_natom,dtset%natom,nspden,dtset%ntypat,&
2016 : : & pawang,pawfgrtab,pawrhoij,pawtab,rprimd,&
2017 [ - + - - : 2850 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
- - - + -
- - - - +
- - - - -
- - - - -
- - - - -
- - - -
+ ]
2018 : : end if
2019 : : call pawmkrho(1,arg,cplex,gprimd,idir1,indsy1,ipert1,&
2020 : : & mpi_enreg,my_natom,dtset%natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,&
2021 : : & pawfgr,pawfgrtab,-10001,pawdrhoij1(:,idir1),pawdrhoij1_unsym(:,idir1),pawtab,&
2022 : : & dtset%qptn,drho1wfg,drho1wfr(:,:,idir1),drhor1,rprimd,symaf1,symrc1,dtset%typat,&
2023 [ - + - - : 2850 : & ucvol,dtset%usewvl,xred,pawang_sym=pawang1,pawnhat=dnhat1,pawrhoij0=pawrhoij)
- - - + -
- - - - +
- - - - -
- - - - -
- - - - -
- - - -
+ ]
2024 [ - + ]: 2850 : ABI_FREE(drho1wfg)
2025 : :
2026 : : ! Compute plane-wave contribution to overlap contribution
2027 : : ! This is subtle as it is a mix of Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
2028 : : ! Details:
2029 : : ! The VH(tild_nZc)^(1) term of Eq(79) is:
2030 : : ! <VH(tild_nZc)^(j2)|delta_tild_rho^(j1)> = <vpsp1|drhor1-dnhat1>
2031 : : ! The first term of Eq(80) is:
2032 : : ! <VHxc^(j2)|delta_tild_rho^(j1)+delta_hat_rho^(j1)> = <vtrial1-vpsp1|drhor1>
2033 : : ! The addition of these two terms gives:
2034 : : ! <vtrial1|drhor1>-<vpsp1|dnhat1>
2035 : : ! And this is more subtle when usexcnhat=0
2036 : 2850 : call dotprod_vn(cplex,drhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vtrial1,ucvol)
2037 [ + + ]: 2850 : if (usexcnhat/=0) then
2038 : 2043 : call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,1 ,2,vpsp1,ucvol)
2039 : : else
2040 [ - + + - : 2421 : ABI_MALLOC(vtmp1,(cplex*nfftf,nspden))
- + + - -
+ - + ]
2041 [ + + ]: 1740 : do ispden=1,nspden
2042 [ + + ]: 38357748 : vtmp1(:,ispden)=vtrial1(:,ispden)-vhartr1(:)
2043 : : end do
2044 : 807 : call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,nspden,2,vtmp1,ucvol)
2045 [ - + ]: 807 : ABI_FREE(vtmp1)
2046 : : end if
2047 : 2850 : dotr=dot1r-dot2r;doti=dot1i-dot2i
2048 : :
2049 : : ! Compute on-site contributions to overlap contribution
2050 : : ! (two last terms of Eq(80) of PRB 78, 035105 (2008)) [[cite:Audouze2008]]
2051 : : ! (note: Dij^(j2) and Vxc^(j2) are computed for ipert at first call)
2052 : : call pawdfptenergy(epawnst,ipert,ipert1,dtset%ixc,my_natom,dtset%natom,dtset%ntypat,&
2053 : : & nzlmopt_ipert,nzlmopt_ipert1,paw_an,paw_an1,paw_ij1,pawang,dtset%pawprtvol,&
2054 : : & pawrad,pawrhoij1,pawdrhoij1(:,idir1),pawtab,dtset%pawxcdev,dtset%xclevel,&
2055 [ - + - - : 2850 : & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
- - - + -
- - - - +
- - - - -
+ - - - -
- + - + -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
+ - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + ]
2056 : :
2057 : : ! Accumulate in 2nd-order matrix:
2058 : : ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2
2059 : : ! has to take the complex conjugate because we want here Int[VHxc^(j1)^*.delta_rho^(j2)]
2060 : 2850 : dotr=dotr+epawnst(1);doti=-(doti+epawnst(2))
2061 : 2850 : d2ovl_drho(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr
2062 : 3800 : d2ovl_drho(2,idir1,ipert1,idir,ipert)=elfd_fact*doti
2063 : :
2064 : : !if(kpert1==5) then
2065 : : ! ABI_BUG("toto")
2066 : : !end if
2067 : : end do ! End loop over perturbation directions
2068 : :
2069 : : ! Free no more needed memory
2070 [ - + ]: 950 : ABI_FREE(drhor1)
2071 [ - + ]: 950 : ABI_FREE(dnhat1)
2072 [ - + ]: 950 : ABI_FREE(drho1wfr)
2073 [ + + ]: 2995 : do iatom=1,my_natom
2074 [ + - ]: 2045 : if (pawfgrtab(iatom)%nhatfr_allocated>0) then
2075 [ - + ]: 2045 : ABI_FREE(pawfgrtab(iatom)%nhatfr)
2076 : : end if
2077 : 2995 : pawfgrtab(iatom)%nhatfr_allocated=0
2078 : : end do
2079 [ + + ]: 950 : if (paral_atom) then
2080 [ + + ]: 2192 : do kdir1=1,mdir1
2081 : 1644 : idir1=jdir1(kdir1)
2082 : 2192 : call pawrhoij_free(pawdrhoij1_unsym(:,idir1))
2083 : : end do
2084 [ - + ]: 548 : ABI_FREE(pawdrhoij1_unsym)
2085 : : end if
2086 [ + + ]: 3800 : do kdir1=1,mdir1
2087 : 2850 : idir1=jdir1(kdir1)
2088 : 3800 : call pawrhoij_free(pawdrhoij1(:,idir1))
2089 : : end do
2090 [ + - + + : 7085 : ABI_FREE(pawdrhoij1)
- + - + -
+ - + - +
- + - + ]
2091 : : end if ! has_drho
2092 : :
2093 : : ! End loop over perturbations (j1)
2094 : : end do
2095 : :
2096 : : #ifdef HAVE_OPENMP_OFFLOAD
2097 : : !$OMP TARGET EXIT DATA MAP(delete:ch1c,ch1c_tmp) IF(gpu_option==ABI_GPU_OPENMP)
2098 : : #endif
2099 : : !Final deallocations
2100 [ - + ]: 361 : ABI_FREE(ch1c)
2101 [ - + ]: 361 : ABI_FREE(ch1c_tmp)
2102 [ + - + + ]: 361 : if (usepaw==1.and.is_metal_or_qne0) then
2103 [ - + ]: 252 : ABI_FREE(cs1c)
2104 [ - + ]: 252 : ABI_FREE(cs1c_tmp)
2105 : : end if
2106 : 361 : call gs_hamkq%free()
2107 [ + + ]: 361 : if(allocated(vectornd_pac)) then
2108 : 57 : ABI_FREE(vectornd_pac)
2109 : : end if
2110 [ + + ]: 361 : if(allocated(vxctaulocal)) then
2111 : 6 : ABI_FREE(vxctaulocal)
2112 : : end if
2113 : :
2114 : : !In case of parallelism, sum over processors
2115 : : if (xmpi_paral==1)then
2116 : 361 : call timab(161,1,tsec)
2117 : 361 : call xmpi_barrier(spaceworld)
2118 : 361 : call timab(161,2,tsec)
2119 [ + - - + : 1083 : ABI_MALLOC(buffer,(6*mpert*(1+usepaw)))
+ - + - +
- - + - +
- + ]
2120 [ + + ]: 722 : buffer(1:6*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/6*mpert/))
2121 [ + - + + ]: 722 : if (usepaw==1) buffer(6*mpert+1:6*mpert+6*mpert)=reshape(d2ovl(:,:,:,idir,ipert),(/6*mpert/))
2122 : 361 : call timab(48,1,tsec)
2123 : 361 : call xmpi_sum(buffer,6*mpert*(1+usepaw),spaceworld,ierr)
2124 : 361 : call timab(48,2,tsec)
2125 [ + + ]: 1444 : d2nl (:,:,:,idir,ipert)=reshape(buffer(1:6*mpert),(/2,3,mpert/))
2126 [ + - + + ]: 1444 : if (usepaw==1) d2ovl(:,:,:,idir,ipert)=reshape(buffer(6*mpert+1:6*mpert+6*mpert),(/2,3,mpert/))
2127 [ - + ]: 361 : ABI_FREE(buffer)
2128 : : end if
2129 : :
2130 : : !Build complete d2ovl matrix
2131 [ + - + + : 36081 : if (usepaw==1) d2ovl(:,:,:,idir,ipert)=d2ovl(:,:,:,idir,ipert)+d2ovl_drho(:,:,:,idir,ipert)
+ + + + ]
2132 : :
2133 [ + - ]: 361 : if (usepaw==1) then
2134 [ - + ]: 361 : ABI_FREE(d2ovl_drho)
2135 : : end if
2136 : :
2137 : : !Close the ddk WF files
2138 [ + + ]: 361 : if (has_ddk_file) then
2139 [ + + ]: 536 : do kdir1=1,mdir1
2140 : 402 : idir1=jdir1(kdir1)
2141 [ + - ]: 536 : if (ddkfil(idir1)/=0) call ddks(idir1)%close()
2142 : : end do
2143 [ - + ]: 134 : ABI_FREE(cg_ddk)
2144 : : end if
2145 [ - + ]: 361 : ABI_FREE(jpert1)
2146 : 361 : ABI_FREE(jdir1)
2147 : :
2148 : : !Symmetrize the phonons contributions, as was needed for the forces in a GS calculation
2149 [ - + + - : 1083 : ABI_MALLOC(work,(2,3,dtset%natom))
+ - + - -
+ - + ]
2150 [ + + ]: 1406 : do ipert1=1,dtset%natom
2151 [ + + ]: 4541 : do idir1=1,3
2152 [ + + ]: 10450 : work(:,idir1,ipert1)=d2nl(:,idir1,ipert1,idir,ipert)
2153 : : end do
2154 : : end do
2155 : 361 : call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1)
2156 [ + - ]: 361 : if (usepaw==1) then
2157 [ + + ]: 1406 : do ipert1=1,dtset%natom
2158 [ + + ]: 4541 : do idir1=1,3
2159 [ + + ]: 10450 : work(:,idir1,ipert1)=d2ovl(:,idir1,ipert1,idir,ipert)
2160 : : end do
2161 : : end do
2162 : 361 : call dfpt_sygra(dtset%natom,d2ovl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1)
2163 : : end if
2164 [ - + ]: 361 : ABI_FREE(work)
2165 : :
2166 : : !In the case of the strain perturbation time-reversal symmetry will always
2167 : : !be true so imaginary part of d2nl will be must be set to zero here since
2168 : : !the symmetry-reduced kpt set will leave a non-zero imaginary part.
2169 [ + + + + ]: 361 : if(ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
2170 [ + + + + ]: 3096 : d2nl(2,:,:,idir,ipert)=zero
2171 [ + - + + : 3096 : if (usepaw==1) d2ovl(2,:,:,idir,ipert)=zero
+ + ]
2172 : : else
2173 [ + + + + ]: 2601 : d2nl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero
2174 [ + - + + : 2601 : if (usepaw==1) d2ovl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero
+ + ]
2175 : : end if
2176 : :
2177 : :
2178 : : !Symmetrize the strain perturbation contributions, as was needed for the stresses in a GS calculation
2179 : : !if (ipert==dtset%natom+3.or.ipert==dtset%natom+4)then
2180 [ + + ]: 361 : if (nsym1>1) then
2181 [ - + ]: 124 : ABI_MALLOC(work,(6,1,1))
2182 : 124 : ii=0
2183 [ + + ]: 372 : do ipert1=dtset%natom+3,dtset%natom+4
2184 [ + + ]: 1116 : do idir1=1,3
2185 : 744 : ii=ii+1
2186 : 992 : work(ii,1,1)=d2nl(1,idir1,ipert1,idir,ipert)
2187 : : end do
2188 : : end do
2189 : 124 : call stresssym(gprimd,nsym1,work(:,1,1),symrc1)
2190 : 124 : ii=0
2191 [ + + ]: 372 : do ipert1=dtset%natom+3,dtset%natom+4
2192 [ + + ]: 1116 : do idir1=1,3
2193 : 744 : ii=ii+1
2194 : 992 : d2nl(1,idir1,ipert1,idir,ipert)=work(ii,1,1)
2195 : : end do
2196 : : end do
2197 [ + - ]: 124 : if (usepaw==1) then
2198 : 124 : ii=0
2199 [ + + ]: 372 : do ipert1=dtset%natom+3,dtset%natom+4
2200 [ + + ]: 1116 : do idir1=1,3
2201 : 744 : ii=ii+1
2202 : 992 : work(ii,1,1)=d2ovl(1,idir1,ipert1,idir,ipert)
2203 : : end do
2204 : : end do
2205 : 124 : call stresssym(gprimd,nsym1,work(:,1,1),symrc1)
2206 : 124 : ii=0
2207 [ + + ]: 372 : do ipert1=dtset%natom+3,dtset%natom+4
2208 [ + + ]: 1116 : do idir1=1,3
2209 : 744 : ii=ii+1
2210 : 992 : d2ovl(1,idir1,ipert1,idir,ipert)=work(ii,1,1)
2211 : : end do
2212 : : end do
2213 [ - + ]: 124 : ABI_FREE(work)
2214 : : end if
2215 : : end if
2216 : : !end if
2217 : :
2218 : : !Must also symmetrize the electric field perturbation response !
2219 : : !Note: d2ovl is not symetrized because it is zero for electric field perturbation
2220 [ + + ]: 361 : if (has_ddk_file) then
2221 [ - + ]: 134 : ABI_MALLOC(d2nl_elfd,(2,3))
2222 : : ! There should not be any imaginary part, but stay general (for debugging)
2223 [ + + + + ]: 1340 : d2nl_elfd (:,:)=d2nl(:,:,dtset%natom+2,idir,ipert)
2224 [ + + ]: 536 : do ii=1,3
2225 : 402 : sumelfd(:)=zero
2226 [ + + ]: 1572 : do ia=1,nsym1
2227 [ + + ]: 4680 : do jj=1,3
2228 [ + + ]: 4680 : if(symrl1(ii,jj,ia)/=0)then
2229 [ - + ]: 1306 : if(ddkfil(jj)==0)then
2230 : 0 : blkflg(ii,dtset%natom+2,idir,ipert)=0
2231 : : end if
2232 : : end if
2233 : : end do
2234 : 1170 : symfact(1)=dble(symrl1(ii,1,ia))
2235 : 1170 : symfact(2)=dble(symrl1(ii,2,ia))
2236 : 1170 : symfact(3)=dble(symrl1(ii,3,ia))
2237 : : sumelfd(:)=sumelfd(:)+symfact(1)*d2nl_elfd(:,1) &
2238 [ + + ]: 3912 : & +symfact(2)*d2nl_elfd(:,2)+symfact(3)*d2nl_elfd(:,3)
2239 : : end do
2240 [ + + ]: 1340 : d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1)
2241 : : end do
2242 [ - + ]: 134 : ABI_FREE(d2nl_elfd)
2243 : : end if
2244 : :
2245 : : !Overlap: store the diagonal part of the matrix in the
2246 : : ! 2nd-order energy non-stationnary expression
2247 [ + - ]: 361 : eovl1=zero;if (usepaw==1) eovl1=d2ovl(1,idir,ipert,idir,ipert)
2248 : :
2249 [ - + ]: 361 : ABI_FREE(bands_treated_now_ndat)
2250 : :
2251 : : #ifdef HAVE_OPENMP_OFFLOAD
2252 : : !$OMP TARGET EXIT DATA MAP(delete:cg,cg1) IF(gpu_option==ABI_GPU_OPENMP)
2253 : : #endif
2254 : :
2255 : 361 : call destroy_mpi_enreg(mpi_enreg_seq)
2256 : 361 : call timab(566,2,tsec)
2257 : :
2258 : : ABI_NVTX_END_RANGE()
2259 : :
2260 : : DBG_EXIT("COLL")
2261 : :
2262 [ - + - + : 1444 : end subroutine dfpt_nstpaw
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - + - +
- - - - -
- - + - +
- + ]
2263 : : !!***
2264 : :
2265 : : !!****f* ABINIT/dfpt_nstwf
2266 : : !! NAME
2267 : : !! dfpt_nstwf
2268 : : !!
2269 : : !! FUNCTION
2270 : : !! This routine computes the non-local contribution to the
2271 : : !! 2DTE matrix elements, in the non-stationary formulation
2272 : : !! Only for norm-conserving pseudopotentials (no PAW)
2273 : : !!
2274 : : !! COPYRIGHT
2275 : : !! Copyright (C) 1999-2026 ABINIT group (XG,AR,MB,MVer,MT, MVeithen)
2276 : : !! This file is distributed under the terms of the
2277 : : !! GNU General Public License, see ~abinit/COPYING
2278 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2279 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
2280 : : !!
2281 : : !! INPUTS
2282 : : !! cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions at k
2283 : : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
2284 : : !! ddkfil(3)=unit numbers for the three possible ddk files for ipert1
2285 : : !! equal to 0 if no dot file is available for this direction
2286 : : !! dtset <type(dataset_type)>=all input variables for this dataset
2287 : : !! eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
2288 : : !! eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
2289 : : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
2290 : : !! icg=shift to be applied on the location of data in the array cg
2291 : : !! icg1=shift to be applied on the location of data in the array cg1
2292 : : !! idir=direction of the current perturbation
2293 : : !! ikpt=number of the k-point
2294 : : !! ipert=type of the perturbation
2295 : : !! isppol=1 for unpolarized, 2 for spin-polarized
2296 : : !! istwf_k=parameter that describes the storage of wfs
2297 : : !! kg_k(3,npw_k)=reduced planewave coordinates.
2298 : : !! kg1_k(3,npw1_k)=reduced planewave coordinates at k+q, with RF k points
2299 : : !! kpt(3)=reduced coordinates of k point
2300 : : !! kpq(3)=reduced coordinates of k+q point
2301 : : !! mkmem =number of k points treated by this node
2302 : : !! mk1mem =number of k points treated by this node (RF data)
2303 : : !! mpert =maximum number of ipert
2304 : : !! mpi_enreg=information about MPI parallelization
2305 : : !! mpw=maximum dimensioned size of npw or wfs at k
2306 : : !! mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
2307 : : !! nband_k=number of bands at this k point for that spin polarization
2308 : : !! npw_k=number of plane waves at this k point
2309 : : !! npw1_k=number of plane waves at this k+q point
2310 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
2311 : : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
2312 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
2313 : : !! rmet(3,3)=real space metric (bohr**2)
2314 : : !! ddks(3)<wfk_t>=struct info for for the three possible DDK files for ipert1
2315 : : !! wtk_k=weight assigned to the k point.
2316 : : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
2317 : : !! ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
2318 : : !!
2319 : : !! OUTPUT
2320 : : !! d2bbb_k(2,3,mband,mband*prtbbb)=band by band decomposition of the second
2321 : : !! order derivatives, for the present k point, and perturbation idir, ipert
2322 : : !! d2nl_k(2,3,mpert)=non-local contributions to
2323 : : !! non-stationary 2DTE, for the present k point, and perturbation idir, ipert
2324 : : !!
2325 : : !! TODO
2326 : : !! XG 20141103 The localization tensor cannot be defined in the metallic case. It should not be computed.
2327 : : !!
2328 : : !! SOURCE
2329 : :
2330 : 99323 : subroutine dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,&
2331 : 99323 : & icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpt,kpq,&
2332 : : & mband_mem_rbz,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,&
2333 : 99323 : & occ_k,psps,rmet,ddks,wtk_k,ylm,ylm1)
2334 : :
2335 : : !Arguments ------------------------------------
2336 : : !scalars
2337 : : integer,intent(in) :: icg,icg1,idir,ikpt,ipert,isppol,istwf_k
2338 : : integer,intent(in) :: mkmem,mk1mem,mpert,mpw,mpw1,nsppol
2339 : : integer,intent(in) :: mband_mem_rbz
2340 : : integer,intent(inout) :: nband_k,npw1_k,npw_k
2341 : : real(dp),intent(in) :: wtk_k
2342 : : type(MPI_type),intent(in) :: mpi_enreg
2343 : : type(dataset_type),intent(in) :: dtset
2344 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
2345 : : type(pseudopotential_type),intent(in) :: psps
2346 : : !arrays
2347 : : integer,intent(in) :: ddkfil(3),kg1_k(3,npw1_k)
2348 : : integer,intent(in) :: kg_k(3,npw_k)
2349 : : real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem_rbz*mkmem*nsppol)
2350 : : real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*nsppol)
2351 : : real(dp),intent(in) :: eig_k(dtset%mband*nsppol),kpt(3),kpq(3),occ_k(nband_k),rmet(3,3)
2352 : : real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
2353 : : real(dp),intent(in) :: ylm1(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)
2354 : : real(dp),intent(inout) :: eig1_k(2*nsppol*dtset%mband**2)
2355 : : real(dp),intent(out) :: d2bbb_k(2,3,dtset%mband,dtset%mband*dtset%prtbbb)
2356 : : real(dp),intent(inout) :: d2nl_k(2,3,mpert)
2357 : : type(wfk_t),intent(inout) :: ddks(3)
2358 : :
2359 : : !Local variables-------------------------------
2360 : : !scalars
2361 : : integer :: berryopt,dimffnl,dimffnl1,dimph3d
2362 : : integer :: iband,ider,idir1,ipert1,ipw,jband,nband_kocc,nkpg,nkpg1
2363 : : integer :: ierr, iband_me, jband_me
2364 : : integer :: npw_disk,nsp,optlocal,optnl,opt_gvnlx1,sij_opt,tim_getgh1c,usevnl
2365 : : integer :: nddk_needed, startband, endband
2366 : : logical :: ddk
2367 : : real(dp) :: aa,dot1i,dot1r,dot2i,dot2r,dot_ndiagi,dot_ndiagr,doti,dotr,lambda
2368 : : character(len=500) :: msg
2369 : 99323 : type(rf_hamiltonian_type) :: rf_hamkq
2370 : : !arrays
2371 : : integer :: ik_ddks(3)
2372 : 198646 : integer :: band_procs(nband_k)
2373 : 198646 : logical :: distrb_cycle(nband_k)
2374 : : real(dp) :: dum_grad_berry(1,1),dum_gvnlx1(1,1),dum_gs1(1,1),dum_ylmgr(1,3,1),tsec(2)
2375 : 99323 : real(dp),allocatable :: cg_k(:,:),cwave0(:,:),cwavef(:,:),cwavef_da(:,:)
2376 : 99323 : real(dp),allocatable :: cwaveddk(:,:,:)
2377 : 99323 : real(dp),allocatable :: cg_ddk(:,:,:) !2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*nsppol,3) ==
2378 : :
2379 : 99323 : real(dp),allocatable :: cwavef_db(:,:),dkinpw(:),eig2_k(:),ffnl1(:,:,:,:),ffnlk(:,:,:,:)
2380 : 99323 : real(dp),allocatable :: eig2_ddk(:,:)
2381 : 99323 : real(dp),allocatable :: gvnlx1(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),ph3d(:,:,:)
2382 : 99323 : type(pawcprj_type),allocatable :: dum_cwaveprj(:,:)
2383 : :
2384 : : ! *********************************************************************
2385 : :
2386 : : DBG_ENTER("COLL")
2387 : :
2388 : : ABI_NVTX_START_RANGE(NVTX_DFPT_NSTWF)
2389 : :
2390 : : !Not valid for PAW
2391 [ - + ]: 99323 : if (psps%usepaw==1) then
2392 : 0 : msg=' This routine cannot be used for PAW (use pawnst3 instead) !'
2393 : 0 : ABI_BUG(msg)
2394 : : end if
2395 : :
2396 : : !Keep track of total time spent in dfpt_nstwf
2397 : 99323 : call timab(112,1,tsec)
2398 : 99323 : tim_getgh1c=2
2399 : :
2400 : : !Miscellaneous inits
2401 [ - + - + ]: 99323 : ABI_MALLOC(dum_cwaveprj,(0,0))
2402 [ + + - + : 99323 : ddk=(ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11)
- + ]
2403 : :
2404 : : ! filter for bands on this cpu for cg cg1 etc.
2405 [ + + ]: 678003 : distrb_cycle = (mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol) /= mpi_enreg%me_kpt)
2406 : : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,nband_k,&
2407 : 99323 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
2408 : :
2409 : : !Additional allocations
2410 [ + + ]: 99323 : if (.not.ddk) then
2411 [ - + + - : 245859 : ABI_MALLOC(dkinpw,(npw_k))
+ - + - -
+ - + -
+ ]
2412 [ - + + - : 245859 : ABI_MALLOC(kinpw1,(npw1_k))
+ - + - -
+ - + -
+ ]
2413 [ + + + + ]: 26222013 : kinpw1=zero;dkinpw=zero
2414 : : else
2415 [ - + - + ]: 17370 : ABI_MALLOC(dkinpw,(0))
2416 [ - + - + ]: 17370 : ABI_MALLOC(kinpw1,(0))
2417 : : end if
2418 [ - + + - : 297969 : ABI_MALLOC(gvnlx1,(2,npw1_k*dtset%nspinor))
+ - + - -
+ - + ]
2419 [ + - - + : 397292 : ABI_MALLOC(eig2_k,(2*nsppol*dtset%mband**2))
+ - + - +
- - + -
+ ]
2420 [ - + + - : 297969 : ABI_MALLOC(cwave0,(2,npw_k*dtset%nspinor))
+ - + - -
+ - + ]
2421 [ - + + - : 198646 : ABI_MALLOC(cwavef,(2,npw1_k*dtset%nspinor))
- + - + ]
2422 : :
2423 : 99323 : nddk_needed = 0
2424 [ + + ]: 397292 : do idir1=1,3
2425 [ + + ]: 397292 : if (ddkfil(idir1)/=0) nddk_needed = nddk_needed+1
2426 : : end do
2427 [ + + ]: 99323 : if (nddk_needed > 0) then
2428 [ - + + - : 104720 : ABI_MALLOC(cwaveddk,(2,npw1_k*dtset%nspinor,3))
+ - + - +
- - + -
+ ]
2429 : : !TODO: for the moment avoid indirect indexing of the ddk directions in case not all are present. Here all are allocated and read in
2430 [ - + + - : 104720 : ABI_MALLOC(cg_ddk,(2,mpw1*dtset%nspinor*mband_mem_rbz,3))
+ - + - +
- - + -
+ ]
2431 [ + + + + : 137915825 : cg_ddk = zero ! not all may be initialized below if only certain ddk directions are provided
+ + ]
2432 : :
2433 [ + - - + : 104720 : ABI_MALLOC(eig2_ddk,(2*dtset%mband**2,3))
+ - + - +
- + - - +
- + ]
2434 [ + + + + ]: 4768418 : eig2_ddk = zero
2435 : : end if
2436 : :
2437 : : !Compute (k+G) vectors
2438 [ + + ]: 99323 : nkpg=0;if (.not.ddk) nkpg=3*gs_hamkq%nloalg(3)
2439 [ - + + - : 397292 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ - - - +
- - + - +
- + - + ]
2440 [ - + ]: 99323 : if (nkpg>0) then
2441 : 0 : call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
2442 : : end if
2443 : :
2444 : : !Compute (k+q+G) vectors
2445 [ + + ]: 99323 : nkpg1=0;if (.not.ddk) nkpg1=3*gs_hamkq%nloalg(3)
2446 [ - + + - : 397292 : ABI_MALLOC(kpg1_k,(npw1_k,nkpg1))
+ - - - +
- - + - +
- + - + ]
2447 [ - + ]: 99323 : if (nkpg1>0) then
2448 : 0 : call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k)
2449 : : end if
2450 : :
2451 : : !Compute nonlocal form factors ffnl at (k+G)
2452 [ + + ]: 99323 : dimffnl=0;if (.not.ddk) dimffnl=1
2453 [ - + + - : 595938 : ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
+ + + - -
+ + - - +
+ - + - +
+ - + - +
- + ]
2454 [ + + ]: 99323 : if (.not.ddk) then
2455 : 81953 : ider=0
2456 : : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,&
2457 : : & gs_hamkq%gprimd,ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,&
2458 : : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,&
2459 : 81953 : & psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,dum_ylmgr)
2460 : : end if
2461 : :
2462 : : !Compute nonlocal form factors ffnl1 at (k+q+G)
2463 : 99323 : dimffnl1=0;if (.not.ddk) dimffnl1=1
2464 [ - + + - : 595938 : ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
+ + + - -
+ + - - +
+ - + - +
+ - + - +
- + ]
2465 [ + + ]: 99323 : if (.not.ddk) then
2466 : 81953 : ider=0
2467 : : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,&
2468 : : & gs_hamkq%gprimd,ider,ider,psps%indlmn,kg1_k,kpg1_k,kpq,&
2469 : : & psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,&
2470 : 81953 : & psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1,dum_ylmgr)
2471 : : end if
2472 : :
2473 : : !Load k-dependent part in the Hamiltonian datastructure
2474 : : call gs_hamkq%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,&
2475 : 99323 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk)
2476 : :
2477 : : !Load k+q-dependent part in the Hamiltonian datastructure
2478 : : ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
2479 [ + + ]: 99323 : dimph3d=0;if (.not.ddk) dimph3d=gs_hamkq%matblk
2480 [ - + + - : 397292 : ABI_MALLOC(ph3d,(2,npw1_k,dimph3d))
+ + + - +
- + + - +
- + - + ]
2481 : : call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
2482 : : & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,&
2483 : 99323 : & ph3d_kp=ph3d,compute_ph3d=(.not.ddk))
2484 : :
2485 : : !Load k-dependent part in the 1st-order Hamiltonian datastructure
2486 : 99323 : call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw)
2487 : :
2488 : : !Take care of the npw and kg records
2489 : : !NOTE : one should be able to modify the rwwf routine to take care
2490 : : !of the band parallelism, which is not the case yet ...
2491 : 99323 : ik_ddks = 0
2492 [ + + ]: 397292 : do idir1=1,3
2493 [ + + ]: 397292 : if (ddkfil(idir1)/=0)then
2494 : : ! Read npw record
2495 : 75546 : nsp=dtset%nspinor
2496 : 75546 : ik_ddks(idir1) = ddks(idir1)%findk(kpt)
2497 [ - + ]: 75546 : ABI_CHECK(ik_ddks(idir1) /= -1, "Cannot find kpt")
2498 : 75546 : npw_disk = ddks(idir1)%hdr%npwarr(ik_ddks(idir1))
2499 [ - + ]: 75546 : if (npw_k /= npw_disk) then
2500 : : write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')&
2501 : 0 : & 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,&
2502 : 0 : & 'the number of plane waves in the ddk file is equal to', npw_disk,ch10,&
2503 : 0 : & 'while it should be ',npw_k
2504 : 0 : ABI_BUG(msg)
2505 : : end if
2506 : :
2507 : : ! NB: this will fail if the bands are not contiguous.
2508 : 75546 : startband = nband_k
2509 : 75546 : endband = 1
2510 [ + + ]: 406077 : do iband=1,nband_k
2511 [ + - ]: 406077 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then
2512 : 330531 : if (iband < startband) startband = iband
2513 : 330531 : if (iband > endband) endband = iband
2514 : : end if
2515 : : end do
2516 : : ! NB: eig_k is band distributed in call to read_band_block, though array has full size,
2517 : : ! only certain columns for my iband are filled, then used below
2518 : : call ddks(idir1)%read_band_block((/startband,endband/),ik_ddks(idir1),isppol,xmpio_collective, &
2519 [ + + ]: 226638 : & cg_k=cg_ddk(:,:,idir1), eig_k=eig2_ddk(:,idir1))
2520 : : end if ! ddk file is already present
2521 : : end do ! idir1
2522 : :
2523 [ + + ]: 99323 : if (ipert==dtset%natom+1) then
2524 : 17370 : nband_kocc = 0
2525 [ + + ]: 129761 : do iband = 1,nband_k
2526 [ + + ]: 112391 : if (abs(occ_k(iband)) > tol8) nband_kocc = nband_kocc + 1
2527 : 129761 : nband_kocc = max (nband_kocc, 1)
2528 : : end do
2529 : : end if
2530 : :
2531 [ + + ]: 99323 : if(dtset%prtbbb==1)then
2532 [ - + + - : 456 : ABI_MALLOC(cwavef_da,(2,npw1_k*dtset%nspinor))
+ - + - -
+ - + ]
2533 [ - + + - : 304 : ABI_MALLOC(cwavef_db,(2,npw1_k*dtset%nspinor))
- + - + ]
2534 [ - + + - : 456 : ABI_MALLOC(cg_k,(2,npw_k*dtset%nspinor*mband_mem_rbz))
+ - + - -
+ - + ]
2535 : : if ((ipert == dtset%natom + 1).or.(ipert <= dtset%natom).or. &
2536 [ + + + + : 152 : & (ipert == dtset%natom + 2).or.(ipert == dtset%natom + 5)) then
- + - - ]
2537 [ + + + + ]: 533912 : cg_k(:,:) = cg(:,1+icg:icg+mband_mem_rbz*npw_k*dtset%nspinor)
2538 : : end if
2539 [ + + + + : 98648 : d2bbb_k(:,:,:,:) = zero
+ + + + ]
2540 : : end if
2541 : :
2542 : : !Loop over ALL bands
2543 : 99323 : iband_me = 0
2544 [ + + ]: 678003 : do iband=1,nband_k
2545 : : ABI_NVTX_START_RANGE(NVTX_DFPT_NSTWF_BAND)
2546 : :
2547 : : ! if band is mine, retrieve it and then broadcast it
2548 [ + + ]: 578680 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then
2549 : 578464 : iband_me = iband_me + 1
2550 : :
2551 : : ! Read ground-state wavefunction for iband
2552 [ + + + + ]: 578464 : if (dtset%prtbbb==0 .or. ipert==dtset%natom+2) then
2553 [ + + + + ]: 365415724 : cwave0(:,:)=cg(:,1+(iband_me-1)*npw_k*dtset%nspinor+icg:iband_me*npw_k*dtset%nspinor+icg)
2554 : : else ! prtbbb==1 and ipert<=natom , already in cg_k
2555 [ + + + + ]: 422400 : cwave0(:,:)=cg_k(:,1+(iband_me-1)*npw_k*dtset%nspinor:iband_me*npw_k*dtset%nspinor)
2556 : : end if
2557 : :
2558 : : ! Get first-order wavefunctions for iband
2559 [ + + + + ]: 365760304 : cwavef(:,:)=cg1(:,1+(iband_me-1)*npw1_k*dtset%nspinor+icg1:iband_me*npw1_k*dtset%nspinor+icg1)
2560 : : ! Get ddk wavefunctions for iband
2561 [ + + ]: 578464 : if(nddk_needed > 0) then
2562 [ + + + + : 132778161 : cwaveddk(:,:,:)=cg_ddk(:,1+(iband_me-1)*npw1_k*dtset%nspinor:iband_me*npw1_k*dtset%nspinor,:)
+ + ]
2563 : : end if
2564 : : end if
2565 : 578680 : call xmpi_bcast(cwave0, band_procs(iband), mpi_enreg%comm_band, ierr)
2566 : 578680 : call xmpi_bcast(cwavef, band_procs(iband), mpi_enreg%comm_band, ierr)
2567 [ + + ]: 578680 : if(nddk_needed > 0) then
2568 : 114975 : call xmpi_bcast(cwaveddk, band_procs(iband), mpi_enreg%comm_band, ierr)
2569 : : end if
2570 : :
2571 : : ! In case non ddk perturbation
2572 [ + + ]: 578680 : if (ipert /= dtset%natom + 1) then
2573 : :
2574 [ + + ]: 4693048 : do ipert1=1,mpert
2575 : :
2576 [ + + + + ]: 4693048 : if( ipert1<=dtset%natom .or. ipert1==dtset%natom+2 )then
2577 : :
2578 : : ! Initialize data for NL 1st-order hamiltonian
2579 : 1328225 : call rf_hamkq%init(1,gs_hamkq,ipert1)
2580 : :
2581 : : if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)) &
2582 [ + + + + : 1328225 : & .and.(ipert1 == dtset%natom+2).and. dtset%prtbbb==1) then
+ + + + ]
2583 : : call gaugetransfo(cg_k,cwavef,cwavef_db,mpi_enreg%comm_band,distrb_cycle,eig_k,eig1_k,&
2584 : : & dtset%ggtrcut,iband,nband_k,dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,&
2585 : 448 : & nsppol,mpi_enreg%nproc_band,occ_k)
2586 [ + + + + ]: 197248 : cwavef(:,:) = cwavef_db(:,:)
2587 : : end if
2588 : :
2589 : : ! Define the direction along which to move the atom :
2590 : : ! the polarisation (ipert1,idir1) is referred as j1.
2591 [ + + ]: 5312900 : do idir1=1,3
2592 [ + + + - : 5312900 : if (ipert1<=dtset%natom.or.(ipert1==dtset%natom+2.and.ddkfil(idir1)/=0)) then
+ + ]
2593 : :
2594 : : ! Get |Vnon-locj^(1)|u0> :
2595 : : ! First-order non-local, applied to zero-order wavefunction
2596 : : ! This routine gives MINUS the non-local contribution
2597 : :
2598 : : ! ==== Atomic displ. perturbation
2599 [ + + ]: 2915187 : if( ipert1<=dtset%natom )then
2600 : 2585808 : lambda=eig_k((isppol-1)*nband_k+iband)
2601 : 2585808 : berryopt=1;optlocal=0;optnl=1;usevnl=0;opt_gvnlx1=0;sij_opt=0
2602 : : call getgh1c(berryopt,cwave0,dum_cwaveprj,gvnlx1,dum_grad_berry,&
2603 : : & dum_gs1,gs_hamkq,dum_gvnlx1,idir1,ipert1,(/lambda/),mpi_enreg,1,optlocal,&
2604 [ + + ]: 5171616 : & optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
2605 : :
2606 : : ! ==== Electric field perturbation
2607 [ - + ]: 329379 : else if( ipert1==dtset%natom+2 )then
2608 : : ! TODO: Several tests fail here ifdef HAVE_MPI_IO_DEFAULT
2609 : : ! The problem is somehow related to the use of MPI-IO file views!.
2610 : : !TODO MJV: Check if it works with HAVE_MPI_IO_DEFAULT now.
2611 : :
2612 [ + - + - : 128343861 : gvnlx1 = cwaveddk(:,:,idir1)
- + - - +
+ + + ]
2613 [ + + ]: 4699101 : eig2_k(1+(iband-1)*2*nband_k:iband*2*nband_k) = eig2_ddk(1+(iband-1)*2*nband_k:iband*2*nband_k,idir1)
2614 : :
2615 : : !write(777,*)"eig2_k, gvnlx1 for band: ",iband,", ikpt: ",ikpt
2616 : : !do ii=1,2*nband_k
2617 : : ! write(777,*)eig2_k(ii+(iband-1))
2618 : : !end do
2619 : : !write(777,*)gvnlx1
2620 : :
2621 : : ! In case of band-by-band,
2622 : : ! construct the first-order wavefunctions in the diagonal gauge
2623 [ + + + + : 329379 : if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)).and.(dtset%prtbbb==1)) then
+ + ]
2624 : : call gaugetransfo(cg_k,gvnlx1,cwavef_da,mpi_enreg%comm_band,distrb_cycle,eig_k,eig2_k,&
2625 : : & dtset%ggtrcut,iband,nband_k,dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,&
2626 : 1344 : & nsppol,mpi_enreg%nproc_band,occ_k)
2627 [ + + + + ]: 591744 : gvnlx1(:,:) = cwavef_da(:,:)
2628 : : end if
2629 : : ! Multiplication by -i
2630 [ + + ]: 42891080 : do ipw=1,npw1_k*dtset%nspinor
2631 : 42561701 : aa=gvnlx1(1,ipw)
2632 : 42561701 : gvnlx1(1,ipw)=gvnlx1(2,ipw)
2633 : 42891080 : gvnlx1(2,ipw)=-aa
2634 : : end do
2635 : :
2636 : : ! MRoyo 030925 :
2637 : : ! 3) Case ipert1=natom+2 and ipert=natom+5 or ipert=natom+12:2*natom+11
2638 : : ! the computation of mixed derivatives wrt to Zeeman and electric fields needs $i \frac{d}{dk}.
2639 [ + + + + : 7450899 : if (ipert==dtset%natom+5.or.(ipert>=dtset%natom+12.and.ipert<=2*dtset%natom+11)) gvnlx1(:,:) = -gvnlx1(:,:)
+ - + + +
+ ]
2640 : :
2641 : : end if
2642 : :
2643 : : ! at this stage if iband is not mine I can cycle
2644 [ + + ]: 2915187 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) then
2645 : : cycle
2646 : : end if
2647 : :
2648 : : ! MVeithen 021212 :
2649 : : ! 1) Case ipert1 = natom + 2 and ipert = natom + 2:
2650 : : ! the second derivative of the energy with respect to an electric
2651 : : ! field is computed from Eq. (38) of X. Gonze, PRB 55 ,10355 (1997) [[cite:Gonze1997a]].
2652 : : ! The evaluation of this formula needs the operator $i \frac{d}{dk}.
2653 : : ! 2) Case ipert1 = natom + 2 and ipert < natom:
2654 : : ! the computation of the Born effective charge tensor uses
2655 : : ! the operator $-i \frac{d}{dk}.
2656 [ + + + + : 184161492 : if (ipert==dtset%natom+2) gvnlx1(:,:) = -gvnlx1(:,:)
+ + ]
2657 : :
2658 : :
2659 : : ! <G|Vnl1|Cnk> is contained in gvnlx1
2660 : : ! construct the matrix element (<uj2|vj1|u0>)complex conjug and add it to the 2nd-order matrix
2661 : 2909355 : call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,cwavef,gvnlx1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2662 : 2909355 : d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*two*dotr
2663 : 2909355 : d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*two*doti
2664 : :
2665 : : ! Band by band decomposition of the Born effective charges
2666 : : ! calculated from a phonon perturbation
2667 [ + + ]: 2909355 : if(dtset%prtbbb==1) then ! .and. mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt)then
2668 : 4032 : d2bbb_k(1,idir1,iband,iband) = wtk_k*occ_k(iband)*two*dotr
2669 : 4032 : d2bbb_k(2,idir1,iband,iband) = -one*wtk_k*occ_k(iband)*two*doti
2670 : : end if
2671 : :
2672 : : end if
2673 : : end do ! idir
2674 : :
2675 : 1328225 : call rf_hamkq%free()
2676 : : end if ! ipert1<=dtset%natom .or. ipert1==dtset%natom+2
2677 : : end do ! ipert1
2678 : : end if ! ipert /= natom +1
2679 : :
2680 : : ! Compute the localization tensor
2681 : :
2682 [ + + ]: 1256683 : if (ipert==dtset%natom+1) then
2683 : :
2684 : 112391 : ipert1=dtset%natom+1
2685 [ + + ]: 112391 : if(dtset%prtbbb==1)then
2686 : : call gaugetransfo(cg_k,cwavef,cwavef_db,mpi_enreg%comm_band,distrb_cycle,eig_k,eig1_k,&
2687 : : & dtset%ggtrcut,iband,nband_k,dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,&
2688 : 768 : & nsppol,mpi_enreg%nproc_band,occ_k)
2689 [ + + + + ]: 337728 : cwavef(:,:) = cwavef_db(:,:)
2690 : : end if
2691 : :
2692 [ + + ]: 449564 : do idir1 = 1,3
2693 [ + + ]: 143637471 : eig2_k(:) = zero
2694 [ + + + + ]: 116095911 : gvnlx1(:,:) = zero
2695 [ + + ]: 337173 : if (idir == idir1) then
2696 [ + + + + ]: 38698637 : gvnlx1(:,:) = cwavef(:,:)
2697 [ + + ]: 47879157 : eig2_k(:) = eig1_k(:)
2698 : : else
2699 [ + + ]: 224782 : if (ddkfil(idir1) /= 0) then
2700 [ + - + - : 356544 : gvnlx1 = cwaveddk(:,:,idir1)
- + - - +
+ + + ]
2701 [ + + ]: 16512 : eig2_k(1+(iband-1)*2*nband_k:iband*2*nband_k) = eig2_ddk(1+(iband-1)*2*nband_k:iband*2*nband_k,idir1)
2702 : :
2703 : : !write(778,*)"eig2_k, gvnlx1 for band: ",iband,", ikpt: ",ikpt
2704 : : !do ii=1,2*nband_k
2705 : : ! write(778,*)eig2_k(ii+(iband-1))
2706 : : !end do
2707 : : !write(778,*)gvnlx1
2708 : :
2709 [ + + ]: 1152 : if(dtset%prtbbb==1)then
2710 : : call gaugetransfo(cg_k,gvnlx1,cwavef_da,mpi_enreg%comm_band,distrb_cycle,eig_k,eig2_k,&
2711 : : & dtset%ggtrcut,iband,nband_k,dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,&
2712 : 768 : & nsppol,mpi_enreg%nproc_band,occ_k)
2713 : :
2714 [ + + + + ]: 337728 : gvnlx1(:,:) = cwavef_da(:,:)
2715 : : end if
2716 : :
2717 : : end if !ddkfil(idir1)
2718 : : end if !idir == idir1
2719 : :
2720 : : ! <G|du/dqa> is contained in gvnlx1 and <G|du/dqb> in cwavef
2721 : : ! construct the matrix elements <du/dqa|du/dqb> -> dot
2722 : : ! <u|du/dqa> -> dot1
2723 : : ! <du/dqb|u> -> dot2
2724 : : ! and add them to the 2nd-order matrix
2725 : :
2726 : 337173 : call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,gvnlx1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2727 : 337173 : d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two)
2728 : 337173 : d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two)
2729 : :
2730 : :
2731 : : ! XG 020216 : Marek, could you check the next forty lines
2732 : : ! In the parallel gauge, dot1 and dot2 vanishes
2733 [ + + ]: 449564 : if(dtset%prtbbb==1)then
2734 [ + - ]: 2304 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then
2735 : 2304 : d2bbb_k(1,idir1,iband,iband)=d2bbb_k(1,idir1,iband,iband)+dotr
2736 : 2304 : d2bbb_k(2,idir1,iband,iband)=d2bbb_k(2,idir1,iband,iband)+doti
2737 : : end if
2738 : 2304 : dot_ndiagr=zero ; dot_ndiagi=zero
2739 : 2304 : jband_me = 0
2740 [ + + ]: 20736 : do jband = 1,nband_k !compute dot1 and dot2
2741 [ + - ]: 18432 : if (mpi_enreg%proc_distrb(ikpt,jband,isppol) /= mpi_enreg%me_kpt) then
2742 : : cycle
2743 : : end if
2744 : 18432 : jband_me = jband_me + 1
2745 : :
2746 [ + - ]: 20736 : if (abs(occ_k(jband)) > tol8) then
2747 : 18432 : dot1r=zero ; dot1i=zero
2748 : 18432 : dot2r=zero ; dot2i=zero
2749 [ + + + + ]: 8105472 : cwave0(:,:)=cg_k(:,1+(jband_me-1)*npw_k*dtset%nspinor:jband_me*npw_k*dtset%nspinor)
2750 : :
2751 : 18432 : call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*dtset%nspinor,2,cwave0,gvnlx1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2752 : 18432 : call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*dtset%nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2753 : :
2754 : 18432 : dot_ndiagr = dot_ndiagr + dot1r*dot2r - dot1i*dot2i
2755 : 18432 : dot_ndiagi = dot_ndiagi + dot1r*dot2i + dot1i*dot2r
2756 : : ! this should fill all of the iband but only the local cpu jband indices
2757 : : d2bbb_k(1,idir1,iband,jband) = d2bbb_k(1,idir1,iband,jband) - &
2758 : 18432 : & (dot1r*dot2r - dot1i*dot2i)
2759 : : d2bbb_k(2,idir1,iband,jband) = d2bbb_k(2,idir1,iband,jband) - &
2760 : 18432 : & (dot1r*dot2i + dot1i*dot2r)
2761 : : end if ! occ_k
2762 : : end do !jband
2763 : :
2764 [ + + + + ]: 57600 : d2bbb_k(:,idir1,iband,:)=d2bbb_k(:,idir1,iband,:)*wtk_k*occ_k(iband)*half
2765 : : d2nl_k(1,idir1,ipert1)= &
2766 : 2304 : & d2nl_k(1,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagr/(nband_kocc*two)
2767 : : d2nl_k(2,idir1,ipert1)=&
2768 : 2304 : & d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagi/(nband_kocc*two)
2769 : : end if ! prtbbb==1
2770 : :
2771 : : end do ! idir1
2772 : : end if ! Compute localization tensor, ipert=natom+1
2773 : :
2774 : : ABI_NVTX_END_RANGE()
2775 : : end do ! End loop over iband
2776 : :
2777 : :
2778 : : ! if(dtset%prtbbb==1)then
2779 : : ! ! complete over jband index
2780 : : ! call xmpi_sum(d2bbb_k, mpi_enreg%comm_band, ierr)
2781 : : ! end if
2782 : :
2783 : : !Final deallocations
2784 [ - + ]: 99323 : ABI_FREE(cwave0)
2785 [ - + ]: 99323 : ABI_FREE(cwavef)
2786 [ - + ]: 99323 : ABI_FREE(eig2_k)
2787 [ - + ]: 99323 : ABI_FREE(gvnlx1)
2788 [ - + ]: 99323 : ABI_FREE(ffnlk)
2789 [ - + ]: 99323 : ABI_FREE(ffnl1)
2790 [ - + ]: 99323 : ABI_FREE(dkinpw)
2791 [ - + ]: 99323 : ABI_FREE(kinpw1)
2792 [ - + ]: 99323 : ABI_FREE(kpg_k)
2793 [ - + ]: 99323 : ABI_FREE(kpg1_k)
2794 [ - + ]: 99323 : ABI_FREE(ph3d)
2795 [ + - - + : 99323 : ABI_FREE(dum_cwaveprj)
- - - - -
+ ]
2796 [ + + ]: 99323 : if(dtset%prtbbb==1) then
2797 [ - + ]: 152 : ABI_FREE(cg_k)
2798 [ - + ]: 152 : ABI_FREE(cwavef_da)
2799 [ - + ]: 152 : ABI_FREE(cwavef_db)
2800 : : end if
2801 [ + + ]: 99323 : if (nddk_needed > 0) then
2802 [ - + ]: 26180 : ABI_FREE(cwaveddk)
2803 [ - + ]: 26180 : ABI_FREE(cg_ddk)
2804 [ - + ]: 26180 : ABI_FREE(eig2_ddk)
2805 : : end if
2806 : :
2807 : 99323 : call timab(112,2,tsec)
2808 : :
2809 : : ABI_NVTX_END_RANGE()
2810 : :
2811 : : DBG_EXIT("COLL")
2812 : :
2813 [ - + - + : 99323 : end subroutine dfpt_nstwf
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - + -
+ - + - +
- + - + ]
2814 : : !!***
2815 : :
2816 : : !!****f* ABINIT/gaugetransfo
2817 : : !! NAME
2818 : : !! gaugetransfo
2819 : : !!
2820 : : !! FUNCTION
2821 : : !! This routine allows the passage from the parallel-transport gauge
2822 : : !! to the diagonal gauge for the first-order wavefunctions
2823 : : !!
2824 : : !! INPUTS
2825 : : !! cg_k(2,mpw*nspinor*mband_mem*nsppol)=planewave coefficients of wavefunctions
2826 : : !! for a particular k point.
2827 : : !! cwavef(2,npw1_k*nspinor)=first order wavefunction for a particular k point
2828 : : !! in the parallel gauge
2829 : : !! comm=mpi communicator for bands
2830 : : !! distrb_cycle=array of logical flags to skip certain bands in parallelization scheme
2831 : : !! eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
2832 : : !! eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
2833 : : !! ggtrcut=cutoff to reject unstable band pairs
2834 : : !! iband=band index of the 1WF for which the transformation has to be applied
2835 : : !! mband=maximum number of bands
2836 : : !! mband_mem_rbz=maximum number of bands on this cpu
2837 : : !! nband_k=number of bands for this k point
2838 : : !! npw_k=maximum dimensioned size of npw or wfs at k
2839 : : !! npw1_k=number of plane waves at this k+q point
2840 : : !! nspinor=number of spinorial components of the wavefunctions
2841 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
2842 : : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k
2843 : : !!
2844 : : !! OUTPUT
2845 : : !! cwavef_d(2,npw1_k*nspinor)=first order wavefunction for a particular k point
2846 : : !! in the diagonal gauge
2847 : : !!
2848 : : !! SOURCE
2849 : :
2850 : 3328 : subroutine gaugetransfo(cg_k,cwavef,cwavef_d,comm,distrb_cycle,eig_k,eig1_k,ggtrcut,iband,nband_k, &
2851 : 3328 : & mband,mband_mem_rbz,npw_k,npw1_k,nspinor,nsppol,nproc_band,occ_k)
2852 : :
2853 : : !Arguments ------------------------------------
2854 : : !scalars
2855 : : integer,intent(in) :: iband,mband,mband_mem_rbz,nband_k,npw1_k,npw_k,nspinor,nsppol
2856 : : integer,intent(in) :: comm, nproc_band
2857 : : real(dp),intent(in) :: ggtrcut
2858 : : !arrays
2859 : : logical, intent(in) :: distrb_cycle(nband_k)
2860 : : real(dp),intent(in) :: cg_k(2,npw_k*nspinor*mband_mem_rbz),cwavef(2,npw1_k*nspinor)
2861 : : real(dp),intent(in) :: eig1_k(2*nsppol*mband**2),eig_k(mband*nsppol)
2862 : : real(dp),intent(in) :: occ_k(nband_k)
2863 : : real(dp),intent(out) :: cwavef_d(2,npw1_k*nspinor)
2864 : :
2865 : : !Local variables-------------------------------
2866 : : !tolerance for non degenerated levels
2867 : : !scalars
2868 : : integer :: ierr, jband,jband_me
2869 : : !arrays
2870 : 3328 : real(dp) :: cwave0(2,npw1_k*nspinor),eig1(2)
2871 : :
2872 : : ! *********************************************************************
2873 : :
2874 [ + + + + ]: 1464448 : cwavef_d(:,:) = cwavef(:,:)
2875 : :
2876 : : jband_me = 0
2877 [ + + ]: 29952 : do jband = 1,nband_k !loop over bands
2878 [ + - ]: 26624 : if (distrb_cycle(jband)) cycle
2879 : 26624 : jband_me = jband_me + 1
2880 : :
2881 : : ! the default value of ggtrcut is 0.001, was formerly a parameter in this routine
2882 : : ! called etol=0.001_dp
2883 [ + + + - ]: 29952 : if ((abs(eig_k(iband)-eig_k(jband)) > ggtrcut).and.(abs(occ_k(jband)) > tol8 )) then
2884 : :
2885 [ + + + + ]: 9948064 : cwave0(:,:) = cg_k(:,1+(jband_me-1)*npw_k*nspinor:jband_me*npw_k*nspinor)
2886 : :
2887 : 22624 : eig1(1) = eig1_k(2*jband-1+(iband-1)*2*nband_k)
2888 : 22624 : eig1(2) = eig1_k(2*jband +(iband-1)*2*nband_k)
2889 : :
2890 : : cwavef_d(1,:)=cwavef_d(1,:) &
2891 [ + + ]: 3331104 : & - (eig1(1)*cwave0(1,:)-eig1(2)*cwave0(2,:))/(eig_k(jband)-eig_k(iband))
2892 : : cwavef_d(2,:)=cwavef_d(2,:) &
2893 [ + + ]: 3331104 : & - (eig1(1)*cwave0(2,:)+eig1(2)*cwave0(1,:))/(eig_k(jband)-eig_k(iband))
2894 : :
2895 : : end if
2896 : :
2897 : : end do !loop over bands
2898 : 3328 : call xmpi_sum(cwavef_d, comm, ierr)
2899 : : ! here we have summed the cwavef N times (N-1 too many), but the correction is completed over bands
2900 [ + + + + ]: 1464448 : cwavef_d = cwavef_d - dble(nproc_band-1)*cwavef
2901 : :
2902 : 3328 : end subroutine gaugetransfo
2903 : : !!***
2904 : :
2905 : :
2906 : : end module m_dfpt_nstwf
2907 : : !!***
|
