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