Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfpt_nstwf
2 : : !! NAME
3 : : !! m_dfpt_nstwf
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group ()
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : ! nvtx related macro definition
22 : : #include "nvtx_macros.h"
23 : :
24 : : module m_dfpt_nstwf
25 : :
26 : : use defs_basis
27 : : use m_xmpi
28 : : use m_errors
29 : : use m_abicore
30 : : use m_wfk
31 : : use m_hamiltonian
32 : : use m_cgtools
33 : : use m_nctk
34 : : use m_dtset
35 : : use m_dtfil
36 : : use m_gputk
37 : : use m_abi_linalg
38 : :
39 : : use defs_datatypes, only : pseudopotential_type
40 : : use defs_abitypes, only : MPI_type
41 : : use m_time, only : timab
42 : : use m_io_tools, only : file_exists
43 : : use m_fourier_interpol, only : transgrid
44 : : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, gemm_nonlop_use_gemm
45 : : use m_geometry, only : stresssym
46 : : use m_dynmat, only : dfpt_sygra
47 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle, proc_distrb_band, proc_distrb_nband
48 : : use m_hdr, only : hdr_skip
49 : : use m_occ, only : occeig
50 : : use m_pawang, only : pawang_type
51 : : use m_pawrad, only : pawrad_type
52 : : use m_pawtab, only : pawtab_type
53 : : use m_paw_an, only : paw_an_type, paw_an_reset_flags
54 : : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
55 : : use m_pawfgrtab,only : pawfgrtab_type
56 : : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, pawrhoij_nullify, &
57 : : pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim
58 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_copy, pawcprj_output
59 : : use m_pawdij, only : pawdijfr
60 : : use m_pawfgr, only : pawfgr_type
61 : : use m_paw_mkrho,only : pawmkrho
62 : : use m_paw_nhat, only : pawnhatfr
63 : : use m_paw_dfpt, only : pawdfptenergy
64 : : use m_kg, only : mkkin, kpgstr, mkkpg
65 : : use m_fft, only : fftpac
66 : : use m_spacepar, only : hartrestr, symrhg
67 : : use m_initylmg, only : initylmg
68 : : use m_mkffnl, only : mkffnl
69 : : use m_getgh1c, only : getgh1c, getdc1
70 : : use m_dfpt_mkrho, only : dfpt_accrho
71 : : use m_atm2fft, only : dfpt_atm2fft
72 : : use m_mkcore, only : dfpt_mkcore
73 : : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll
74 : : use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
75 : : use m_mklocl, only : dfpt_vlocal, vlocalstr
76 : : use m_cgprj, only : getcprj
77 : :
78 : : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
79 : :
80 : : #if defined(HAVE_GPU_MARKERS)
81 : : use m_nvtx_data
82 : : #endif
83 : :
84 : : implicit none
85 : :
86 : : private
87 : : !!***
88 : :
89 : : public :: dfpt_nstpaw
90 : : public :: dfpt_nstwf
91 : : public :: gaugetransfo
92 : : !!***
93 : :
94 : : contains
95 : : !!***
96 : :
97 : : !!****f* ABINIT/dfpt_nstpaw
98 : : !! NAME
99 : : !! dfpt_nstpaw
100 : : !!
101 : : !! FUNCTION
102 : : !! Initially designed for PAW approach, but works also for NCPP.
103 : : !! This routine compute the non-stationary expression for the
104 : : !! second derivative of the total energy, for a whole row of
105 : : !! mixed derivatives (including diagonal terms contributing
106 : : !! to non-stationary 2nd-order total energy).
107 : : !! Compared with NC-pseudopotentials, PAW contributions include:
108 : : !! - changes of the overlap between 0-order wave-functions,
109 : : !! - on-site contributions.
110 : : !!
111 : : !! COPYRIGHT
112 : : !! Copyright (C) 2010-2026 ABINIT group (MT, AM)
113 : : !! This file is distributed under the terms of the
114 : : !! GNU General Public License, see ~abinit/COPYING
115 : : !! or http://www.gnu.org/copyleft/gpl.txt .
116 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
117 : : !!
118 : : !! INPUTS
119 : : !! cg (2,mpw *nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions at k
120 : : !! cgq(2,mpw1*nspinor*mband_mem*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
121 : : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
122 : : !! cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
123 : : !! cprj(natom,nspinor*mband_mem*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors
124 : : !! cprjq(natom,nspinor*mband_mem*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors
125 : : !! docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
126 : : !! doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
127 : : !! dtfil <type(datafiles_type)>=variables related to files
128 : : !! dtset <type(dataset_type)>=all input variables for this dataset
129 : : !! eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
130 : : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
131 : : !! eigen1(2*mband*mband*nkpt_rbz*nsppol)=1st-order eigenvalues at k,q (hartree)
132 : : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
133 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
134 : : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
135 : : !! idir=direction of the perturbation
136 : : !! indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points
137 : : !! indsy1(4,nsym1,natom)=indirect indexing array for atom labels
138 : : !! ipert=type of the perturbation
139 : : !! irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for RF symmetries
140 : : !! istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
141 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
142 : : !! kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
143 : : !! kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ
144 : : !! kxc(nfftf,nkxc)=exchange and correlation kernel
145 : : !! mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90)
146 : : !! mkmem =number of k points treated by this node.
147 : : !! mkqmem =number of k+q points treated by this node (GS data).
148 : : !! mk1mem =number of k points treated by this node (RF data)
149 : : !! mpert =maximum number of ipert
150 : : !! mpi_enreg=information about MPI parallelization
151 : : !! mpw=maximum dimensioned size of npw or wfs at k
152 : : !! mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
153 : : !! nattyp(ntypat)= # atoms of each type.
154 : : !! nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
155 : : !! mband_mem=number of bands per processor
156 : : !! ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
157 : : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid
158 : : !! ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid
159 : : !! nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
160 : : !! nhat1(cplex*nfftf,nspden*usepaw)=1st-order compensation charge density (PAW)
161 : : !! nkpt_rbz=number of k points in the reduced BZ for this perturbation
162 : : !! nkxc=second dimension of the kxc array
163 : : !! npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
164 : : !! npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
165 : : !! nspden=number of spin-density components
166 : : !! nspinor=number of spinorial components of the wavefunctions
167 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
168 : : !! nsym1=number of symmetry elements in space group consistent with i perturbation
169 : : !! n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, cplex*nfftf
170 : : !! occkq(mband*nkpt_rbz*nsppol)=occupation number for each band at each k+q point of the reduced BZ
171 : : !! occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced BZ
172 : : !! paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS
173 : : !! paw_an1(natom) <type(paw_an_type)>=1st-order paw arrays given on angular mesh for the perturbation (j1)
174 : : !! paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
175 : : !! paw_ij1(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels
176 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
177 : : !! pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
178 : : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
179 : : !! pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
180 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
181 : : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
182 : : !! pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
183 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
184 : : !! phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic transl. phases, for RF symmetries
185 : : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
186 : : !! ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information for the "fine" grid
187 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
188 : : !! rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
189 : : !! rhor1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3.
190 : : !! rmet(3,3)=real space metric (bohr**2)
191 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
192 : : !! symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations
193 : : !! symrc1(3,3,nsym1)=symmetry operations in reciprocal space
194 : : !! symrl1(3,3,nsym1)=symmetry operations in real space in terms
195 : : !! ucvol=unit cell volume in bohr**3.
196 : : !! usecprj= 1 if cprj, cprjq arrays are stored in memory
197 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
198 : : !! usevxctau=1 if if XC functional depends on kinetic energy density
199 : : !! usexcnhat= -PAW only- flag controlling use of compensation density in Vxc
200 : : !! useylmgr1= 1 if ylmgr1 array is allocated
201 : : !! vectornd(with_vectornd*nfftf,3)=nuclear dipole moment vector potential
202 : : !! vhartr1(cplex*nfft)=1-order Hartree potential
203 : : !! vpsp1(cplex*nfftf)=first-order derivative of the ionic potential
204 : : !! vtrial(nfftf,nspden)=GS potential (Hartree).
205 : : !! vtrial1(cplex*nfftf,nspden)= RF 1st-order potential (Hartree).
206 : : !! vxc(nfftf,nspden)=XC GS potential
207 : : !! vxctau(nfftf,nspden,4*usevxctau)=derivative of e_xc with respect to kinetic energy density, for mGGA
208 : : !! with_vectornd = 1 if vectornd allocated
209 : : !! wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
210 : : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
211 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
212 : : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
213 : : !! ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
214 : : !! ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q
215 : : !!
216 : : !! OUTPUT
217 : : !! blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
218 : : !! d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
219 : : !! d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
220 : : !! d2ovl(2,3,mpert,3,mpert*usepaw)=overlap contributions to the 2DTEs (PAW only)
221 : : !! eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
222 : : !! PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
223 : : !!
224 : : !! NOTES
225 : : !! We perform here the computation of
226 : : !! delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
227 : : !! see PRB 78, 035105 (2008), Eq. (42) [[cite:Audouze2008]]
228 : : !!
229 : : !! SOURCE
230 : :
231 : 1077 : subroutine dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,&
232 : 718 : & eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,&
233 : 359 : & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,&
234 : 1077 : & mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,mband_mem_rbz,ncpgr,nfftf,ngfftf,nhat,nhat1,&
235 : 718 : & nkpt_rbz,nkxc,npwarr,npwar1,nspden,nspinor,nsppol,nsym1,n3xccc,occkq,occ_rbz,&
236 [ - + - + ]: 359 : & paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,&
237 [ - + ]: 359 : & pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,&
238 [ - + ]: 359 : & ucvol,usecprj,usepaw,usevxctau,usexcnhat,useylmgr1,vectornd,vhartr1,vpsp1,vtrial,vtrial1,vxc,vxctau,&
239 [ - + - + : 359 : & 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 : 359 : type(gs_hamiltonian_type) :: gs_hamkq
325 : 359 : type(rf_hamiltonian_type) :: rf_hamkq
326 : 359 : type(MPI_type) :: mpi_enreg_seq
327 : : !arrays
328 : : integer :: ddkfil(3),my_spintab(2),nband_tmp(1),npwar1_tmp(1)
329 : 359 : integer,allocatable,target :: bands_treated_now_ndat(:,:),band_procs(:)
330 : 359 : integer,pointer :: bands_treated_now(:)
331 : 359 : integer,allocatable :: jpert1(:),jdir1(:),kg1_k(:,:),kg_k(:,:)
332 : 359 : integer,pointer :: my_atmtab(:)
333 : : real(dp) :: dum1(1,1),dum2(1,1),dum3(1,1),epawnst(2),kpoint(3),kpq(3)
334 : 359 : real(dp),allocatable :: vdotr(:),vdoti(:),vdot1r(:),vdot1i(:),vdot2r(:),vdot2i(:)
335 : : real(dp) :: sumelfd(2),symfact(3),tsec(2),ylmgr_dum(1,3,1)
336 : 359 : real(dp),allocatable :: buffer(:),ch1c(:,:,:,:),cs1c(:,:,:,:)
337 : 359 : real(dp),allocatable,target :: ch1c_tmp(:,:)
338 : 359 : real(dp),allocatable :: cs1c_tmp(:,:)
339 : 359 : real(dp),allocatable :: dcwavef(:,:)
340 : 359 : real(dp), contiguous, pointer :: cwave0(:,:),cwavef(:,:)!,dcwavef(:,:)
341 : 359 : real(dp),allocatable :: cg_ddk(:,:,:)
342 : 359 : real(dp),allocatable :: doccde_k(:),doccde_kq(:)
343 : 359 : real(dp),allocatable :: dnhat1(:,:),drhoaug1(:,:,:,:)
344 : 359 : real(dp),allocatable :: drhor1(:,:),drho1wfg(:,:),drho1wfr(:,:,:)
345 : 718 : real(dp),allocatable :: d2nl_elfd(:,:),dkinpw(:)
346 : 359 : real(dp),allocatable :: d2nl_k(:,:),d2ovl_drho(:,:,:,:,:),d2ovl_k(:,:)
347 : 359 : real(dp),allocatable :: eig_k(:),eig_kq(:),eig1_k(:)
348 : 359 : real(dp),allocatable,target :: e1kbfr_spin(:,:,:,:,:,:),ffnlk(:,:,:,:),ffnl1(:,:,:,:)
349 : 359 : real(dp),allocatable,target :: gh1(:,:),gs1(:,:),gvnlx1(:,:),gvnlx2(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:)
350 : 359 : real(dp),allocatable,target :: gvnlx1_tmp(:,:)
351 : 359 : real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:),rhotmp(:,:),rocceig(:,:)
352 : 359 : real(dp),allocatable :: vectornd_pac(:,:,:,:,:),vectornd_pac_idir(:,:,:,:),vlocal(:,:,:,:),vtrial_(:,:)
353 : 359 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
354 : 359 : real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:),vtmp1(:,:),vxc10(:,:)
355 : 359 : real(dp),allocatable,target :: work(:,:,:),e1kb_work(:,:,:,:)
356 : 359 : real(dp),pointer :: e1kbfr(:,:,:,:,:),e1kb_ptr(:,:,:,:)
357 : 359 : real(dp), contiguous, pointer :: ffnl1_idir1(:,:,:,:)
358 : 359 : real(dp),pointer :: vhartr01(:),vpsp1_idir1(:),xccc3d1_idir1(:)
359 : 359 : type(pawcprj_type),allocatable :: dcwaveprj(:,:)
360 : 359 : type(pawcprj_type),allocatable,target :: cwaveprj0(:,:)
361 : 359 : type(pawcprj_type),pointer :: cwaveprj0_idir1(:,:)
362 : 359 : type(paw_ij_type),allocatable :: paw_ij10(:,:)
363 : 359 : type(pawrhoij_type),target,allocatable :: pawdrhoij1(:,:)
364 : 359 : type(pawrhoij_type),pointer :: pawdrhoij1_unsym(:,:)
365 [ + + + + : 35182 : 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 : 359 : call timab(566,1,tsec)
375 : :
376 : : !Not valid for PrintBandByBand
377 [ - + ]: 359 : if (dtset%prtbbb/=0) then
378 : 0 : ABI_BUG('not yet valid for prtbbb/=0!')
379 : : end if
380 : :
381 : : !NCPP restrictions
382 [ - + ]: 359 : 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 [ + - ]: 359 : if (usepaw==1) then
395 : : ! Test on FFT grid sizes
396 [ - + ]: 359 : if (pawfgr%nfft/=nfftf) then
397 : 0 : ABI_BUG('PAW: wrong values for nfft, nfftf!')
398 : : end if
399 : : ! Test gradients of cprj
400 [ + + - + ]: 359 : if (ipert<=dtset%natom.and.ncpgr/=3) then
401 : 0 : ABI_BUG('PAW: wrong value of ncpgr for ipert<=natom!')
402 : : end if
403 [ + + + + : 359 : 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 [ + + - + ]: 359 : 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 [ + + + + : 359 : 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 [ + + + + ]: 359 : 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 : 359 : gpu_option=dtset%gpu_option
422 : 359 : ndat=dtset%bandpp
423 : 359 : master=0;me=mpi_enreg%me_kpt
424 : 359 : spaceworld=mpi_enreg%comm_cell
425 : 359 : paral_atom=(mpi_enreg%my_natom/=dtset%natom)
426 : 359 : my_comm_atom=mpi_enreg%comm_atom
427 : 359 : my_natom=mpi_enreg%my_natom
428 : 359 : my_nproc_band=mpi_enreg%nproc_band
429 : 359 : my_atmtab=>mpi_enreg%my_atmtab
430 [ + + ]: 1077 : my_spintab=mpi_enreg%my_isppoltab
431 [ + + + + ]: 1077 : my_nsppol=count(my_spintab==1)
432 : 359 : use_ompgpu=(gpu_option==ABI_GPU_OPENMP)
433 : :
434 : : !Fake MPI data to be used in sequential calls to parallel routines
435 : 359 : call initmpi_seq(mpi_enreg_seq)
436 : 359 : 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 : 359 : nkpt_me=0
442 [ + + ]: 751 : do isppol=1,nsppol
443 [ + + ]: 15522 : do ikpt=1,nkpt_rbz
444 : 14771 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
445 [ + + ]: 15163 : 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 : 359 : mcgq=mpw1*nspinor*mband_mem_rbz*mkqmem*nsppol
452 : 359 : mcprjq=nspinor*mband_mem_rbz*mkqmem*nsppol*usecprj
453 : :
454 [ + + - + : 16207 : 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 : 359 : ddkfil(:)=0
462 [ + + ]: 1436 : do idir1=1,3
463 : 1077 : ddkcase=idir1+dtset%natom*3
464 : 1077 : call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(idir1))
465 : 1077 : t_exist = file_exists(fiwfddk(idir1))
466 : :
467 [ + - ]: 1077 : if (.not. t_exist) then
468 : : ! Try netcdf file.
469 : 1077 : t_exist = file_exists(nctk_ncify(fiwfddk(idir1)))
470 [ + + ]: 1077 : 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 : 359 : 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 [ + + + + ]: 1034 : 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 [ + + ]: 359 : 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 : 141 : mpert1=1
511 [ - + ]: 141 : ABI_MALLOC(jpert1,(mpert1))
512 : 141 : jpert1(1)=dtset%natom+1
513 : : end if
514 : 359 : mdir1=3
515 [ - + ]: 359 : ABI_MALLOC(jdir1,(mdir1))
516 [ + + ]: 1436 : jdir1(1:3)= (/ (idir1,idir1=1,3) /)
517 : 359 : jdir1(1)=idir;jdir1(idir)=1
518 : :
519 : : !Index of strain perturbation, if any
520 [ + + ]: 359 : istr=idir;if (ipert==dtset%natom+4) istr=idir+3
521 : :
522 : : !Open ddk WF file(s) in sequential mode
523 [ + + ]: 359 : 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 [ + + + + : 25089 : d2nl(:,:,1:dtset%natom+4,idir,ipert)=zero
+ + ]
539 [ + - + + : 25089 : if (usepaw==1) d2ovl(:,:,1:dtset%natom+4,idir,ipert)=zero
+ + + + ]
540 : :
541 : : !Update list of computed matrix elements
542 [ + + ]: 1576 : do kpert1=1,mpert1
543 : 1217 : ipert1=jpert1(kpert1)
544 [ + + ]: 5227 : do kdir1=1,mdir1
545 : 3651 : 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 [ + + + + : 3651 : & (ipert1==dtset%natom+2.and.ddkfil(idir1)/=0).or.&
+ + + + +
- + + + +
+ + + - ]
549 : 1217 : & ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.ddkfil(idir1)/=0)) then
550 : 3498 : 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 : 359 : & usecprj=usecprj,nucdipmom=dtset%nucdipmom,gpu_option=gpu_option)
562 : 359 : has_vectornd = (with_vectornd .EQ. 1)
563 [ + + ]: 359 : if(has_vectornd) then
564 : : ! vlocal is needed when vectornd is present, for zora
565 [ - + + - : 342 : ABI_MALLOC(vlocal,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
566 [ - + + - : 228 : ABI_MALLOC(vtrial_,(nfftf,nspden))
- + + - +
- + - - +
- + ]
567 : : ! need a mutable version of vtrial
568 [ - + + + : 16185390 : vtrial_=vtrial
+ + ]
569 [ - + + - : 399 : ABI_MALLOC(vectornd_pac,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc,3))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + - + ]
570 [ - + + - : 285 : ABI_MALLOC(vectornd_pac_idir,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - - + -
+ - + ]
571 : : end if
572 : 359 : with_vxctau = ( usevxctau > 0 )
573 [ + + ]: 359 : 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 [ + - - + : 171973 : arg=maxval(occ_rbz)-minval(occ_rbz)
- - + + +
- - + - -
+ + ]
579 : 359 : qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14)
580 [ + + + + ]: 359 : is_metal=((dtset%occopt>=3.and.dtset%occopt<=8).or.(abs(arg)>tol8))
581 : 359 : is_metal_or_qne0=((is_metal).or.(qne0))
582 [ + - + - ]: 359 : nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
583 [ - + + - : 1795 : ABI_MALLOC(ch1c,(2,mband_mem_rbz,dtset%mband,nkpt_me))
- + + - -
+ + - + -
+ - - + -
+ ]
584 [ - + + - : 1077 : ABI_MALLOC(ch1c_tmp,(2,mband_mem_rbz))
+ - + - -
+ - + ]
585 [ + + + + : 2175766 : 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 : 359 : nzlmopt_ipert=0;nzlmopt_ipert1=0
591 [ + - ]: 359 : if (usepaw==1) then
592 [ - + + - : 1795 : ABI_MALLOC(d2ovl_drho,(2,3,mpert,3,mpert))
+ - - + +
- + - + -
- + - + ]
593 [ + + + + : 1094979 : d2ovl_drho=zero
+ + + + +
+ ]
594 [ + - ]: 359 : if (dtset%pawnzlm/=0) then
595 [ + + ]: 359 : nzlmopt_ipert=1;if (dtset%nstep<2) nzlmopt_ipert=-1
596 : 359 : nzlmopt_ipert1=-1
597 : : end if
598 [ + + ]: 359 : 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 [ + - ]: 359 : if (usepaw==1) then
607 : 359 : call paw_an_reset_flags(paw_an1)
608 : 359 : call paw_ij_reset_flags(paw_ij1,dijhartree=.true.)
609 : : end if
610 : :
611 : : !LOOP OVER PERTURBATION TYPES (j1)
612 [ + + ]: 1576 : do kpert1=1,mpert1
613 : 1217 : ipert1=jpert1(kpert1)
614 : :
615 : : ! Flag for use of DDK file
616 [ + + + - : 1217 : 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 : 1217 : elfd_fact=one
620 [ + + + + : 1217 : if ((ipert <=dtset%natom.or.ipert ==dtset%natom+3.or.ipert ==dtset%natom+4).and. &
+ + + + ]
621 : 1217 : (ipert1==dtset%natom+2)) elfd_fact=-one
622 [ + + + + : 1217 : if ((ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and. &
+ + + + ]
623 : 1217 : (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 : 1217 : has_dcwf=.false.;has_dcwf2=.false.;has_drho=.false.
628 [ + - ]: 1217 : if (usepaw==1) then
629 : 1217 : has_dcwf =(ipert1/=dtset%natom+2)
630 : 1217 : has_dcwf2=(ipert /=dtset%natom+2)
631 [ + + + + ]: 1217 : has_drho =(has_dcwf.and.ipert1/=dtset%natom+1)
632 : : end if
633 : :
634 : : ! Select which WF are needed
635 : 1217 : need_wfk=.true.
636 : 1217 : need_wf1=.true.
637 : :
638 : : ! Initialize data for NL 1st-order (j1) hamiltonian
639 : 1217 : 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 [ + + + - ]: 1217 : 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 [ + + + + ]: 1410 : d2lo(1:2,1:mdir1,ipert1,idir,ipert)=zero
806 : 141 : need_pawij10=.false.
807 [ - + ]: 141 : 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 : 1217 : iorder_cprj=0
812 : :
813 : : ! Allocate arrays used to accumulate density change due to overlap
814 [ + + ]: 1217 : 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 : 1217 : bdtot_index=0
843 : 1217 : bd2tot_index=0
844 : 1217 : ibg=0;icg=0
845 : 1217 : ibg1=0;icg1=0
846 : 1217 : 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 [ + + ]: 1217 : 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 [ + + ]: 2507 : do isppol=1,nsppol
872 : :
873 : 1290 : ikpt_me=0
874 : :
875 : : ! Rewind (k+G) data if needed
876 : 1290 : ikg=0;ikg1=0
877 : :
878 : : ! Continue to initialize the GS/RF Hamiltonian
879 : 1290 : call gs_hamkq%load_spin(isppol,with_nonlocal=.true.)
880 [ + + ]: 1290 : 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 [ + + ]: 1290 : if(has_vectornd) then
889 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
890 : 60 : & 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 : 60 : & dtset%nspden, gs_hamkq%nvloc, 3, pawfgr, mpi_enreg, vectornd,vectornd_pac)
893 : 60 : call gs_hamkq%load_spin(isppol, vlocal=vlocal,vectornd=vectornd_pac)
894 [ + + + + : 6586098 : vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,idir)
+ + + + ]
895 : 60 : 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 [ + + ]: 1290 : 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 [ + + + + : 35279712 : if (has_drho) drhoaug1(:,:,:,:)=zero
+ + + + +
+ ]
907 : :
908 : : ! LOOP OVER K-POINTS
909 [ + + ]: 32835 : do ikpt=1,nkpt_rbz
910 : :
911 : : ! Load dimensions for this k-point
912 : 31545 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
913 : 31545 : istwf_k=istwfk_rbz(ikpt)
914 : 31545 : npw_k=npwarr(ikpt)
915 : 31545 : npw1_k=npwar1(ikpt)
916 : 31545 : 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 [ + + ]: 31545 : 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 [ - + ]: 23129 : ABI_MALLOC(d2nl_k,(2,3))
927 [ - + ]: 23129 : ABI_MALLOC(d2ovl_k,(2,3))
928 [ - + + - : 69387 : ABI_MALLOC(eig_k,(nband_k))
+ - + - -
+ - + ]
929 [ - + + - : 46258 : ABI_MALLOC(eig_kq,(nband_k))
- + - + ]
930 [ + - - + : 69387 : ABI_MALLOC(eig1_k,(2*nband_k**2))
+ - + - +
- - + -
+ ]
931 [ - + + - : 46258 : ABI_MALLOC(occ_k,(nband_k))
- + - + ]
932 [ - + + - : 46258 : ABI_MALLOC(vdotr,(nband_k))
- + - + ]
933 [ - + + - : 46258 : ABI_MALLOC(vdoti,(nband_k))
- + - + ]
934 [ - + + - : 46258 : ABI_MALLOC(vdot1r,(nband_k))
- + - + ]
935 [ - + + - : 46258 : ABI_MALLOC(vdot1i,(nband_k))
- + - + ]
936 [ - + + - : 46258 : ABI_MALLOC(vdot2r,(nband_k))
- + - + ]
937 [ - + + - : 46258 : 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 [ + + + + ]: 231290 : d2nl_k(:,:)=zero
942 [ + + + + ]: 231290 : d2ovl_k(:,:)=zero
943 [ + + ]: 160529 : eig_k (:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
944 [ + + ]: 160529 : eig_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
945 [ + + ]: 2443409 : eig1_k(:)=eigen1(1+bd2tot_index:2*nband_k**2+bd2tot_index)
946 [ + + ]: 160529 : occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
947 [ + + + + ]: 160529 : nband_kocc=count(abs(occ_k(:))>tol8)
948 [ + + ]: 92516 : kpoint(:)=kpt_rbz(:,ikpt)
949 [ + + + + ]: 86960 : kpq(:)=kpoint(:);if (ipert1<dtset%natom+3) kpq(:)=kpq(:)+dtset%qptn(1:3)
950 : 23129 : wtk_k=wtk_rbz(ikpt)
951 : 23129 : need_ylmgr1=0;dimylmgr1=0
952 : 23129 : nkpg=0;nkpg1=0
953 : 23129 : ikpt_me=ikpt_me+1
954 [ + + ]: 23129 : 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 [ + + ]: 23129 : 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 [ - + + - : 69387 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + -
+ ]
1002 [ - + + - : 69387 : ABI_MALLOC(kg1_k,(3,npw1_k))
+ - + - -
+ - + -
+ ]
1003 [ - + + - : 92516 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
- + + - +
- + - - +
- + - + ]
1004 [ - + + - : 92516 : ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
- + + - +
- + - - +
- + - + ]
1005 [ - + + + : 23129 : if (psps%useylm==1.and.(need_ddk_file.or.ipert1==dtset%natom+1)) need_ylmgr1=1
+ + ]
1006 [ - + + + : 23129 : if (psps%useylm==1.and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)) need_ylmgr1=1
+ + ]
1007 : 23129 : dimylmgr1=max(useylmgr1,need_ylmgr1)
1008 [ - + + - : 92516 : 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 [ + + + + ]: 13696949 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
1012 [ + - ]: 23129 : if (psps%useylm==1) then
1013 [ + + ]: 121155 : do jj=1,psps%mpsang*psps%mpsang
1014 [ + + ]: 14551265 : 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 [ + + + + ]: 13696949 : kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
1020 [ + - ]: 23129 : if (psps%useylm==1) then
1021 [ + + ]: 121155 : do jj=1,psps%mpsang*psps%mpsang
1022 [ + + ]: 14551265 : ylm1_k(1:npw1_k,jj)=ylm1(1+ikg1:npw1_k+ikg1,jj)
1023 : : end do
1024 [ + + + + ]: 23129 : if (need_ylmgr1==1.and.useylmgr1/=0) then
1025 [ + + ]: 35985 : do jj=1,psps%mpsang*psps%mpsang
1026 [ + + ]: 124023 : do ia=1,3
1027 [ + + ]: 17815470 : 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 [ + - + + ]: 7391 : 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 [ + + ]: 23129 : nkpg=0;if(ipert1<=dtset%natom) nkpg=3*dtset%nloalg(3)
1044 [ - + + - : 92516 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ + + - +
- + + - +
- + - + ]
1045 [ + + ]: 23129 : 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 [ + + + - ]: 23129 : nkpg1=0;if(ipert1<=dtset%natom.or.need_ylmgr1==1) nkpg1=3*dtset%nloalg(3)
1051 [ - + + - : 92516 : ABI_MALLOC(kpg1_k,(npw1_k,nkpg1))
- + + - +
- + - - +
- + - + ]
1052 [ + - ]: 23129 : if (nkpg1>0) then
1053 : 23129 : call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k)
1054 : : end if
1055 : :
1056 : : ! Allocate kinetic contributions
1057 [ - + + - : 69387 : ABI_MALLOC(dkinpw,(npw_k))
+ - + - -
+ - + -
+ ]
1058 [ - + + - : 69387 : ABI_MALLOC(kinpw1,(npw1_k))
+ - + - -
+ - + -
+ ]
1059 [ + + ]: 3441584 : 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 : 23129 : 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 : 23129 : ider=0;idir0=0
1066 [ + + ]: 23129 : dimffnl=0;if (ipert1<=dtset%natom) dimffnl=1
1067 [ - + + - : 138774 : ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
+ + + - -
+ + - - +
+ - + - +
+ - + - +
- + ]
1068 [ + + ]: 23129 : 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 [ + + ]: 23129 : ider=1;if (ipert1<=dtset%natom) ider=0
1080 [ + + + + ]: 23129 : 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 [ + + ]: 21277 : dimffnl1=1;if (ider>=1) dimffnl1=2+2*psps%useylm
1085 [ + + + - ]: 21277 : idir0=0;if (ider>0.and.psps%useylm==1) idir0=4
1086 : : end if
1087 [ - + + - : 138774 : 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 : 23129 : & 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 [ + + ]: 23129 : if (ipert1<=dtset%natom) then
1098 : 15738 : ffnl1_idir1 => ffnl1(:,:,:,:)
1099 : 15738 : dimffnl1_idir1=dimffnl1
1100 : : else
1101 : 7391 : dimffnl1_idir1=1+ider
1102 [ - + + - : 44346 : ABI_MALLOC(ffnl1_idir1,(npw1_k,dimffnl1_idir1,psps%lmnmax,psps%ntypat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1103 [ - + ]: 7391 : ii=1;if (psps%useylm==0) ii=1+ider
1104 [ + + ]: 20842 : do itypat=1,psps%ntypat
1105 [ + + ]: 137436 : do ilmn=1,psps%lmnmax
1106 [ + - + + : 42771229 : 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 [ - + + - : 92516 : 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 : 23129 : & 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 [ + + + + ]: 115645 : if (size(ffnlk)>0) then
1122 : 15738 : call gs_hamkq%load_k(ffnl_k=ffnlk)
1123 : : else
1124 : 7391 : 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 : 23129 : & compute_gbound=.true.)
1132 [ + + ]: 23129 : 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 : 23129 : call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw)
1139 : :
1140 : : ! Allocate memory space for one band
1141 [ - + + - : 69387 : 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 : 23129 : nullify(cwaveprj0_idir1)
1146 [ + - ]: 23129 : if (usecprj==1) then
1147 [ - + + - : 167598 : ABI_MALLOC(cwaveprj0,(dtset%natom,nspinor*ndat))
- + + - +
- + - - +
- + - + +
+ + + ]
1148 : 23129 : call pawcprj_alloc(cwaveprj0,ncpgr,gs_hamkq%dimcprj)
1149 : : end if
1150 [ + + ]: 23129 : if (has_dcwf) then
1151 [ - + + - : 43154 : 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 [ - + + - : 69387 : ABI_MALLOC(band_procs, (nband_k))
+ - + - -
+ - + ]
1160 : : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,nband_k,&
1161 : 23129 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
1162 : :
1163 [ - + + - : 46258 : ABI_MALLOC(gvnlx1, (2,npw1_k*nspinor*ndat))
- + - + -
+ ]
1164 [ - + + - : 69387 : 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 [ + + ]: 23129 : 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 : 23129 : iband_me = 0
1178 [ - + - - : 160529 : 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 [ + + ]: 137400 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle
1183 : : end if
1184 : 135420 : iband_me = iband_me + ndat
1185 : :
1186 [ + - ]: 135420 : if(gpu_option==ABI_GPU_DISABLED) then
1187 [ + + + + ]: 3590340 : 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 : 135420 : & cg(:,1+(iband_me-ndat)*npw_k*nspinor+icg:iband_me*npw_k*nspinor+icg)
1196 [ + - ]: 135420 : if (usecprj==1) then
1197 [ + + ]: 270840 : 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 : 270840 : & 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 : 135420 : & 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 [ + + ]: 564809 : do kdir1=1,mdir1
1214 : 406260 : idir1=jdir1(kdir1)
1215 [ + + ]: 406260 : 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 [ + + + - ]: 406260 : if (ipert1==dtset%natom+2.and.ddkfil(idir1)==0) cycle
1219 : :
1220 [ + - ]: 406260 : if(gpu_option==ABI_GPU_DISABLED) then
1221 [ + + + + ]: 215650998 : 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 [ + + + - : 406260 : if (dimffnl1_idir1>=2.and.psps%useylm==1.and.ipert1>dtset%natom) then
+ - ]
1228 [ + + ]: 447684 : do itypat=1,psps%ntypat
1229 [ + + ]: 3306000 : do ilmn=1,psps%lmnmax
1230 [ + - + + : 1100353956 : 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 [ - + ]: 406260 : 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 [ + - ]: 406260 : if (usecprj==1) then
1246 : 406260 : cpopt=1
1247 : :
1248 : : ! === Atomic displ. perturbation
1249 [ + + ]: 406260 : 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 [ + + ]: 203364 : else if (ipert==dtset%natom+1) then
1273 : 94608 : 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 [ + + ]: 406260 : if (ipert1==dtset%natom+1) then
1331 : 94608 : 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 [ + + ]: 406260 : if (need_pawij10) rf_hamkq%e1kbfr => e1kbfr(:,:,:,:,idir1)
1338 : :
1339 : : ! Read DDK wave function (if ipert1=electric field)
1340 [ + + ]: 406260 : 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 : 375972 : 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 : 406260 : lambda=eig_k(iband);berryopt=0;optlocal=0
1370 [ + + + + ]: 406260 : optnl=0;if (ipert1/=dtset%natom+1.or.idir==idir1) optnl=1
1371 [ + + ]: 406260 : opt_gvnlx1=0;if (ipert1==dtset%natom+2) opt_gvnlx1=2
1372 [ + + ]: 406260 : sij_opt=-1;if (has_dcwf) sij_opt=1
1373 [ - + ]: 406260 : 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 : 406260 : & eig_k(iband:iband+ndat-1),mpi_enreg,ndat,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
1381 [ + + + + ]: 406260 : if (sij_opt==1.and.optnl==1) then
1382 [ + - ]: 312900 : if(gpu_option==ABI_GPU_DISABLED) then
1383 [ + + ]: 625800 : do idat=1,ndat
1384 : : gh1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor) = &
1385 [ + + + + ]: 164954010 : & 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 [ + + + + : 406260 : if (has_dcwf.or.(ipert==ipert1.and.idir==idir1.and.usepaw==1)) then
+ + + - ]
1415 : : ! note: gvnlx1 used as temporary space
1416 [ + + + + ]: 379892 : if (ipert==ipert1.and.idir==idir1) then
1417 : 70252 : option=0;
1418 [ + - ]: 70252 : if(gpu_option==ABI_GPU_DISABLED) then
1419 [ + - + - : 30848798 : 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 : 309640 : option=1;
1427 [ + - ]: 309640 : if(gpu_option==ABI_GPU_DISABLED) then
1428 [ + + + + ]: 168856984 : 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 [ + + ]: 759784 : do idat=1,ndat
1437 : :
1438 : 379892 : bands_treated_now => bands_treated_now_ndat(:,idat)
1439 [ + + ]: 3656152 : bands_treated_now = 0
1440 : 379892 : bands_treated_now(iband+idat-1) = 1
1441 : 379892 : call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr)
1442 [ + + ]: 4036044 : do iband_=1, nband_k
1443 [ + + ]: 3276260 : if (bands_treated_now(iband_) == 0) cycle
1444 : :
1445 : : ! distribute gvnlx1 to my subcomm
1446 [ + - ]: 385832 : if(gpu_option==ABI_GPU_DISABLED) then
1447 [ + + + + ]: 232999538 : 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 [ + + ]: 385832 : if (iband_ == iband+idat-1) then
1452 [ + - ]: 379892 : if(gpu_option==ABI_GPU_DISABLED) then
1453 [ + - + - : 200015422 : 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 [ + + ]: 385832 : if(xmpi_comm_size(mpi_enreg%comm_band)>1) then
1466 : 7920 : call xmpi_sum(gvnlx1_tmp,mpi_enreg%comm_band,ierr,use_omp_map=use_ompgpu)
1467 : : end if
1468 : :
1469 [ + + ]: 385832 : if (option == 1) then
1470 : : ! in case I need to reuse the ch1c (option 1) then load them here
1471 [ + - ]: 315328 : if(gpu_option==ABI_GPU_DISABLED) then
1472 [ + + + + ]: 8050840 : 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 : 385832 : & 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 [ + + ]: 385832 : if(xmpi_comm_size(mpi_enreg%comm_band)>1) then
1488 : 7920 : 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 [ + + ]: 385832 : 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 [ + - ]: 379892 : if(gpu_option==ABI_GPU_DISABLED) then
1498 : : gvnlx1(:,1+(idat-1)*npw1_k*nspinor:idat*npw1_k*nspinor) = &
1499 [ + + + + ]: 199635530 : 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 [ + + ]: 765724 : 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 [ + - ]: 70504 : if(gpu_option==ABI_GPU_DISABLED) then
1516 [ + + + + ]: 2163772 : 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 [ + + ]: 379892 : if (has_dcwf) then
1528 [ + + + + ]: 375972 : if (ipert==ipert1.and.idir==idir1) then
1529 [ + - ]: 66332 : if(gpu_option==ABI_GPU_DISABLED) then
1530 [ - + + + : 29027230 : 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 : 375972 : do_scprod=0
1537 [ + + ]: 751944 : do idat=1,ndat
1538 [ + + ]: 751944 : if (abs(occ_k(iband+idat-1))>tol8) do_scprod=1
1539 : : end do !idat
1540 [ + + ]: 375972 : 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 : 309450 : & 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 [ + + ]: 618900 : do idat=1,ndat
1553 [ + - ]: 618900 : if (abs(occ_k(iband+idat-1))>tol8) then
1554 : 309450 : d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband+idat-1)*elfd_fact*vdotr(idat)
1555 : 309450 : 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 [ + + ]: 402340 : 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 [ + + ]: 812520 : do idat=1,ndat
1595 : :
1596 : 406260 : 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 [ + + ]: 812520 : 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 [ + + ]: 406260 : 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 [ + + ]: 94608 : if (idir==idir1) then
1719 [ + - ]: 31536 : if(gpu_option==ABI_GPU_DISABLED) then
1720 [ + - + + : 33145236 : 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 [ - + - - ]: 63072 : 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 [ + - ]: 63072 : if(gpu_option==ABI_GPU_DISABLED) then
1737 [ + + + + ]: 33176772 : 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 : 94608 : do_scprod=0
1744 [ + + ]: 189216 : do idat=1,ndat
1745 [ + + ]: 189216 : if (abs(occ_k(iband+idat-1))>tol8) then
1746 : 63522 : do_scprod=1
1747 : : end if
1748 : : end do
1749 [ + + ]: 94608 : 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 : 63522 : & 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 : 63522 : & 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 : 63522 : & 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 [ + + ]: 127044 : do idat=1,ndat
1766 [ + - ]: 127044 : if (abs(occ_k(iband+idat-1))>tol8) then
1767 : 63522 : vdotr(idat)=vdotr(idat)-(vdot1r(idat)*vdot2r(idat)-vdot1i(idat)*vdot2i(idat))
1768 : 63522 : vdoti(idat)=vdoti(idat)-(vdot1r(idat)*vdot2i(idat)+vdot1i(idat)*vdot2r(idat))
1769 : 63522 : d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband+idat-1)*vdotr(idat)/(nband_kocc*two)
1770 : 63522 : 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 [ + + ]: 406260 : 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 [ - + + + ]: 543660 : 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 [ - + ]: 23129 : ABI_FREE(band_procs)
1828 : :
1829 : : ! Accumulate contribution of this k-point
1830 [ + + + + ]: 231290 : d2nl (:,:,ipert1,idir,ipert)=d2nl (:,:,ipert1,idir,ipert)+d2nl_k (:,:)
1831 [ + - + + : 231290 : 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 [ + + ]: 23129 : 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 [ - + ]: 23129 : ABI_FREE(gvnlx1)
1849 [ - + ]: 23129 : ABI_FREE(gvnlx1_tmp)
1850 [ - + ]: 23129 : ABI_FREE(gh1)
1851 [ - + ]: 23129 : ABI_FREE(gs1)
1852 : : if (need_wfk) then
1853 : 23129 : nullify(cwave0)
1854 : : end if
1855 : : if (need_wf1) then
1856 : 23129 : nullify(cwavef)
1857 : : end if
1858 [ - + ]: 23129 : ABI_FREE(kg_k)
1859 [ - + ]: 23129 : ABI_FREE(kg1_k)
1860 [ - + ]: 23129 : ABI_FREE(ylm_k)
1861 [ - + ]: 23129 : ABI_FREE(ylm1_k)
1862 [ - + ]: 23129 : ABI_FREE(ylmgr1_k)
1863 [ - + ]: 23129 : ABI_FREE(kpg_k)
1864 [ - + ]: 23129 : ABI_FREE(kpg1_k)
1865 [ - + ]: 23129 : ABI_FREE(d2nl_k)
1866 [ - + ]: 23129 : ABI_FREE(d2ovl_k)
1867 [ - + ]: 23129 : ABI_FREE(eig_k)
1868 [ - + ]: 23129 : ABI_FREE(eig_kq)
1869 [ - + ]: 23129 : ABI_FREE(eig1_k)
1870 [ - + ]: 23129 : 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 [ - + ]: 23129 : ABI_FREE(vdotr)
1875 [ - + ]: 23129 : ABI_FREE(vdoti)
1876 [ - + ]: 23129 : ABI_FREE(vdot1r)
1877 [ - + ]: 23129 : ABI_FREE(vdot1i)
1878 [ - + ]: 23129 : ABI_FREE(vdot2r)
1879 [ - + ]: 23129 : ABI_FREE(vdot2i)
1880 [ + + ]: 23129 : 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 [ - + ]: 23129 : ABI_FREE(dkinpw)
1887 [ - + ]: 23129 : 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 [ - + ]: 23129 : ABI_FREE(ph3d)
1896 [ + + ]: 23129 : if (allocated(ph3d1)) then
1897 : 11360 : ABI_FREE(ph3d1)
1898 : : end if
1899 [ - + ]: 23129 : ABI_FREE(ffnlk)
1900 [ - + ]: 23129 : ABI_FREE(ffnl1)
1901 [ + + ]: 23129 : 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 [ - + ]: 7391 : ABI_FREE(ffnl1_idir1)
1906 : : end if
1907 : 23129 : nullify(ffnl1_idir1)
1908 [ + - ]: 23129 : if (usecprj==1) then
1909 : 23129 : call pawcprj_free(cwaveprj0)
1910 [ + - + + : 74680 : ABI_FREE(cwaveprj0)
- + - + -
+ ]
1911 : : end if
1912 : 23129 : nullify(cwaveprj0_idir1)
1913 : : ! Shift arrays
1914 : 23129 : bdtot_index=bdtot_index+nband_k
1915 : 23129 : bd2tot_index=bd2tot_index+2*nband_k**2
1916 [ + - ]: 23129 : if (mkmem/=0) then
1917 : 23129 : ibg=ibg+nspinor*nband_me
1918 : 23129 : icg=icg+npw_k*nspinor*nband_me
1919 : 23129 : ikg=ikg+npw_k
1920 : : end if
1921 [ + - ]: 23129 : if (mkqmem/=0) then
1922 : 23129 : ibgq=ibgq+nspinor*nband_me
1923 : 23129 : icgq=icgq+npw1_k*nspinor*nband_me
1924 : : end if
1925 [ + - ]: 24419 : if (mk1mem/=0) then
1926 : 23129 : ibg1=ibg1+nspinor*nband_me
1927 : 23129 : icg1=icg1+npw1_k*nspinor*nband_me
1928 : 23129 : 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 [ + + ]: 2507 : 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 : 1217 : call rf_hamkq%free()
1950 [ + + ]: 1217 : if (has_vectornd) then
1951 [ - + ]: 57 : ABI_FREE(vectornd_pac_idir)
1952 : : end if
1953 [ + + ]: 1217 : if (allocated(vlocal)) then
1954 : 57 : ABI_FREE(vlocal)
1955 : : end if
1956 [ + + ]: 1217 : if (allocated(vtrial_)) then
1957 : 57 : ABI_FREE(vtrial_)
1958 : : end if
1959 [ + + ]: 1217 : if (has_drho) then
1960 [ - + ]: 942 : ABI_FREE(drhoaug1)
1961 : : end if
1962 [ + + ]: 1217 : 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 [ + - + + : 8030 : ABI_FREE(paw_ij10)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
1970 : :
1971 : : ! In case of parallelism, sum 1st-order density and occupation matrix over processors
1972 [ + + ]: 1217 : 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 : 359 : 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 [ - + ]: 359 : ABI_FREE(ch1c)
2094 [ - + ]: 359 : ABI_FREE(ch1c_tmp)
2095 [ + - + + ]: 359 : if (usepaw==1.and.is_metal_or_qne0) then
2096 [ - + ]: 250 : ABI_FREE(cs1c)
2097 [ - + ]: 250 : ABI_FREE(cs1c_tmp)
2098 : : end if
2099 : 359 : call gs_hamkq%free()
2100 [ + + ]: 359 : if(allocated(vectornd_pac)) then
2101 : 57 : ABI_FREE(vectornd_pac)
2102 : : end if
2103 [ + + ]: 359 : 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 : 359 : call timab(161,1,tsec)
2110 : 359 : call xmpi_barrier(spaceworld)
2111 : 359 : call timab(161,2,tsec)
2112 [ + - - + : 1077 : ABI_MALLOC(buffer,(6*mpert*(1+usepaw)))
+ - + - +
- - + - +
- + ]
2113 [ + + ]: 718 : buffer(1:6*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/6*mpert/))
2114 [ + - + + ]: 718 : if (usepaw==1) buffer(6*mpert+1:6*mpert+6*mpert)=reshape(d2ovl(:,:,:,idir,ipert),(/6*mpert/))
2115 : 359 : call timab(48,1,tsec)
2116 : 359 : call xmpi_sum(buffer,6*mpert*(1+usepaw),spaceworld,ierr)
2117 : 359 : call timab(48,2,tsec)
2118 [ + + ]: 1436 : d2nl (:,:,:,idir,ipert)=reshape(buffer(1:6*mpert),(/2,3,mpert/))
2119 [ + - + + ]: 1436 : if (usepaw==1) d2ovl(:,:,:,idir,ipert)=reshape(buffer(6*mpert+1:6*mpert+6*mpert),(/2,3,mpert/))
2120 [ - + ]: 359 : ABI_FREE(buffer)
2121 : : end if
2122 : :
2123 : : !Build complete d2ovl matrix
2124 [ + - + + : 35859 : if (usepaw==1) d2ovl(:,:,:,idir,ipert)=d2ovl(:,:,:,idir,ipert)+d2ovl_drho(:,:,:,idir,ipert)
+ + + + ]
2125 : :
2126 [ + - ]: 359 : if (usepaw==1) then
2127 [ - + ]: 359 : ABI_FREE(d2ovl_drho)
2128 : : end if
2129 : :
2130 : : !Close the ddk WF files
2131 [ + + ]: 359 : 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 [ - + ]: 359 : ABI_FREE(jpert1)
2139 : 359 : ABI_FREE(jdir1)
2140 : :
2141 : : !Symmetrize the phonons contributions, as was needed for the forces in a GS calculation
2142 [ - + + - : 1077 : ABI_MALLOC(work,(2,3,dtset%natom))
+ - + - -
+ - + ]
2143 [ + + ]: 1396 : do ipert1=1,dtset%natom
2144 [ + + ]: 4507 : do idir1=1,3
2145 [ + + ]: 10370 : work(:,idir1,ipert1)=d2nl(:,idir1,ipert1,idir,ipert)
2146 : : end do
2147 : : end do
2148 : 359 : call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1)
2149 [ + - ]: 359 : if (usepaw==1) then
2150 [ + + ]: 1396 : do ipert1=1,dtset%natom
2151 [ + + ]: 4507 : do idir1=1,3
2152 [ + + ]: 10370 : work(:,idir1,ipert1)=d2ovl(:,idir1,ipert1,idir,ipert)
2153 : : end do
2154 : : end do
2155 : 359 : call dfpt_sygra(dtset%natom,d2ovl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1)
2156 : : end if
2157 [ - + ]: 359 : 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 [ + + + + ]: 359 : 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 [ + + + + ]: 2583 : d2nl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero
2167 [ + - + + : 2583 : 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 [ + + ]: 359 : 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 [ + + ]: 359 : 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 [ + - ]: 359 : eovl1=zero;if (usepaw==1) eovl1=d2ovl(1,idir,ipert,idir,ipert)
2241 : :
2242 [ - + ]: 359 : 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 : 359 : call destroy_mpi_enreg(mpi_enreg_seq)
2249 : 359 : call timab(566,2,tsec)
2250 : :
2251 : : ABI_NVTX_END_RANGE()
2252 : :
2253 : : DBG_EXIT("COLL")
2254 : :
2255 [ - + - + : 1436 : 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 : : !!***
|
