Line data Source code
1 : !!****m* ABINIT/m_forstr
2 : !! NAME
3 : !! m_forstr
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, AF, AR, MB, MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_forstr
26 :
27 : use defs_basis
28 : use defs_wvltypes
29 : use m_abicore
30 : use m_efield
31 : use m_errors
32 : use m_xmpi
33 : use m_fock
34 : use m_cgtools
35 : use m_xcdata
36 : use m_dtset
37 : use m_extfpmd
38 : use m_ompgpu_utils
39 : use m_xg
40 : use m_xg_nonlop
41 : use m_xgTransposer
42 :
43 : use defs_datatypes, only : pseudopotential_type
44 : use defs_abitypes, only : MPI_type
45 : use m_time, only : timab
46 : use m_geometry, only : xred2xcart, metric, stresssym
47 : use m_energies, only : energies_type
48 : use m_pawang, only : pawang_type
49 : use m_pawrad, only : pawrad_type
50 : use m_pawtab, only : pawtab_type
51 : use m_paw_ij, only : paw_ij_type
52 : use m_pawfgrtab, only : pawfgrtab_type
53 : use m_pawrhoij, only : pawrhoij_type
54 : use m_pawfgr, only : pawfgr_type
55 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_reorder, pawcprj_getdim
56 : use m_paw_dfpt, only : pawgrnl
57 : use libxc_functionals, only : libxc_functionals_is_hybrid
58 : use m_stress, only : stress
59 : use m_forces, only : forces
60 : use m_initylmg, only : initylmg
61 : use m_xchybrid, only : xchybrid_ncpp_cc
62 : use m_kg, only : mkkpg
63 : use m_hamiltonian, only : gs_hamiltonian_type, gs_hamiltonian_type, gspot_transgrid_and_pack, K_H_K, KPRIME_H_KPRIME !,K_H_KPRIME
64 : use m_electronpositron, only : electronpositron_type, electronpositron_calctype
65 : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_type, prep_bandfft_tabs, &
66 : & bandfft_kpt_savetabs, bandfft_kpt_restoretabs
67 : use m_spacepar, only : meanvalue_g, hartre
68 : use m_mkffnl, only : mkffnl
69 : use m_mpinfo, only : proc_distrb_cycle
70 : use m_nonlop, only : nonlop
71 : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, reset_gemm_nonlop, gemm_nonlop_use_gemm, &
72 : gemm_nonlop_block_size, gemm_nonlop_is_distributed
73 : use m_common, only : get_gemm_nonlop_ompgpu_blocksize
74 : use m_fock_getghc, only : fock_getghc
75 : use m_prep_kgb, only : prep_nonlop
76 : use m_paw_nhat, only : pawmknhat
77 : use m_rhotoxc, only : rhotoxc
78 : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll
79 : use m_cgprj, only : ctocprj,xg_cprj_copy,CPRJ_TO_XG
80 : use m_psolver, only : psolver_hartree
81 : use m_wvl_psi, only : wvl_nl_gradient
82 : use m_fft, only : fourdp,fourwf
83 : use m_alloc_hamilt_gpu, only : hamilt_gpu_nfft_blocks
84 : use, intrinsic :: iso_c_binding, only : c_loc,c_f_pointer,c_double,c_size_t
85 :
86 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
87 : use gator_mod
88 : use m_gpu_toolbox, only : CPU_DEVICE_ID, gpu_device_synchronize, gpu_data_prefetch_async
89 : #endif
90 :
91 : #if defined(HAVE_GPU_MARKERS)
92 : use m_nvtx_data
93 : #endif
94 :
95 : implicit none
96 :
97 : private
98 : !!***
99 :
100 : public :: forstr ! Drives the computation of forces and/or stress tensor
101 : public :: nres2vres ! Convert a density residual into a potential residual
102 : !!***
103 :
104 : contains
105 : !!***
106 :
107 : !----------------------------------------------------------------------
108 :
109 : !!****f* ABINIT/forstr
110 : !! NAME
111 : !! forstr
112 : !!
113 : !! FUNCTION
114 : !! Drives the computation of forces and/or stress tensor
115 : !!
116 : !! INPUTS
117 : !! atindx1(natom)=index table for atoms, inverse of atindx
118 : !! cg(2,mcg)=wavefunctions (may be read from disk instead of input)
119 : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
120 : !! dtefield <type(efield_type)> = variables related to Berry phase
121 : !! dtset <type(dataset_type)>=all input variables in this dataset
122 : !! | berryopt = 4: electric field is on -> add the contribution of the
123 : !! | - \Omega E.P term to the total energy
124 : !! | /= 4: electric field is off
125 : !! | from Etot(npw) data (at fixed geometry), used for making
126 : !! | Pulay correction to stress tensor (hartree). Should be <=0.
127 : !! | ecut=cut-off energy for plane wave basis sphere (Ha)
128 : !! | ecutsm=smearing energy for plane wave kinetic energy (Ha)
129 : !! | effmass_free=effective mass for electrons (1. in common case)
130 : !! | efield = cartesian coordinates of the electric field in atomic units
131 : !! | ionmov=governs the movement of atoms (see help file)
132 : !! | densfor_pred=governs the mixed electronic-atomic part of the preconditioner
133 : !! | istwfk(nkpt)=input option parameter that describes the storage of wfs
134 : !! | kptns(3,nkpt)=reduced coordinates of k points in Brillouin zone
135 : !! | mband=maximum number of bands
136 : !! | mgfft=maximum size of 1D FFTs
137 : !! | mkmem=maximum number of k points in core memory
138 : !! | mpw = maximum number of plane waves
139 : !! | natom=number of atoms in cell
140 : !! | nband(nkpt*nsppol)=number of bands to be included in summation at each k point
141 : !! | nfft=(effective) number of FFT grid points (for this processor)
142 : !! | ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
143 : !! | nkpt=number of k points in Brillouin zone
144 : !! | nloalg(3)=governs the choice of the algorithm for non-local operator.
145 : !! | nspden=number of spin-density components
146 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
147 : !! | pawprtvol=control print volume and debugging output for PAW
148 : !! | prtvol=integer controlling volume of printed output
149 : !! | symafm(nsym)=(anti)ferromagnetic part of symmetry operations
150 : !! | tfkinfunc=1 if use of Thomas-Fermi kinetic functional
151 : !! | =2 if use of recursion method
152 : !! | typat(natom)=type integer for each atom in cell
153 : !! | wtk(nkpt)=weights associated with various k points
154 : !! | nsym=number of symmetries in space group
155 : !! energies <type(energies_type)>=all part of total energy.
156 : !! | e_localpsp(IN)=local psp energy (hartree)
157 : !! | e_hartree(IN)=Hartree part of total energy (hartree units)
158 : !! | e_corepsp(IN)=psp core-core energy
159 : !! | e_kinetic(IN)=kinetic energy part of total energy.
160 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
161 : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
162 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
163 : !! grchempottn(3,natom)=d(E_chemical potential)/d(xred) (hartree)
164 : !! grcondft(3,natom)=d(E_constrainedDFT)/d(xred) (hartree)
165 : !! grewtn(3,natom)=d(Ewald)/d(xred) (hartree)
166 : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree)
167 : !! gsqcut=cutoff value on G**2 for (large) sphere inside FFT box.
168 : !! gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)
169 : !! indsym(4,nsym,natom)=index showing transformation of atom labels
170 : !! under symmetry operations (computed in symatm)
171 : !! kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere
172 : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
173 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
174 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
175 : !! mgfftf= -PAW ONLY- maximum size of 1D FFTs for the fine grid
176 : !! (mgfftf=mgfft for norm-conserving potential runs)
177 : !! mpi_enreg=information about MPI parallelization
178 : !! my_natom=number of atoms treated by current processor
179 : !! n3xccc=dimension of the xccc3d array (0 or nfftf).
180 : !! nattyp(ntypat)=number of atoms of each type
181 : !! nfftf= -PAW ONLY- number of FFT grid points for the fine grid
182 : !! (nfftf=nfft for norm-conserving potential runs)
183 : !! ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
184 : !! (ngfftf=ngfft for norm-conserving potential runs)
185 : !! ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
186 : !! nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density
187 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
188 : !! npwarr(nkpt)=number of planewaves in basis and on boundary for each k
189 : !! ntypat=number of types of atoms
190 : !! nvresid(nfftf,nspden)=array for the residual of the density/potential
191 : !! occ(mband*nkpt*nsppol)=occupancies of bands at various k points
192 : !! optfor=1 if computation of forces is required
193 : !! optres=0 if the potential residual has to be used for forces corrections
194 : !! =1 if the density residual has to be used for forces corrections
195 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
196 : !! pawang <type(pawang_type)>=paw angular mesh and related data
197 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
198 : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
199 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
200 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
201 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
202 : !! ph1df(2,3*(2*mgfftf+1)*natom)=-PAW only- 1-dim structure factor phases for the fine grid
203 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
204 : !! | mpsang= 1+maximum angular momentum
205 : !! rhog(2,nfftf)=Fourier transform of charge density (bohr^-3)
206 : !! rhor(nfftf,nspden)=array for electron density in electrons/bohr**3.
207 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
208 : !! strscondft(6)=cDFT correction to stress
209 : !! strsxc(6)=xc correction to stress
210 : !! stress_needed=1 if computation of stress tensor is required
211 : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
212 : !! ucvol=unit cell volume in bohr**3
213 : !! usecprj=1 if cprj datastructure is stored in memory
214 : !! usevxctau=1 if kinetic energy density contribution has to be included (mGGA)
215 : !! vhartr(nfftf)=array for holding Hartree potential
216 : !! vpsp(nfftf)=array for holding local psp
217 : !! vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space
218 : !! vxctau(nfftf,nspden,4*usevxctau)=(only for meta-GGA): derivative of XC energy density
219 : !! wrt kinetic energy density (depsxcdtau)
220 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
221 : !! xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction
222 : !! xred(3,natom)=reduced dimensionless atomic coordinates
223 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
224 : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
225 : !!
226 : !! OUTPUT
227 : !! ==== if (optfor==1) ====
228 : !! diffor=maximal absolute value of changes in the components of
229 : !! force between the input and the output.
230 : !! favg(3)=mean of the forces before correction for translational symmetry
231 : !! fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
232 : !! at input, previous value of forces,
233 : !! at output, new value.
234 : !! Note : unlike gred, this array has been corrected by enforcing
235 : !! the translational symmetry, namely that the sum of force
236 : !! on all atoms is zero.
237 : !! forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
238 : !! gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
239 : !! gresid(3,natom)=forces due to the residual of the density/potential
240 : !! grhf(3,natom)=Hellman-Feynman derivatives of the total energy
241 : !! grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used
242 : !! maxfor=maximal absolute value of the output array force.
243 : !! synlgr(3,natom)=symmetrized gradients of energy due to nonlocal contributions
244 : !! ==== if (stress_needed==1) ====
245 : !! strten(6)=components of the stress tensor (hartree/bohr^3) for the
246 : !! 6 unique components of this symmetric 3x3 tensor:
247 : !! Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
248 : !!
249 : !! SIDE EFFECTS
250 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
251 : !! forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
252 : !! ===== if psps%usepaw==1
253 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
254 : !! (gradients of rhoij for each atom with respect to atomic positions are computed here)
255 : !! wvl <type(wvl_data)>=all wavelets data.
256 : !!
257 : !! NOTES
258 : !! Be careful to the meaning of nfft (size of FFT grids):
259 : !! - In case of norm-conserving calculations the FFT grid is the usual FFT grid.
260 : !! - In case of PAW calculations:
261 : !! Two FFT grids are used; one with nfft points (coarse grid) for
262 : !! the computation of wave functions ; one with nfftf points
263 : !! (fine grid) for the computation of total density.
264 : !!
265 : !! SOURCE
266 :
267 6333 : subroutine forstr(atindx1,cg,cprj,diffor,dtefield,dtset,eigen,electronpositron,energies,favg,fcart,fock,&
268 6333 : & forold,gred,grchempottn,grcondft,gresid,grewtn,grhf,grvdw,grxc,gsqcut,extfpmd,indsym,&
269 6333 : & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,&
270 6333 : & nfftf,ngfftf,ngrvdw,nhat,nkxc,npwarr,&
271 6333 : & ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,&
272 6333 : & pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,psps,rhog,rhor,rprimd,stress_needed,&
273 6333 : & strscondft,strsxc,strten,symrec,synlgr,ucvol,usecprj,usevxctau,vhartr,vpsp,&
274 6333 : & vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero,xg_nonlop)
275 :
276 : !Arguments ------------------------------------
277 : !scalars
278 : integer,intent(in) :: mcg,mcprj,mgfftf,my_natom,n3xccc,nfftf,ngrvdw,nkxc,ntypat,optfor,optres
279 : integer,intent(in) :: stress_needed,usecprj,usevxctau
280 : real(dp),intent(in) :: gsqcut,qvpotzero,ucvol
281 : real(dp),intent(inout) :: diffor,maxfor
282 : type(electronpositron_type),pointer :: electronpositron
283 : type(MPI_type),intent(inout) :: mpi_enreg
284 : type(efield_type),intent(in) :: dtefield
285 : type(dataset_type),intent(in) :: dtset
286 : type(energies_type),intent(in) :: energies
287 : type(extfpmd_type),pointer,intent(inout) :: extfpmd
288 : type(pawang_type),intent(in) :: pawang
289 : type(pawfgr_type),intent(in) :: pawfgr
290 : type(pseudopotential_type),intent(in) :: psps
291 : type(wvl_data),intent(inout) :: wvl
292 : type(fock_type),pointer, intent(inout) :: fock
293 : type(xg_nonlop_t), intent(inout) :: xg_nonlop
294 : !arrays
295 : integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
296 : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfftf(18)
297 : integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym)
298 : real(dp),intent(in) :: cg(2,mcg)
299 : real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
300 : real(dp),intent(in) :: grchempottn(3,dtset%natom),grcondft(3,dtset%natom),grewtn(3,dtset%natom)
301 : real(dp),intent(in) :: grvdw(3,ngrvdw),kxc(dtset%nfft,nkxc)
302 : real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
303 : real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
304 : real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom)
305 : real(dp),intent(in) :: rhog(2,nfftf),strscondft(6),strsxc(6),vhartr(nfftf)
306 : real(dp),intent(in) :: vpsp(nfftf),vxc(nfftf,dtset%nspden),vxctau(nfftf,dtset%nspden,4*usevxctau)
307 : real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
308 : real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
309 : real(dp),intent(inout) :: forold(3,dtset%natom)
310 : real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw),rhor(nfftf,dtset%nspden),rprimd(3,3)
311 : real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom)
312 : real(dp),intent(inout),target :: nvresid(nfftf,dtset%nspden)
313 : real(dp),intent(out) :: favg(3)
314 : real(dp),intent(inout) :: fcart(3,dtset%natom),gred(3,dtset%natom)
315 : real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
316 : real(dp),intent(inout) :: grxc(3,dtset%natom),strten(6),synlgr(3,dtset%natom)
317 : type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj)
318 : type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
319 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
320 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
321 : type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
322 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
323 :
324 : !Local variables-------------------------------
325 : !scalars
326 : integer :: comm_grid,ifft,ispden,ncpgr,occopt_,optgr,optgr2,option,optnc,optstr,optstr2,iorder_cprj,ctocprj_choice
327 : integer :: idir,iatom,unpaw,mcgbz,usexg
328 6333 : integer,allocatable :: dimcprj(:)
329 : real(dp) ::dum,dum1,dum2,ucvol_
330 : logical :: apply_residual
331 : !arrays
332 : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
333 : real(dp) :: kinstr(6),mggastr(6),nlstr(6),tsec(2),strdum(6),gmet(3,3),gprimd(3,3),rmet(3,3)
334 : real(dp) :: dummy(0)
335 6333 : real(dp),allocatable :: grnl(:),vlocal(:,:),vxc_hf(:,:),xcart(:,:),ylmbz(:,:),ylmgrbz(:,:,:)
336 6333 : real(dp), contiguous, pointer :: resid(:,:)
337 : ! *************************************************************************
338 :
339 6333 : call timab(910,1,tsec)
340 : ABI_NVTX_START_RANGE(NVTX_FORSTR)
341 6333 : call timab(911,1,tsec)
342 :
343 : !Do nothing if nothing is required
344 6333 : if (optfor==0.and.stress_needed==0) return
345 :
346 : !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
347 6333 : if (dtset%usewvl==0) then
348 6333 : if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
349 0 : ABI_BUG('Wrong values for nfft, nfftf!')
350 : end if
351 6333 : if ((psps%usepaw==1.and.pawfgr%mgfft/=mgfftf).or.(psps%usepaw==0.and.dtset%mgfft/=mgfftf)) then
352 0 : ABI_BUG('Wrong values for mgfft, mgfftf!')
353 : end if
354 : end if
355 :
356 : !==========================================================================
357 : !Here compute terms common to forces and stresses
358 : !==========================================================================
359 :
360 : !output only if (optfor==1) but we have to allocate it
361 18999 : ABI_MALLOC(grnl,(3*dtset%natom*optfor))
362 40632 : grnl(:)=zero
363 :
364 : !Compute nonlocal psp + potential Fock ACE parts of forces and stress tensor
365 : !-involves summations over wavefunctions at all k points
366 6333 : if (dtset%tfkinfunc>0.and.stress_needed==1) then
367 301 : kinstr(1:3)=-two/three*energies%e_kinetic/ucvol ; kinstr(4:6)=zero
368 43 : mggastr(1:6)=zero ;nlstr(1:6)=zero
369 6290 : else if (dtset%usewvl==0) then
370 6290 : occopt_=0 ! This means that occ are now fixed
371 6290 : if(dtset%usefock==1 .and. associated(fock)) then
372 : ! if((dtset%optstress/=0).and.(psps%usepaw==1)) then
373 122 : if((psps%usepaw==1).and.((dtset%optstress/=0).or.(dtset%optforces==2))) then
374 26 : if(dtset%optstress==0) then
375 1 : ctocprj_choice=2
376 1 : ncpgr=3
377 : end if
378 26 : if(dtset%optstress/=0) then
379 25 : ctocprj_choice=20*optfor+3*dtset%optstress
380 25 : ncpgr=6*dtset%optstress+3*optfor
381 : end if
382 26 : if (allocated(fock%fock_BZ%cwaveocc_prj)) then
383 26 : call pawcprj_free(fock%fock_BZ%cwaveocc_prj)
384 1238 : ABI_FREE(fock%fock_BZ%cwaveocc_prj)
385 1922 : ABI_MALLOC(fock%fock_BZ%cwaveocc_prj,(dtset%natom,fock%fock_BZ%mcprj))
386 78 : ABI_MALLOC(dimcprj,(dtset%natom))
387 26 : call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
388 26 : call pawcprj_alloc(fock%fock_BZ%cwaveocc_prj,ncpgr,dimcprj)
389 26 : ABI_FREE(dimcprj)
390 : end if
391 26 : iatom=-1;idir=0;iorder_cprj=0;unpaw=26
392 26 : call metric(gmet,gprimd,-1,rmet,rprimd,dum)
393 26 : if (fock%fock_BZ%mkpt/=dtset%mkmem.or.(fock%fock_BZ%mpi_enreg%paral_hf ==1)) then
394 40 : ABI_MALLOC(ylmbz,(dtset%mpw*fock%fock_BZ%mkpt,psps%mpsang*psps%mpsang*psps%useylm))
395 40 : ABI_MALLOC(ylmgrbz,(dtset%mpw*fock%fock_BZ%mkpt,3,psps%mpsang*psps%mpsang*psps%useylm))
396 10 : option=1; mcgbz=dtset%mpw*fock%fock_BZ%mkptband*fock%fock_common%my_nsppol
397 : call initylmg(gprimd,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,&
398 : & psps%mpsang,dtset%mpw,fock%fock_BZ%nbandocc_bz,fock%fock_BZ%mkpt,&
399 10 : & fock%fock_BZ%npwarr,dtset%nsppol,option,rprimd,ylmbz,ylmgrbz)
400 : call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,&
401 : & iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcgbz,&
402 : & fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,psps%mpsang,&
403 : & dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,&
404 : & dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,&
405 : & dtset%nsppol,fock%fock_common%my_nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,&
406 10 : & xred,ylmbz,ylmgrbz)
407 10 : ABI_FREE(ylmbz)
408 10 : ABI_FREE(ylmgrbz)
409 : else
410 : call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,&
411 : & iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcg,&
412 : & fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,mpi_enreg,psps%mpsang,&
413 : & dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,&
414 : & dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,&
415 : & dtset%nsppol,fock%fock_common%my_nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,&
416 16 : & xred,ylm,ylmgr)
417 : end if
418 : end if
419 : end if
420 6290 : usexg = 0
421 6290 : if (dtset%cprj_in_memory==1) usexg = 1
422 : call forstrnps(cg,cprj,dtset%ecut,dtset%ecutsm,dtset%effmass_free,eigen,electronpositron,fock,grnl,&
423 : & dtset%istwfk,kg,kinstr,nlstr,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,mggastr,dtset%mkmem,&
424 : & mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,nfftf,dtset%ngfft,&
425 : & dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,&
426 : & dtset%nucdipmom,occ,optfor,paw_ij,pawfgr,pawtab,ph1d,psps,dtset%qgbt,rprimd,stress_needed,symrec,dtset%typat,&
427 : & dtset%use_gbt,usecprj,dtset%usefock,usevxctau,vxctau,usexg,dtset%gpu_option,dtset%gpu_nl_distrib,&
428 6290 : & dtset%gpu_nl_splitsize,dtset%gpu_nfft_blocks,dtset%wtk,xred,ylm,ylmgr,xg_nonlop)
429 0 : else if (optfor>0) then !WVL
430 0 : ABI_MALLOC(xcart,(3, dtset%natom))
431 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
432 0 : call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart)
433 0 : ABI_FREE(xcart)
434 : end if
435 :
436 6333 : call timab(911,2,tsec)
437 6333 : call timab(912,1,tsec)
438 :
439 : !PAW: add gradients due to Dij derivatives to non-local term
440 6333 : if (psps%usepaw==1) then
441 5668 : ABI_MALLOC(vlocal,(nfftf,dtset%nspden))
442 :
443 : !$OMP PARALLEL DO COLLAPSE(2)
444 3086 : do ispden=1,min(dtset%nspden,2)
445 58227172 : do ifft=1,nfftf
446 58225755 : vlocal(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
447 : end do
448 : end do
449 1417 : if (dtset%nspden==4) then
450 : !$OMP PARALLEL DO COLLAPSE(2)
451 162 : do ispden=3,4
452 1319648 : do ifft=1,nfftf
453 1319594 : vlocal(ifft,ispden)=vxc(ifft,ispden)
454 : end do
455 : end do
456 : end if
457 1417 : ucvol_=ucvol
458 : #if defined HAVE_BIGDFT
459 : if (dtset%usewvl==1) ucvol_=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp)
460 : #endif
461 1417 : optgr=optfor;optgr2=0;optstr=stress_needed;optstr2=0
462 1417 : comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl
463 : call pawgrnl(atindx1,dtset%nspden,dummy,1,dummy,grnl,gsqcut,mgfftf,my_natom,dtset%natom,&
464 : & nattyp,nfftf,ngfftf,nhat,nlstr,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,&
465 : & pawang,pawfgrtab,pawrhoij,pawtab,ph1df,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred,&
466 1417 : & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=comm_grid)
467 : !DEBUG
468 : ! write(6,*)' after pawgrnl, nlstr=',nlstr(1:6)
469 : !ENDDEBUG
470 1417 : ABI_FREE(vlocal)
471 :
472 : end if
473 6333 : call timab(912,2,tsec)
474 6333 : call timab(913,1,tsec)
475 :
476 : !==========================================================================
477 : !Here compute forces (if required)
478 : !==========================================================================
479 6333 : if (optfor==1) then
480 : apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. &
481 3269 : & abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5)
482 :
483 : ! If residual is a density residual (and forces from residual asked),
484 : ! has to convert it into a potential residual before calling forces routine
485 : if (apply_residual) then
486 2852 : ABI_MALLOC(resid,(nfftf,dtset%nspden))
487 713 : option=0; if (dtset%densfor_pred<0) option=1
488 713 : optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2
489 : call nres2vres(dtset,gsqcut,psps%usepaw,kxc,mpi_enreg,my_natom,nfftf,ngfftf,nhat,&
490 : & nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,&
491 713 : & rhor,rprimd,psps%usepaw,resid,xccc3d,xred,vxc,xcctau3d)
492 : else
493 2556 : resid => nvresid
494 : end if
495 :
496 : call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,gred,grchempottn,grcondft,gresid,grewtn,&
497 : & grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfftf,&
498 : & mpi_enreg,psps%n1xccc,n3xccc,nattyp,&
499 : & nfftf,ngfftf,ngrvdw,ntypat,pawrad,pawtab,ph1df,psps,rhog,&
500 : & rhor,rprimd,symrec,synlgr,dtset%usefock,usevxctau,resid,vxc,vxctau,wvl%descr,wvl%den,xred,&
501 3269 : & electronpositron=electronpositron)
502 :
503 3269 : if (apply_residual) then
504 713 : ABI_FREE(resid)
505 : end if
506 : end if
507 :
508 6333 : call timab(913,2,tsec)
509 6333 : call timab(914,1,tsec)
510 :
511 : !==========================================================================
512 : !Here compute stress tensor (if required)
513 : !==========================================================================
514 :
515 6333 : if (stress_needed==1.and.dtset%usewvl==0) then
516 : ! if (dtset%usefock==1 .and. associated(fock).and.fock%fock_common%optstr.and.psps%usepaw==0) then
517 6302 : if (dtset%usefock==1 .and. associated(fock)) then
518 121 : if (fock%fock_common%optstr) then
519 484 : fock%fock_common%stress(1:3)=fock%fock_common%stress(1:3)-(two*energies%e_fock-energies%e_fock0)/ucvol
520 121 : if (n3xccc>0.and.psps%usepaw==0 .and. &
521 : & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) then
522 120 : ABI_MALLOC(vxc_hf,(nfftf,dtset%nspden))
523 : !compute Vxc^GGA(rho_val)
524 30 : call xchybrid_ncpp_cc(dtset,dum,dum2,mpi_enreg,nfftf,ngfftf,n3xccc,rhor,rprimd,strdum,dum1,xccc3d,vxc=vxc_hf,optstr=1)
525 : end if
526 : end if
527 : end if
528 : call stress(atindx1,dtset%berryopt,dtefield,energies%e_localpsp,dtset%efield,&
529 : & energies%e_hartree,energies%e_corepsp,fock,gsqcut,extfpmd,dtset%ixc,kinstr,mggastr,mgfftf,&
530 : & mpi_enreg,psps%mqgrid_vl,psps%n1xccc,n3xccc,dtset%natom,nattyp,&
531 : & nfftf,ngfftf,nlstr,dtset%nspden,dtset%nsym,ntypat,psps,pawrad,pawtab,ph1df,&
532 : & dtset%prtvol,psps%qgrid_vl,dtset%red_efieldbar,rhog,rprimd,strten,strscondft,strsxc,symrec,&
533 : & dtset%typat,dtset%usefock,dtset%usekden,psps%usepaw,usevxctau,&
534 : & dtset%vdw_tol,dtset%vdw_tol_3bt,dtset%vdw_xc,psps%vlspl,vxc,vxctau,vxc_hf,&
535 : & psps%xccc1d,xccc3d,xcctau3d,psps%xcccrc,xred,psps%ziontypat,psps%znucltypat,qvpotzero,&
536 6302 : & electronpositron=electronpositron)
537 : end if
538 :
539 : !Memory deallocation
540 6333 : ABI_FREE(grnl)
541 6333 : if (allocated(vxc_hf)) then
542 30 : ABI_FREE(vxc_hf)
543 : end if
544 :
545 :
546 6333 : call timab(914,2,tsec)
547 : ABI_NVTX_END_RANGE()
548 6333 : call timab(910,2,tsec)
549 :
550 12666 : end subroutine forstr
551 : !!***
552 :
553 : !----------------------------------------------------------------------
554 :
555 : !!****f* ABINIT/forstrnps
556 : !! NAME
557 : !! forstrnps
558 : !!
559 : !! FUNCTION
560 : !! Compute nonlocal pseudopotential energy contribution to forces and/or stress tensor
561 : !! as well as kinetic energy contribution to stress tensor.
562 : !!
563 : !! INPUTS
564 : !! cg(2,mcg)=wavefunctions (may be read from disk file)
565 : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
566 : !! ecut=cut-off energy for plane wave basis sphere (Ha)
567 : !! ecutsm=smearing energy for plane wave kinetic energy (Ha)
568 : !! effmass_free=effective mass for electrons (1. in common case)
569 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
570 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
571 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
572 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
573 : !! kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis
574 : !! kpt(3,nkpt)=k points in reduced coordinates
575 : !! mband=maximum number of bands
576 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
577 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
578 : !! mgfft=maximum size of 1D FFTs
579 : !! mkmem=number of k points treated by this node.
580 : !! mpi_enreg=information about MPI parallelization
581 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
582 : !! mpw= maximum number of plane waves
583 : !! my_natom=number of atoms treated by current processor
584 : !! natom=number of atoms in cell.
585 : !! nband(nkpt)=number of bands at each k point
586 : !! nfft=number of FFT grid points
587 : !! nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for NCPP)
588 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
589 : !! nkpt=number of k points in Brillouin zone
590 : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
591 : !! npwarr(nkpt)=number of planewaves in basis and boundary at each k
592 : !! nspden=Number of spin Density components
593 : !! nspinor=number of spinorial components of the wavefunctions
594 : !! nsppol=1 for unpolarized, 2 for spin-polarized
595 : !! nsym=number of elements in symmetry group
596 : !! ntypat=number of types of atoms
597 : !! nucdipmom(3,my_natom)= nuclear dipole moments
598 : !! occ(mband*nkpt*nsppol)=occupation numbers for each band over all k points
599 : !! optfor=1 if computation of forces is required
600 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
601 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
602 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
603 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
604 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
605 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
606 : !! stress_needed=1 if computation of stress tensor is required
607 : !! symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless)
608 : !! typat(natom)=type of each atom
609 : !! usecprj=1 if cprj datastructure has been allocated
610 : !! usefock=1, if Fock contribution has to be included
611 : !! usevxctau=1 if kinetic energy density contribution has to be included (mGGA)
612 : !! vxctau(nfftf,nspden,4*usevxctau)=(only for meta-GGA): derivative of XC energy density
613 : !! wrt kinetic energy density (depsxcdtau)
614 : !! gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
615 : !! wtk(nkpt)=weight associated with each k point
616 : !! xred(3,natom)=reduced dimensionless atomic coordinates
617 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
618 : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
619 : !!
620 : !! OUTPUT
621 : !! if (optfor==1)
622 : !! grnl(3*natom*optfor)=stores grads of nonlocal energy wrt atomic coordinates
623 : !! if (stress_needed==1)
624 : !! kinstr(6)=kinetic energy part of stress tensor (hartree/bohr^3)
625 : !! Store 6 unique components of symmetric 3x3 tensor in the order
626 : !! 11, 22, 33, 32, 31, 21
627 : !! mggastr(6)=meta-GGA part of stress tensor (hartree/bohr^3)
628 : !! Only non-local contribution from Div(V_tau.Grad(Psi))
629 : !! npsstr(6)=nonlocal pseudopotential energy part of stress tensor
630 : !! (hartree/bohr^3)
631 : !!
632 : !! SOURCE
633 :
634 6290 : subroutine forstrnps(cg,cprj,ecut,ecutsm,effmass_free,eigen,electronpositron,fock,&
635 6290 : & grnl,istwfk,kg,kinstr,npsstr,kpt,mband,mcg,mcprj,mgfft,mggastr,mkmem,mpi_enreg,mpsang,&
636 6290 : & mpw,my_natom,natom,nband,nfft,nfftf,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,nsym,&
637 6290 : & ntypat,nucdipmom,occ,optfor,paw_ij,pawfgr,pawtab,ph1d,psps,qgbt,rprimd,&
638 6290 : & stress_needed,symrec,typat,use_gbt,usecprj,usefock,usevxctau,vxctau,usexg,&
639 6290 : & gpu_option,gpu_nl_distrib,gpu_nl_splitsize,gpu_nfft_blocks,wtk,xred,ylm,ylmgr,xg_nonlop)
640 :
641 : !Arguments ------------------------------------
642 : !scalars
643 : integer,intent(in) :: mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,nfft,nfftf,nkpt
644 : integer,intent(in) :: nspden,nsppol,nspinor,nsym,ntypat,optfor,stress_needed
645 : integer,intent(in) :: use_gbt,usecprj,usefock,usevxctau,usexg,gpu_option
646 : integer,intent(in) :: gpu_nl_distrib,gpu_nl_splitsize,gpu_nfft_blocks
647 : real(dp),intent(in) :: ecut,ecutsm,effmass_free
648 : type(electronpositron_type),pointer :: electronpositron
649 : type(MPI_type),intent(inout) :: mpi_enreg
650 : type(pawfgr_type), intent(in) :: pawfgr
651 : type(pseudopotential_type),intent(in) :: psps
652 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
653 : !arrays
654 : integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
655 : integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt)
656 : integer,intent(in) :: symrec(3,3,nsym),typat(natom)
657 : real(dp),intent(in) :: cg(2,mcg)
658 : real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kpt(3,nkpt),nucdipmom(3,my_natom)
659 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom),qgbt(3)
660 : real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom)
661 : real(dp),intent(in),target :: vxctau(nfftf,nspden,4*usevxctau)
662 : real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
663 : real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm)
664 : real(dp),intent(out) :: grnl(3*natom*optfor),kinstr(6),mggastr(6),npsstr(6)
665 : type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj)
666 : type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
667 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
668 : type(fock_type),pointer, intent(inout) :: fock
669 : !Local variables-------------------------------
670 : !scalars
671 : integer,parameter :: tim_rwwf=7
672 : integer :: bandpp,bdtot_index,choice,cpopt,dimffnl,dimffnl_str,iband,iband_cprj,iband_last,ibg,icg,ider,ider_str
673 : integer :: idir,idir_str,ierr,ii,ikg,ikpt,ilm,ipositron,ipw,ishift,isppol,istwf_k,npw_k_
674 : integer :: mband_cprj,me_distrb,my_ikpt,my_nspinor,nband_k,nband_cprj_k,ndat,nkpg
675 : integer :: nnlout,npw_k,paw_opt,signs,spaceComm
676 : integer :: tim_nonlop,tim_nonlop_prep,usecprj_local,use_ACE_old
677 : integer :: blocksize,iblock,iblocksize,ibs,nblockbd,nblk_gemm_nonlop
678 : integer :: space,me_g0,ncols_cprj,me_g0_fft
679 : real(dp) :: ar,ar2,renorm_factor,dfsm,ecutsm_inv,fact_kin,fsm,htpisq
680 : real(dp) :: kin,kin_kphq,xx
681 6290 : type(gs_hamiltonian_type),target :: gs_hamk
682 : logical :: compute_gbound,usefock_loc
683 : character(len=500) :: msg
684 : type(fock_common_type),pointer :: fockcommon
685 : type(xgBlock_t) :: xgx0_tr
686 6290 : type(xgTransposer_t) :: xgTransposer
687 : !arrays
688 6290 : integer,allocatable,target :: kg_k(:,:)
689 6290 : integer,pointer :: gbound_k(:,:),kg_k_(:,:)
690 : real(dp) :: kpoint(3),kphq(3),nonlop_dum(1,1),rmet(3,3),tsec(2)
691 : real(dp) :: kgr(3),kgr_kphq(3),kgc(3),kgc_kphq(3)
692 : #if defined HAVE_GPU && defined HAVE_YAKL
693 : real(c_double), contiguous, pointer :: cwavef(:,:) => null()
694 : #else
695 6290 : real(dp),allocatable,target :: cwavef(:,:)
696 : #endif
697 6290 : real(dp),allocatable :: enlout(:),ffnl_sav(:,:,:,:),ffnl_str(:,:,:,:)
698 6290 : real(dp),allocatable :: ghc_dum(:,:),gprimd(:,:),kpg_k(:,:),kpg_kphq(:,:),kpg_k_sav(:,:)
699 6290 : real(dp),allocatable :: kstr1(:),kstr2(:),kstr3(:),kstr4(:),kstr5(:),kstr6(:)
700 6290 : real(dp),allocatable :: kstr1_kphq(:),kstr2_kphq(:),kstr3_kphq(:),kstr4_kphq(:),kstr5_kphq(:),kstr6_kphq(:)
701 6290 : real(dp),allocatable :: lambda(:),occblock(:),ph3d(:,:,:),ph3d_kphq(:,:,:),ph3d_sav(:,:,:)
702 6290 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
703 6290 : real(dp),allocatable :: weight(:),ylm_k(:,:),ylmgr_k(:,:,:)
704 6290 : real(dp),allocatable,target :: ffnl(:,:,:,:),ffnl_kphq(:,:,:,:)
705 6290 : real(dp),pointer :: vxctau_ptr(:,:,:)
706 : type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
707 6290 : type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
708 6290 : type(pawcprj_type),pointer :: cwaveprj_idat(:,:)
709 : type(xgBlock_t) :: xgx0,xgeigen,xgforces,xgstress
710 : type(xg_t) :: cprj_xgx0,cprj_work
711 6290 : real(dp),allocatable :: enlout_2d(:,:),enlout_2d_stress(:,:)
712 6290 : real(dp),allocatable :: cwavef_spin(:,:),enlout_spin(:)
713 6290 : real(dp),pointer :: cwavef_tr(:,:)
714 : !*************************************************************************
715 :
716 : ABI_NVTX_START_RANGE(NVTX_FORSTRNPS)
717 6290 : call timab(920,1,tsec) ; call timab(921,-1,tsec)
718 :
719 : !Init mpicomm and me
720 6290 : if(mpi_enreg%paral_kgb==1)then
721 920 : spaceComm=mpi_enreg%comm_kpt
722 920 : me_distrb=mpi_enreg%me_kpt
723 : else
724 : !* In case of HF calculation
725 5370 : if (mpi_enreg%paral_hf==1) then
726 0 : spaceComm=mpi_enreg%comm_kpt
727 0 : me_distrb=mpi_enreg%me_kpt
728 : else
729 5370 : spaceComm=mpi_enreg%comm_cell
730 5370 : me_distrb=mpi_enreg%me_cell
731 : end if
732 : end if
733 :
734 : !Some constants
735 : ipositron=abs(electronpositron_calctype(electronpositron))
736 6290 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
737 : !Smearing of plane wave kinetic energy
738 6290 : ecutsm_inv=zero;if( ecutsm>1.0d-20) ecutsm_inv=1/ecutsm
739 : !htpisq is (1/2) (2 Pi) **2:
740 6290 : htpisq=0.5_dp*(two_pi)**2
741 :
742 : !Check that fock is present if want to use fock option
743 6290 : compute_gbound=.false.
744 12458 : usefock_loc = (usefock==1 .and. associated(fock))
745 : !Arrays initializations
746 37394 : grnl(:)=zero
747 6290 : if (usefock_loc) then
748 122 : fockcommon =>fock%fock_common
749 122 : fockcommon%optfor=.false.
750 122 : fockcommon%optstr=.false.
751 122 : use_ACE_old=fockcommon%use_ACE
752 122 : fockcommon%use_ACE=0
753 : if (fockcommon%optfor) compute_gbound=.true.
754 : end if
755 6290 : if (stress_needed==1) then
756 6259 : kinstr(:)=zero;mggastr(:)=0;npsstr(:)=zero
757 6259 : if (usevxctau>0) compute_gbound=.true.
758 6259 : if (usefock_loc) then
759 121 : fockcommon%optstr=.TRUE.
760 847 : fockcommon%stress=zero
761 121 : compute_gbound=.true.
762 : end if
763 : end if
764 :
765 6290 : usecprj_local=usecprj
766 :
767 6290 : if ((usefock_loc).and.(psps%usepaw==1)) then
768 26 : usecprj_local=1
769 26 : if(optfor==1)then
770 10 : fockcommon%optfor=.true.
771 10 : if (.not.allocated(fockcommon%forces_ikpt)) then
772 40 : ABI_MALLOC(fockcommon%forces_ikpt,(3,natom,mband))
773 : end if
774 10 : if (.not.allocated(fockcommon%forces)) then
775 30 : ABI_MALLOC(fockcommon%forces,(3,natom))
776 : end if
777 90 : fockcommon%forces=zero
778 10 : compute_gbound=.true.
779 : end if
780 : end if
781 :
782 : !Initialize Hamiltonian (k-independent terms)
783 : call gs_hamk%init(psps,pawtab,nspinor,nsppol,nspden,natom,&
784 : & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj_local,&
785 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
786 : & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,fock=fock,&
787 6290 : & nucdipmom=nucdipmom,gpu_option=gpu_option,nfft_blocks=gpu_nfft_blocks)
788 251600 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
789 :
790 6290 : if (usevxctau>0) then
791 259 : ABI_MALLOC(vxctaulocal,(gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,gs_hamk%nvloc,4))
792 : end if
793 :
794 6290 : if (usexg==1) then
795 166 : call xg_nonlop_update_weight(xg_nonlop,gs_hamk%ucvol)
796 166 : if (xg_nonlop%paw) call xg_nonlop_make_Dij(xg_nonlop,paw_ij,nsppol,gs_hamk%atindx)
797 : end if
798 :
799 : !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted)
800 6290 : if (psps%usepaw==1.and.usecprj_local==1) then
801 106 : call pawcprj_reorder(cprj,gs_hamk%atindx)
802 : end if
803 :
804 : !Common data for "nonlop" routine
805 6290 : signs=1 ; idir=0 ; ishift=0 ; tim_nonlop=4 ; tim_nonlop_prep=12
806 6290 : choice=2*optfor;if (stress_needed==1) choice=10*choice+3
807 6290 : if (optfor==1.and.stress_needed==1) ishift=6
808 6290 : nnlout=max(1,6*stress_needed+3*natom*optfor)
809 6290 : if (psps%usepaw==0) then
810 4873 : paw_opt=0 ; cpopt=-1
811 : else
812 1417 : paw_opt=2 ; cpopt=-1+3*usecprj_local
813 : end if
814 :
815 6290 : call timab(921,2,tsec)
816 :
817 : !LOOP OVER SPINS
818 6290 : bdtot_index=0;ibg=0;icg=0
819 13187 : do isppol=1,nsppol
820 :
821 : ! Continue to initialize the Hamiltonian (PAW DIJ coefficients)
822 6897 : call gs_hamk%load_spin(isppol,with_nonlocal=.true.)
823 6897 : if (usefock_loc) fockcommon%isppol=isppol
824 :
825 : ! If any, set up local potential vtau on the coarse FFT mesh
826 6897 : if (usevxctau>0) then
827 39 : vxctau_ptr => vxctau ! This is to bepass annoying inout attribute
828 : call gspot_transgrid_and_pack(isppol,psps%usepaw,mpi_enreg%paral_kgb,nfft,ngfft,nfftf, &
829 39 : nspden,gs_hamk%nvloc,4,pawfgr,mpi_enreg,vxctau_ptr,vxctaulocal)
830 39 : call gs_hamk%load_spin(isppol,vxctaulocal=vxctaulocal)
831 : end if
832 :
833 6897 : if (usexg==1) then
834 170 : if (xg_nonlop%paw) call xg_nonlop_set_Dij_spin(xg_nonlop,isppol)
835 : end if
836 :
837 : ! Loop over k points
838 6897 : ikg=0
839 73457 : do ikpt=1,nkpt
840 60270 : if (usefock_loc) fockcommon%ikpt=ikpt
841 60270 : nband_k=nband(ikpt+(isppol-1)*nkpt)
842 60270 : istwf_k=istwfk(ikpt)
843 60270 : npw_k=npwarr(ikpt)
844 241080 : kpoint(:)=kpt(:,ikpt)
845 60270 : if (use_gbt/=0) then
846 256 : kpoint(:)=kpt(:,ikpt)-half*qgbt(:)
847 256 : kphq(:)=kpt(:,ikpt)+half*qgbt(:)
848 : end if
849 :
850 60270 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
851 15868 : bdtot_index=bdtot_index+nband_k
852 15868 : cycle
853 : end if
854 :
855 44402 : call timab(922,1,tsec)
856 :
857 44402 : my_ikpt=mpi_enreg%my_kpttab(ikpt)
858 : ! Parallelism over FFT and/or bands: define sizes and tabs
859 44402 : if (mpi_enreg%paral_kgb==1) then
860 3764 : nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
861 3764 : bandpp=mpi_enreg%bandpp
862 3764 : my_bandfft_kpt => bandfft_kpt(my_ikpt)
863 : else
864 40638 : bandpp=mpi_enreg%bandpp
865 40638 : nblockbd=nband_k/bandpp
866 : end if
867 44402 : blocksize=nband_k/nblockbd
868 44402 : mband_cprj=mband/mpi_enreg%nproc_band
869 44402 : nband_cprj_k=nband_k/mpi_enreg%nproc_band
870 :
871 44402 : if(gpu_option == ABI_GPU_KOKKOS) then
872 : #if defined HAVE_GPU && defined HAVE_YAKL
873 : ABI_MALLOC_MANAGED(cwavef,(/2,npw_k*my_nspinor*blocksize/))
874 : #endif
875 : else
876 133206 : ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
877 : end if
878 :
879 44402 : if (psps%usepaw==1.and.usecprj_local==1) then
880 10589 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp))
881 1517 : call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
882 : else
883 42885 : ABI_MALLOC(cwaveprj,(0,0))
884 : end if
885 :
886 44402 : if (stress_needed==1) then
887 132519 : ABI_MALLOC(kstr1,(npw_k))
888 88346 : ABI_MALLOC(kstr2,(npw_k))
889 88346 : ABI_MALLOC(kstr3,(npw_k))
890 88346 : ABI_MALLOC(kstr4,(npw_k))
891 88346 : ABI_MALLOC(kstr5,(npw_k))
892 88346 : ABI_MALLOC(kstr6,(npw_k))
893 44173 : if (use_gbt /= 0) then
894 128 : ABI_MALLOC(kstr1_kphq,(npw_k))
895 128 : ABI_MALLOC(kstr2_kphq,(npw_k))
896 128 : ABI_MALLOC(kstr3_kphq,(npw_k))
897 128 : ABI_MALLOC(kstr4_kphq,(npw_k))
898 128 : ABI_MALLOC(kstr5_kphq,(npw_k))
899 128 : ABI_MALLOC(kstr6_kphq,(npw_k))
900 : end if
901 : end if
902 :
903 133206 : ABI_MALLOC(kg_k,(3,mpw))
904 : !$OMP PARALLEL DO
905 10071837 : do ipw=1,npw_k
906 40154142 : kg_k(:,ipw)=kg(:,ipw+ikg)
907 : end do
908 :
909 177608 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
910 44402 : if (stress_needed==1) then
911 176692 : ABI_MALLOC(ylmgr_k,(npw_k,3,mpsang*mpsang*psps%useylm))
912 : else
913 229 : ABI_MALLOC(ylmgr_k,(0,0,0))
914 : end if
915 44402 : if (psps%useylm==1) then
916 : !$OMP PARALLEL DO COLLAPSE(2)
917 61903 : do ilm=1,mpsang*mpsang
918 12466806 : do ipw=1,npw_k
919 12457453 : ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
920 : end do
921 : end do
922 9353 : if (stress_needed==1) then
923 : !$OMP PARALLEL DO COLLAPSE(2)
924 60608 : do ilm=1,mpsang*mpsang
925 215006 : do ii=1,3
926 37252156 : do ipw=1,npw_k
927 37200690 : ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm)
928 : end do
929 : end do
930 : end do
931 : end if
932 : end if
933 :
934 :
935 : ! Prepare kinetic contribution to stress tensor (Warning : the symmetry
936 : ! has not been broken, like in mkkin.f or kpg3.f . It should be, in order to be coherent).
937 44402 : if (stress_needed==1) then
938 44173 : ABI_MALLOC(gprimd,(3,3))
939 618422 : gprimd=gs_hamk%gprimd
940 : !$OMP PARALLEL DO PRIVATE(fact_kin,ipw,kgr,kgr_kphq,kgc,kgc_kphq,kin,xx,fsm,dfsm) &
941 : !$OMP&SHARED(ecut,ecutsm,ecutsm_inv,gs_hamk,htpisq,kg_k,kpoint,kphq,kstr1,kstr2,kstr3,kstr4,kstr5,kstr6,kstr1_kphq,kstr2_kphq,kstr3_kphq,kstr4_kphq,kstr5_kphq,kstr6_kphq,npw_k)
942 10050876 : do ipw=1,npw_k
943 : ! Compute Cartesian coordinates of (k+G)
944 40026812 : kgr = kpoint + kg_k(:,ipw)
945 130087139 : kgc = MATMUL(gprimd, kgr)
946 40026812 : kin = htpisq * DOT_PRODUCT(kgc, kgc)
947 : ! kgc1=gprimd(1,1)*(kpoint(1)+kg_k(1,ipw))+&
948 : ! & gprimd(1,2)*(kpoint(2)+kg_k(2,ipw))+&
949 : ! & gprimd(1,3)*(kpoint(3)+kg_k(3,ipw))
950 : ! kgc2=gprimd(2,1)*(kpoint(1)+kg_k(1,ipw))+&
951 : ! & gprimd(2,2)*(kpoint(2)+kg_k(2,ipw))+&
952 : ! & gprimd(2,3)*(kpoint(3)+kg_k(3,ipw))
953 : ! kgc3=gprimd(3,1)*(kpoint(1)+kg_k(1,ipw))+&
954 : ! & gprimd(3,2)*(kpoint(2)+kg_k(2,ipw))+&
955 : ! & gprimd(3,3)*(kpoint(3)+kg_k(3,ipw))
956 : ! kin=htpisq* ( kgc1**2 + kgc2**2 + kgc3**2 )
957 10006703 : fact_kin=1.0_dp
958 10006703 : if (kin>ecut-ecutsm) then
959 494750 : if (kin>ecut) then
960 : fact_kin=0.0_dp
961 : else
962 : ! See the routine mkkin.f, for the smearing procedure
963 204208 : xx=(ecut-kin)*ecutsm_inv
964 : ! This kinetic cutoff smoothing function and its xx derivatives
965 : ! were produced with Mathematica and the fortran code has been
966 : ! numerically checked against Mathematica.
967 204208 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
968 204208 : dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
969 : ! d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
970 : ! & (-144+45*xx))))))*fsm**3
971 204208 : fact_kin=fsm+kin*(-ecutsm_inv)*dfsm
972 : end if
973 : end if
974 10006703 : kstr1(ipw) = fact_kin * kgc(1) * kgc(1)
975 10006703 : kstr2(ipw) = fact_kin * kgc(2) * kgc(2)
976 10006703 : kstr3(ipw) = fact_kin * kgc(3) * kgc(3)
977 10006703 : kstr4(ipw) = fact_kin * kgc(3) * kgc(2)
978 10006703 : kstr5(ipw) = fact_kin * kgc(3) * kgc(1)
979 10050876 : kstr6(ipw) = fact_kin * kgc(2) * kgc(1)
980 : ! kstr1(ipw)=fact_kin*kgc1*kgc1
981 : ! kstr2(ipw)=fact_kin*kgc2*kgc2
982 : ! kstr3(ipw)=fact_kin*kgc3*kgc3
983 : ! kstr4(ipw)=fact_kin*kgc3*kgc2
984 : ! kstr5(ipw)=fact_kin*kgc3*kgc1
985 : ! kstr6(ipw)=fact_kin*kgc2*kgc1
986 : end do ! ipw
987 44173 : if (use_gbt/=0) then
988 19358 : do ipw=1,npw_k
989 77176 : kgr_kphq = kphq+kg_k(:,ipw)
990 250822 : kgc_kphq = MATMUL(gprimd, kgr_kphq)
991 77176 : kin_kphq = htpisq * DOT_PRODUCT(kgc_kphq, kgc_kphq)
992 19294 : fact_kin=1.0_dp
993 19294 : if (kin_kphq>ecut-ecutsm) then
994 679 : if (kin_kphq>ecut) then
995 : fact_kin=0.0_dp
996 : else
997 291 : xx=(ecut-kin_kphq)*ecutsm_inv
998 291 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
999 291 : dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
1000 291 : fact_kin=fsm+kin_kphq*(-ecutsm_inv)*dfsm
1001 : end if
1002 : end if
1003 19294 : kstr1_kphq(ipw)=fact_kin*kgc_kphq(1)*kgc_kphq(1)
1004 19294 : kstr2_kphq(ipw)=fact_kin*kgc_kphq(2)*kgc_kphq(2)
1005 19294 : kstr3_kphq(ipw)=fact_kin*kgc_kphq(3)*kgc_kphq(3)
1006 19294 : kstr4_kphq(ipw)=fact_kin*kgc_kphq(3)*kgc_kphq(2)
1007 19294 : kstr5_kphq(ipw)=fact_kin*kgc_kphq(3)*kgc_kphq(1)
1008 19358 : kstr6_kphq(ipw)=fact_kin*kgc_kphq(2)*kgc_kphq(1)
1009 : end do ! ipw
1010 : end if ! GBT
1011 44173 : ABI_FREE(gprimd)
1012 : end if
1013 :
1014 :
1015 : ! Compute (k+G) vectors (only if useylm=1)
1016 44402 : if (usexg/=1) then
1017 43242 : nkpg=3*nloalg(3)
1018 172968 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
1019 43242 : if (nkpg>0) call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1020 43242 : if (use_gbt/=0) then
1021 192 : ABI_MALLOC(kpg_kphq,(npw_k,nkpg))
1022 64 : call mkkpg(kg_k,kpg_kphq,kphq,nkpg,npw_k)
1023 : end if
1024 : else ! cprj_in_memory = 1
1025 1160 : nkpg=3
1026 4640 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
1027 1160 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1028 1160 : if (use_gbt/=0) then
1029 0 : ABI_MALLOC(kpg_kphq,(npw_k,nkpg))
1030 0 : call mkkpg(kg_k,kpg_kphq,kphq,nkpg,npw_k)
1031 : end if
1032 : end if
1033 :
1034 : ! Compute nonlocal form factors ffnl at all (k+G)
1035 44402 : ider=0;idir=0;dimffnl=1
1036 44402 : if (stress_needed==1) then
1037 44173 : ider=1;dimffnl=2+2*psps%useylm
1038 : end if
1039 266412 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
1040 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1041 : & ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1042 44402 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
1043 44402 : if (use_gbt/=0) then
1044 320 : ABI_MALLOC(ffnl_kphq,(npw_k,dimffnl,psps%lmnmax,ntypat)) ! 0
1045 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_kphq,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1046 : & ider,idir,psps%indlmn,kg_k,kpg_kphq,kphq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1047 64 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
1048 : end if
1049 : #ifdef HAVE_OPENMP_OFFLOAD
1050 : !$OMP TARGET ENTER DATA MAP(to:ffnl) IF(gpu_option==ABI_GPU_OPENMP)
1051 : #endif
1052 44402 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
1053 54 : ider_str=1; dimffnl_str=7;idir_str=-7
1054 324 : ABI_MALLOC(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,ntypat))
1055 : call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1056 : & ider_str,idir_str,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1057 54 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
1058 : #ifdef HAVE_OPENMP_OFFLOAD
1059 : !$OMP TARGET ENTER DATA MAP(to:ffnl_str) IF(gpu_option==ABI_GPU_OPENMP)
1060 : #endif
1061 : end if
1062 :
1063 : ! Load k-dependent part in the Hamiltonian datastructure
1064 : ! - Compute 3D phase factors
1065 : ! - Prepare various tabs in case of band-FFT parallelism
1066 : ! - Load k-dependent quantities in the Hamiltonian
1067 177608 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
1068 : call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
1069 44402 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.)
1070 44402 : if (use_gbt/=0) then
1071 256 : ABI_MALLOC(ph3d_kphq,(2,npw_k,gs_hamk%matblk))
1072 : call gs_hamk%load_kprime(kpt_kp=kphq,&
1073 64 : & kpg_kp=kpg_kphq,ffnl_kp=ffnl_kphq,ph3d_kp=ph3d_kphq,compute_gbound=compute_gbound,compute_ph3d=.true.)
1074 : end if
1075 :
1076 : ! Load band-FFT tabs (transposed k-dependent arrays)
1077 44402 : if (mpi_enreg%paral_kgb==1.and.usexg/=1) then
1078 2950 : call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1079 2950 : call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg)
1080 : call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
1081 : & gbound_k =my_bandfft_kpt%gbound, &
1082 : & kg_k =my_bandfft_kpt%kg_k_gather, &
1083 : & kpg_k =my_bandfft_kpt%kpg_k_gather, &
1084 : & ffnl_k =my_bandfft_kpt%ffnl_gather, &
1085 2950 : & ph3d_k =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound)
1086 : end if
1087 :
1088 : ! If OpenMP GPU, load "hamiltonian" on GPU device
1089 44402 : if (gpu_option == ABI_GPU_OPENMP) then
1090 0 : if(mpi_enreg%paral_kgb==0) then
1091 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k)
1092 0 : else if(gs_hamk%istwf_k==1) then
1093 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k,kg_k_gather=my_bandfft_kpt%kg_k_gather)
1094 : else
1095 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k,kg_k_gather=my_bandfft_kpt%kg_k_gather_sym)
1096 : end if
1097 : end if
1098 :
1099 : ! Setup gemm_nonlop
1100 44402 : if (gemm_nonlop_use_gemm) then
1101 : call set_gemm_nonlop_ikpt(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%istwf_k,gs_hamk%indlmn,&
1102 108 : & gs_hamk%ntypat,gs_hamk%nattyp,gs_hamk%gpu_option)
1103 108 : call reset_gemm_nonlop()
1104 :
1105 108 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1106 0 : gemm_nonlop_block_size = gpu_nl_splitsize
1107 : call get_gemm_nonlop_ompgpu_blocksize(my_ikpt,gs_hamk,mpi_enreg%bandpp,nband_k,&
1108 : & nspinor,1,mpi_enreg%paral_kgb,mpi_enreg%nproc_band,&
1109 : & optfor,stress_needed,-1,gs_hamk%gpu_option,(gpu_nl_distrib/=0),&
1110 : & gemm_nonlop_block_size,nblk_gemm_nonlop,gs_hamk%nfft_blocks,&
1111 0 : & warn_on_fail=.true.,disable_output=(ikpt>1))
1112 0 : gemm_nonlop_is_distributed = (gpu_nl_distrib/=0 .and. nblk_gemm_nonlop > 0)
1113 0 : if(nblk_gemm_nonlop==-1) then
1114 0 : gs_hamk%gpu_option=ABI_GPU_DISABLED
1115 0 : ABI_WARNING("GPU has been disabled for forces and stress computation due to memory constraints.")
1116 : end if
1117 : end if
1118 : end if
1119 :
1120 44402 : if (usexg==1) then
1121 1160 : call xg_nonlop_make_k(xg_nonlop,my_ikpt,istwf_k,mpi_enreg%me_g0,mpi_enreg%me_g0_fft,npw_k,ffnl,ph3d,kpg_k,.true.)
1122 : end if
1123 :
1124 : ! Loop over (blocks of) bands; accumulate forces and/or stresses
1125 : ! The following is now wrong. In sequential, nblockbd=nband_k/bandpp
1126 : ! blocksize= bandpp (JB 2016/04/16)
1127 : ! Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
1128 133206 : ABI_MALLOC(lambda,(blocksize))
1129 88804 : ABI_MALLOC(occblock,(blocksize))
1130 88804 : ABI_MALLOC(weight,(blocksize))
1131 133206 : ABI_MALLOC(enlout,(nnlout*blocksize))
1132 44402 : if (usexg==1) then
1133 4640 : ABI_MALLOC(enlout_2d,(3*natom,blocksize*optfor))
1134 3480 : ABI_MALLOC(enlout_2d_stress,(6,blocksize*stress_needed))
1135 48644 : if (optfor==1) enlout_2d=zero
1136 67604 : if (stress_needed==1) enlout_2d_stress=zero
1137 : end if
1138 1654764 : occblock=zero;weight=zero;enlout(:)=zero
1139 44402 : if (usefock_loc) then
1140 536 : if (fockcommon%optstr) then
1141 1599 : ABI_MALLOC(fockcommon%stress_ikpt,(6,nband_k))
1142 19041 : fockcommon%stress_ikpt=zero
1143 : end if
1144 : end if
1145 536 : if ((usefock_loc).and.(psps%usepaw==1)) then
1146 57 : if (fockcommon%optfor) then
1147 1528 : fockcommon%forces_ikpt=zero
1148 : end if
1149 : end if
1150 :
1151 44402 : call timab(922,2,tsec)
1152 :
1153 44402 : if (usexg==1) then
1154 1160 : ncols_cprj = bandpp*my_nspinor
1155 1160 : call xg_init(cprj_xgx0,xg_nonlop%space_cprj,xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
1156 1160 : call xg_init(cprj_work,xg_nonlop%space_cprj,xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
1157 : end if
1158 :
1159 275026 : do iblock=1,nblockbd
1160 :
1161 230624 : iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
1162 230624 : iband_cprj=(iblock-1)*bandpp+1
1163 230624 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
1164 :
1165 : ! Select occupied bandsddk
1166 506182 : occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
1167 736806 : if( abs(maxval(occblock))>=tol8 ) then
1168 156728 : call timab(923,1,tsec)
1169 350342 : weight(:)=wtk(ikpt)*occblock(:)
1170 :
1171 : ! gs_hamk%ffnl_k is changed in fock_getghc, so that it is necessary to restore it when stresses are to be calculated.
1172 156728 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
1173 202 : call gs_hamk%load_k(ffnl_k=ffnl)
1174 : end if
1175 :
1176 : ! Load contribution from n,k
1177 : cwavef(:,1:npw_k*my_nspinor*blocksize)=&
1178 178332179 : & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
1179 156728 : if (psps%usepaw==1.and.usecprj_local==1) then
1180 : call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,&
1181 : & mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,&
1182 8441 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1183 : end if
1184 :
1185 156728 : call timab(923,2,tsec)
1186 156728 : call timab(924,-1,tsec)
1187 :
1188 350342 : lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
1189 : ABI_NVTX_START_RANGE(NVTX_FORSTR_NONLOP)
1190 156728 : if (mpi_enreg%paral_kgb/=1.and.usexg/=1) then
1191 147846 : if (use_gbt == 0) then
1192 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
1193 146391 : & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
1194 : else
1195 1455 : gs_hamk%nspinor=1
1196 4365 : ABI_MALLOC(cwavef_spin, (2, npw_k*blocksize))
1197 4365 : ABI_MALLOC(enlout_spin,(nnlout*blocksize))
1198 17088 : enlout_spin(:) = zero
1199 1455 : gs_hamk%ispin_gbt = 1
1200 1455 : call cg_copy_spin(1,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1201 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
1202 1455 : & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef_spin,cwavef_spin,select_k=K_H_K)
1203 :
1204 1455 : gs_hamk%ispin_gbt = 2
1205 1455 : call cg_copy_spin(2,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1206 : call nonlop(choice,cpopt,cwaveprj,enlout_spin,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
1207 1455 : & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef_spin,cwavef_spin,select_k=KPRIME_H_KPRIME)
1208 :
1209 17088 : enlout(1:nnlout*blocksize) = enlout(1:nnlout*blocksize) + enlout_spin(1:nnlout*blocksize)
1210 1455 : gs_hamk%nspinor = 2
1211 1455 : ABI_FREE(cwavef_spin)
1212 1455 : ABI_FREE(enlout_spin)
1213 : end if ! GBT
1214 147846 : if((stress_needed==1).and.(usevxctau==1)) then
1215 66 : cwavef_tr=>cwavef
1216 : endif
1217 8882 : else if (usexg/=1) then ! paral_kgb = 1
1218 : ! here we MUST pass option gpu_option=ABI_GPU_DISABLED, as cwavef here is a host memory buffer
1219 7534 : if((stress_needed==1).and.(usevxctau==1).and.mpi_enreg%nproc_band>1) then
1220 192 : ABI_MALLOC(cwavef_tr,(2,my_bandfft_kpt%ndatarecv*my_nspinor*mpi_enreg%bandpp))
1221 : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,&
1222 : & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,cwavef,&
1223 64 : & already_transposed=.False.,gpu_option=ABI_GPU_DISABLED,cwavef_tr=cwavef_tr)
1224 : else
1225 : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,&
1226 : & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,cwavef,&
1227 7470 : & already_transposed=.False.,gpu_option=ABI_GPU_DISABLED)
1228 7470 : if((stress_needed==1).and.(usevxctau==1)) cwavef_tr=>cwavef
1229 : endif
1230 : else ! usexg==1
1231 1348 : if ( istwf_k > 1 ) then ! Real only
1232 529 : space = SPACE_CR
1233 : else ! complex
1234 819 : space = SPACE_C
1235 : end if
1236 1348 : me_g0 = -1
1237 1348 : me_g0_fft = -1
1238 1348 : if (space==SPACE_CR) then
1239 529 : me_g0 = 0
1240 529 : me_g0_fft = 0
1241 529 : if (istwf_k == 2) then
1242 81 : if (mpi_enreg%me_g0 == 1) me_g0 = 1
1243 81 : if (mpi_enreg%me_g0_fft == 1) me_g0_fft = 1
1244 : end if
1245 : end if
1246 :
1247 : call xgBlock_map(xgx0,cwavef,space,npw_k*my_nspinor,blocksize,comm=mpi_enreg%comm_band,me_g0=me_g0,&
1248 1348 : & gpu_option=gpu_option)
1249 :
1250 1348 : if ((stress_needed==1).and.(usevxctau==1)) then
1251 66 : if(mpi_enreg%nproc_band>1) then
1252 : call xgTransposer_constructor(xgTransposer,xgx0,xgx0_tr,my_nspinor,&
1253 64 : STATE_LINALG,TRANS_ALL2ALL,mpi_enreg%comm_spinorfft,mpi_enreg%comm_band,0,0,me_g0_fft)
1254 64 : call xgTransposer_transpose(xgTransposer,STATE_COLSROWS)
1255 64 : call xgBlock_reverseMap(xgx0_tr,cwavef_tr,rows=1,cols=my_nspinor*mpi_enreg%bandpp*my_bandfft_kpt%ndatarecv) !,rows=1,cols=spacedim*blockdim)
1256 : else
1257 2 : cwavef_tr=>cwavef
1258 : endif
1259 : endif
1260 :
1261 :
1262 1348 : call xgBlock_map_1d(xgeigen,lambda,SPACE_R,blocksize)
1263 :
1264 1348 : if (psps%usepaw==1.and.usecprj_local==1) then
1265 0 : call xg_cprj_copy(cwaveprj,cprj_xgx0%self,xg_nonlop,CPRJ_TO_XG)
1266 : else
1267 1348 : call xg_nonlop_getcprj(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self)
1268 : end if
1269 :
1270 1348 : if (optfor==1) call xgBlock_map(xgforces,enlout_2d,SPACE_R,3*natom,blocksize)
1271 1348 : if (stress_needed==1) call xgBlock_map(xgstress,enlout_2d_stress,SPACE_R,6,blocksize)
1272 :
1273 1348 : if (optfor==1.and.stress_needed==0) then
1274 : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,&
1275 208 : forces=xgforces)
1276 : end if
1277 1348 : if (optfor==0.and.stress_needed==1) then
1278 : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,&
1279 800 : stress=xgstress,gprimd=gs_hamk%gprimd)
1280 : end if
1281 1348 : if (optfor==1.and.stress_needed==1) then
1282 : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,&
1283 340 : forces=xgforces,stress=xgstress,gprimd=gs_hamk%gprimd)
1284 : end if
1285 :
1286 : end if ! end usexg==1
1287 : ABI_NVTX_END_RANGE()
1288 :
1289 156728 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
1290 202 : call gs_hamk%load_k(ffnl_k=ffnl_str)
1291 : end if
1292 :
1293 : ! Accumulate non-local contributions from n,k
1294 156728 : if (usexg/=1) then
1295 155380 : if (optfor==1) then
1296 240745 : do iblocksize=1,blocksize
1297 128633 : ibs=nnlout*(iblocksize-1)
1298 2094418 : grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout(ibs+1+ishift:ibs+3*natom+ishift)
1299 : end do
1300 : end if
1301 155380 : if (stress_needed==1) then
1302 335902 : do iblocksize=1,blocksize
1303 180626 : ibs=nnlout*(iblocksize-1)
1304 1419658 : npsstr(1:6)=npsstr(1:6) + weight(iblocksize)*enlout(ibs+1:ibs+6)
1305 : end do
1306 : end if
1307 : else
1308 1348 : if (optfor==1) then
1309 7128 : do iblocksize=1,blocksize
1310 49488 : grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout_2d(1:3*natom,iblocksize)
1311 : end do
1312 : end if
1313 1348 : if (stress_needed==1) then
1314 11208 : do iblocksize=1,blocksize
1315 71616 : npsstr(1:6) = npsstr(1:6) + weight(iblocksize)*enlout_2d_stress(1:6,iblocksize)
1316 : end do
1317 : end if
1318 : end if
1319 :
1320 156728 : call timab(924,2,tsec)
1321 :
1322 : #if defined HAVE_GPU && defined HAVE_YAKL
1323 : if(gpu_option==ABI_GPU_KOKKOS) then
1324 : ! the following is done on CPU, so prefetch wave functions from device to host (for efficiency)
1325 : call gpu_data_prefetch_async(C_LOC(cwavef), INT(2, c_size_t)*npw_k*my_nspinor*blocksize, CPU_DEVICE_ID)
1326 : call gpu_device_synchronize()
1327 : end if
1328 : #endif
1329 :
1330 : ! Accumulate stress tensor kinetic contributions
1331 156728 : if (stress_needed==1) then
1332 156416 : call timab(925,1,tsec)
1333 156416 : if (use_gbt == 0) then
1334 344200 : do iblocksize=1,blocksize
1335 : call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1336 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1337 189239 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1338 189239 : kinstr(1)=kinstr(1)+weight(iblocksize)*ar
1339 : call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1340 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1341 189239 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1342 189239 : kinstr(2)=kinstr(2)+weight(iblocksize)*ar
1343 : call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1344 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1345 189239 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1346 189239 : kinstr(3)=kinstr(3)+weight(iblocksize)*ar
1347 : call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1348 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1349 189239 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1350 189239 : kinstr(4)=kinstr(4)+weight(iblocksize)*ar
1351 : call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1352 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1353 189239 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1354 189239 : kinstr(5)=kinstr(5)+weight(iblocksize)*ar
1355 : call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1356 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1357 189239 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1358 344200 : kinstr(6)=kinstr(6)+weight(iblocksize)*ar
1359 : end do
1360 : else
1361 2910 : do iblocksize=1,blocksize
1362 : call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,1,&
1363 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1364 1455 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1365 : call meanvalue_g(ar2,kstr1_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1366 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1367 1455 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1368 1455 : kinstr(1)=kinstr(1)+weight(iblocksize)*(ar+ar2) ! (1,1)
1369 : call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,1,&
1370 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1371 1455 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1372 : call meanvalue_g(ar2,kstr2_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1373 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1374 1455 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1375 1455 : kinstr(2)=kinstr(2)+weight(iblocksize)*(ar+ar2) ! (2,2)
1376 : call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,1,&
1377 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1378 1455 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1379 : call meanvalue_g(ar2,kstr3_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1380 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1381 1455 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1382 1455 : kinstr(3)=kinstr(3)+weight(iblocksize)*(ar+ar2) ! (3,3)
1383 : call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,1,&
1384 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1385 1455 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1386 : call meanvalue_g(ar2,kstr4_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1387 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1388 1455 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1389 1455 : kinstr(4)=kinstr(4)+weight(iblocksize)*(ar+ar2) ! (3,2)
1390 : call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,1,&
1391 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1392 1455 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1393 : call meanvalue_g(ar2,kstr5_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1394 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1395 1455 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1396 1455 : kinstr(5)=kinstr(5)+weight(iblocksize)*(ar+ar2) ! (3,1)
1397 : call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,1,&
1398 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1399 1455 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1400 : call meanvalue_g(ar2,kstr6_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1401 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1402 1455 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1403 2910 : kinstr(6)=kinstr(6)+weight(iblocksize)*(ar+ar2) ! (2,1)`
1404 : end do
1405 : end if ! GBT
1406 156416 : call timab(925,2,tsec)
1407 : end if
1408 :
1409 : ! Accumulate stress tensor in case meta-GGA using v_tau
1410 156728 : if ((stress_needed==1).and.(usevxctau==1)) then
1411 197 : if(mpi_enreg%paral_kgb==1) then
1412 131 : gbound_k => my_bandfft_kpt%gbound
1413 131 : kg_k_ => my_bandfft_kpt%kg_k_gather
1414 131 : npw_k_=my_bandfft_kpt%ndatarecv
1415 : else
1416 66 : gbound_k => gs_hamk%gbound_k
1417 66 : kg_k_ => kg_k
1418 66 : npw_k_=npw_k
1419 : endif
1420 : call stress_mGGA(mggastr,cwavef_tr,effmass_free,gbound_k,gs_hamk%gprimd,istwf_k, &
1421 : & kg_k_,kpoint,mgfft,mpi_enreg,my_nspinor,mpi_enreg%bandpp,ngfft,npw_k_,gs_hamk%nvloc, &
1422 : & gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occblock(1+mpi_enreg%me_band*mpi_enreg%bandpp:(mpi_enreg%me_band+1)*mpi_enreg%bandpp),gs_hamk%ucvol,vxctaulocal, &
1423 197 : & wtk(ikpt),gpu_option=gpu_option)
1424 197 : if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>1) then
1425 128 : if(usexg==1) then
1426 64 : call xgTransposer_free(xgTransposer)
1427 : else
1428 64 : ABI_FREE(cwavef_tr)
1429 : endif
1430 : else
1431 69 : nullify(cwavef_tr)
1432 : endif
1433 : end if
1434 :
1435 : ! Accumulate stress tensor and forces for the Fock part
1436 156728 : if (usefock_loc) then
1437 2052 : if(fockcommon%optstr.or.fockcommon%optfor) then
1438 2052 : call timab(926,1,tsec)
1439 2052 : if (mpi_enreg%paral_kgb==1) then
1440 0 : msg='forsrtnps: Paral_kgb is not yet implemented for fock stresses'
1441 0 : ABI_BUG(msg)
1442 : end if
1443 2052 : ndat=mpi_enreg%bandpp
1444 2052 : if (gs_hamk%usepaw==0) cwaveprj_idat => cwaveprj
1445 2052 : ABI_MALLOC(ghc_dum,(0,0))
1446 2052 : fockcommon%ieigen=(iblock-1)*blocksize+1
1447 2052 : fockcommon%iband=(iblock-1)*blocksize+1
1448 2052 : if (gs_hamk%usepaw==1) then
1449 214 : cwaveprj_idat => cwaveprj(:,1+blocksize*(my_nspinor-1):blocksize*my_nspinor)
1450 : end if
1451 : call fock_getghc(cwavef(:,1+blocksize*npw_k*(my_nspinor-1):blocksize*npw_k*my_nspinor),cwaveprj_idat,&
1452 2052 : & ghc_dum,gs_hamk,mpi_enreg,blocksize)
1453 4104 : do iblocksize=1,blocksize
1454 2052 : if (fockcommon%optstr) then
1455 14280 : fockcommon%stress(:)=fockcommon%stress(:)+weight(iblocksize)*fockcommon%stress_ikpt(:,fockcommon%ieigen+iblocksize-1)
1456 : end if
1457 4104 : if (fockcommon%optfor) then
1458 1224 : fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen+iblocksize-1)
1459 : end if
1460 : end do
1461 2052 : ABI_FREE(ghc_dum)
1462 2052 : call timab(926,2,tsec)
1463 : end if
1464 : end if ! usefock_loc
1465 : end if
1466 275026 : if ( gpu_option == ABI_GPU_OPENMP) then
1467 0 : call ompgpu_free_hamilt_buffers()
1468 : end if
1469 :
1470 : end do ! End of loop on block of bands
1471 :
1472 44402 : if (usexg==1) then
1473 1160 : call xg_free(cprj_xgx0)
1474 1160 : call xg_free(cprj_work)
1475 : end if
1476 :
1477 44402 : call timab(927,1,tsec)
1478 :
1479 : ! Restore the bandfft tabs
1480 44402 : if (mpi_enreg%paral_kgb==1) then
1481 3764 : call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1482 : end if
1483 :
1484 : ! Increment indexes
1485 44402 : bdtot_index=bdtot_index+nband_k
1486 44402 : if (mkmem/=0) then
1487 44402 : ibg=ibg+my_nspinor*nband_cprj_k
1488 44402 : icg=icg+npw_k*my_nspinor*nband_k
1489 44402 : ikg=ikg+npw_k
1490 : end if
1491 :
1492 44402 : if (usefock_loc) then
1493 536 : if (fockcommon%optstr) then
1494 533 : ABI_FREE(fockcommon%stress_ikpt)
1495 : end if
1496 : end if
1497 :
1498 44402 : if (psps%usepaw==1) then
1499 8997 : call pawcprj_free(cwaveprj)
1500 : end if
1501 47331 : ABI_FREE(cwaveprj)
1502 :
1503 44402 : if(gpu_option == ABI_GPU_KOKKOS) then
1504 : #if defined HAVE_GPU && defined HAVE_YAKL
1505 : ABI_FREE_MANAGED(cwavef)
1506 : #endif
1507 : else
1508 44402 : ABI_FREE(cwavef)
1509 : end if
1510 :
1511 44402 : ABI_FREE(lambda)
1512 44402 : ABI_FREE(occblock)
1513 44402 : ABI_FREE(weight)
1514 44402 : ABI_FREE(enlout)
1515 44402 : if (usexg==1) then
1516 1160 : ABI_FREE(enlout_2d)
1517 1160 : ABI_FREE(enlout_2d_stress)
1518 : end if
1519 : #ifdef HAVE_OPENMP_OFFLOAD
1520 : !$OMP TARGET EXIT DATA MAP(delete:ffnl) IF(gpu_option==ABI_GPU_OPENMP)
1521 : #endif
1522 44402 : ABI_FREE(ffnl)
1523 44402 : ABI_SFREE(ffnl_kphq)
1524 44402 : ABI_FREE(kg_k)
1525 44402 : ABI_FREE(kpg_k)
1526 44402 : ABI_SFREE(kpg_kphq)
1527 44402 : ABI_FREE(ylm_k)
1528 44402 : ABI_FREE(ylmgr_k)
1529 44402 : ABI_FREE(ph3d)
1530 44402 : ABI_SFREE(ph3d_kphq)
1531 44402 : if (stress_needed==1) then
1532 44173 : ABI_FREE(kstr1)
1533 44173 : ABI_FREE(kstr2)
1534 44173 : ABI_FREE(kstr3)
1535 44173 : ABI_FREE(kstr4)
1536 44173 : ABI_FREE(kstr5)
1537 44173 : ABI_FREE(kstr6)
1538 44173 : if (use_gbt /= 0) then
1539 64 : ABI_FREE(kstr1_kphq)
1540 64 : ABI_FREE(kstr2_kphq)
1541 64 : ABI_FREE(kstr3_kphq)
1542 64 : ABI_FREE(kstr4_kphq)
1543 64 : ABI_FREE(kstr5_kphq)
1544 64 : ABI_FREE(kstr6_kphq)
1545 : end if
1546 : end if
1547 44402 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
1548 : #ifdef HAVE_OPENMP_OFFLOAD
1549 : !$OMP TARGET EXIT DATA MAP(delete:ffnl_str) IF(gpu_option==ABI_GPU_OPENMP)
1550 : #endif
1551 54 : ABI_FREE(ffnl_str)
1552 : end if
1553 :
1554 51299 : call timab(927,2,tsec)
1555 :
1556 : end do ! End k point loop
1557 :
1558 : end do ! End loop over spins
1559 :
1560 6290 : call timab(928,1,tsec)
1561 :
1562 6290 : if (usexg==1) then
1563 166 : if (xg_nonlop%paw) call xg_nonlop_destroy_Dij(xg_nonlop)
1564 : end if
1565 :
1566 : !Stress is equal to dE/d_strain * (1/ucvol)
1567 44030 : npsstr(:)=npsstr(:)/gs_hamk%ucvol
1568 :
1569 : !Parallel case: accumulate (n,k) contributions
1570 : if (xmpi_paral==1) then
1571 : ! Forces
1572 6290 : if (optfor==1) then
1573 3226 : call timab(65,1,tsec)
1574 3226 : call xmpi_sum(grnl,spaceComm,ierr)
1575 3226 : call timab(65,2,tsec)
1576 3226 : if ((usefock_loc).and.(psps%usepaw==1)) then
1577 10 : call xmpi_sum(fockcommon%forces,spaceComm,ierr)
1578 : end if
1579 : end if
1580 : ! Stresses
1581 6290 : if (stress_needed==1) then
1582 6259 : call timab(65,1,tsec)
1583 6259 : call xmpi_sum(kinstr,spaceComm,ierr)
1584 6259 : call xmpi_sum(npsstr,spaceComm,ierr)
1585 6259 : if (usevxctau>0) then
1586 37 : call xmpi_sum(mggastr,spaceComm,ierr)
1587 : end if
1588 6259 : if (usefock_loc) then
1589 121 : if (fockcommon%optstr) then
1590 121 : call xmpi_sum(fockcommon%stress,spaceComm,ierr)
1591 : end if
1592 : end if
1593 6259 : call timab(65,2,tsec)
1594 : end if
1595 : end if
1596 :
1597 : !Do final normalizations and symmetrizations of stress tensor contributions
1598 6290 : if (stress_needed==1) then
1599 6259 : renorm_factor=-(two_pi**2)/effmass_free/gs_hamk%ucvol
1600 43813 : kinstr(:)=kinstr(:)*renorm_factor
1601 6259 : if (nsym>1) then
1602 4850 : call stresssym(gs_hamk%gprimd,nsym,kinstr,symrec)
1603 4850 : call stresssym(gs_hamk%gprimd,nsym,npsstr,symrec)
1604 4850 : if (usevxctau>0) then
1605 37 : call stresssym(gs_hamk%gprimd,nsym,mggastr,symrec)
1606 : end if
1607 4850 : if (usefock_loc) then
1608 121 : if (fockcommon%optstr) then
1609 121 : call stresssym(gs_hamk%gprimd,nsym,fockcommon%stress,symrec)
1610 : end if
1611 : end if
1612 : end if
1613 : end if
1614 :
1615 : !Need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted)
1616 6290 : if (psps%usepaw==1.and.usecprj_local==1) then
1617 106 : call pawcprj_reorder(cprj,gs_hamk%atindx1)
1618 : end if
1619 :
1620 : !Deallocate temporary space
1621 6290 : if (usevxctau>0) then
1622 37 : ABI_FREE(vxctaulocal)
1623 : end if
1624 6290 : call gs_hamk%free()
1625 6290 : if (usefock_loc) then
1626 122 : fockcommon%use_ACE=use_ACE_old
1627 : end if
1628 :
1629 6290 : call timab(928,2,tsec) ; call timab(920,-2,tsec)
1630 : ABI_NVTX_END_RANGE()
1631 :
1632 12580 : end subroutine forstrnps
1633 : !!***
1634 :
1635 : !----------------------------------------------------------------------
1636 :
1637 : !!****f* ABINIT/nres2vres
1638 : !!
1639 : !! NAME
1640 : !! nres2vres
1641 : !!
1642 : !! FUNCTION
1643 : !! Convert a density residual into a potential residual
1644 : !! using a first order formula:
1645 : !! V^res(r)=dV/dn.n^res(r)
1646 : !! =V_hartree(n^res)(r) + Kxc.n^res(r)
1647 : !! INPUTS
1648 : !! dtset <type(dataset_type)>=all input variables in this dataset
1649 : !! | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver
1650 : !! | natom= number of atoms in cell
1651 : !! | nspden=number of spin-density components
1652 : !! | ntypat=number of atom types
1653 : !! | typat(natom)=type (integer) for each atom
1654 : !! gsqcut=cutoff value on G**2 for sphere inside fft box
1655 : !! izero=if 1, unbalanced components of Vhartree(g) are set to zero
1656 : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
1657 : !! mpi_enreg=information about MPI parallelization
1658 : !! my_natom=number of atoms treated by current processor
1659 : !! nfft=(effective) number of FFT grid points (for this processor)
1660 : !! ngfft(18)=contain all needed information about 3D FFT
1661 : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
1662 : !! nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
1663 : !! nresid(nfft,nspden)= the input density residual
1664 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
1665 : !! optnc=option for non-collinear magnetism (nspden=4):
1666 : !! 1: the whole 2x2 Vres matrix is computed
1667 : !! 2: only Vres^{11} and Vres^{22} are computed
1668 : !! optxc=0 if LDA part of XC kernel has only to be taken into account (even for GGA)
1669 : !! 1 if XC kernel has to be fully taken into
1670 : !! -1 if XC kernel does not have to be taken into account
1671 : !! pawang <type(pawang_type)>=paw angular mesh and related data
1672 : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
1673 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
1674 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1675 : !! rhor(nfft,nspden)=electron density in real space
1676 : !! (used only if Kxc was not computed before)
1677 : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
1678 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
1679 : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
1680 : !! xcctau3d(n3xccc)=3D core electron kinetic energy density for XC core correction (bohr^-3)
1681 : !! xred(3,natom)=reduced dimensionless atomic coordinates
1682 : !!
1683 : !! === optional inputs ===
1684 : !! vxc(cplex*nfft,nspden)=XC GS potential
1685 : !!
1686 : !! OUTPUT
1687 : !! vresid(nfft,nspden)= the output potential residual
1688 : !!
1689 : !! SOURCE
1690 :
1691 6482 : subroutine nres2vres(dtset,gsqcut,izero,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
1692 6482 : & nkxc,nresid,n3xccc,optnc,optxc,pawang,pawfgrtab,pawrhoij,pawtab,&
1693 6482 : & rhor,rprimd,usepaw,vresid,xccc3d,xred,&
1694 713 : & vxc, xcctau3d)
1695 :
1696 : !Arguments ------------------------------------
1697 : !scalars
1698 : integer,intent(in) :: izero,my_natom,n3xccc,nfft,nkxc,optnc,optxc,usepaw
1699 : real(dp),intent(in) :: gsqcut
1700 : type(MPI_type),intent(in) :: mpi_enreg
1701 : type(dataset_type),intent(in) :: dtset
1702 : type(pawang_type),intent(in) :: pawang
1703 : !arrays
1704 : integer,intent(in) :: ngfft(18)
1705 : real(dp),intent(in) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden)
1706 : real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3),xccc3d(n3xccc),xred(3,dtset%natom)
1707 : real(dp),intent(in),optional :: xcctau3d(n3xccc*dtset%usekden)
1708 : real(dp),intent(inout) :: nhat(nfft,dtset%nspden*usepaw)
1709 : real(dp),intent(out) :: vresid(nfft,dtset%nspden)
1710 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*usepaw)
1711 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*usepaw)
1712 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw)
1713 : real(dp),intent(in) :: vxc(nfft,dtset%nspden) !FR TODO:cplex?
1714 :
1715 : !Local variables-------------------------------
1716 : !scalars
1717 : integer :: cplex,ider,idir,ipert,ispden,nhatgrdim,nkxc_cur,option,me,nproc,comm,usexcnhat
1718 : logical :: has_nkxc_gga,non_magnetic_xc
1719 : real(dp) :: dum,energy,m_norm_min,ucvol,vxcavg
1720 : character(len=500) :: message
1721 : type(xcdata_type) :: xcdata
1722 : !arrays
1723 : integer :: nk3xc
1724 : real(dp) :: gmet(3,3),gprimd(3,3),qq(3),rmet(3,3)
1725 6482 : real(dp),allocatable :: dummy(:),kxc_cur(:,:),nhatgr(:,:,:)
1726 6482 : real(dp),allocatable :: nresg(:,:),rhor0(:,:),vhres(:)
1727 : ! *************************************************************************
1728 :
1729 : !Compatibility tests:
1730 6482 : has_nkxc_gga=(nkxc==7.or.nkxc==19)
1731 :
1732 6482 : if(optxc<-1.or.optxc>1)then
1733 0 : write(message,'(a,i0)')' Wrong value for optxc ',optxc
1734 0 : ABI_BUG(message)
1735 : end if
1736 :
1737 6482 : if((optnc/=1.and.optnc/=2).or.(dtset%nspden/=4.and.optnc/=1))then
1738 0 : write(message,'(a,i0)')' Wrong value for optnc ',optnc
1739 0 : ABI_BUG(message)
1740 : end if
1741 :
1742 6482 : if(dtset%icoulomb==1.and.optxc/=-1)then
1743 0 : write(message,'(a)')' This routine is not compatible with icoulomb==1 and optxc/=-1 !'
1744 0 : ABI_BUG(message)
1745 : end if
1746 :
1747 6482 : if(dtset%nspden==4.and.dtset%xclevel==2.and.optxc==1.and.(.not.has_nkxc_gga))then
1748 0 : ABI_ERROR(' Wrong values for optxc and nkxc !')
1749 : end if
1750 :
1751 6482 : qq=zero
1752 6482 : nkxc_cur=0
1753 6482 : m_norm_min=EPSILON(0.0_dp)**2
1754 16633 : usexcnhat=0;if (usepaw==1) usexcnhat=maxval(pawtab(1:dtset%ntypat)%usexcnhat)
1755 6482 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
1756 6482 : if (dtset%xclevel==1.or.optxc==0) nkxc_cur= 2*min(dtset%nspden,2)-1 ! LDA: nkxc=1,3
1757 6482 : if (dtset%xclevel==2.and.optxc==1)nkxc_cur=12*min(dtset%nspden,2)-5 ! GGA: nkxc=7,19
1758 19446 : ABI_MALLOC(vhres,(nfft))
1759 :
1760 : !Compute different geometric tensor, as well as ucvol, from rprimd
1761 6482 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1762 :
1763 : !Compute density residual in reciprocal space
1764 6482 : if (dtset%icoulomb==0) then
1765 19446 : ABI_MALLOC(nresg,(2,nfft))
1766 12964 : ABI_MALLOC(dummy,(nfft))
1767 332247771 : dummy(:)=nresid(:,1)
1768 6482 : call fourdp(1,nresg,dummy,-1,mpi_enreg,nfft,1,ngfft,0)
1769 6482 : ABI_FREE(dummy)
1770 : end if
1771 :
1772 : !For GGA, has to recompute gradients of nhat
1773 6482 : nhatgrdim=0
1774 6482 : if ((nkxc==nkxc_cur.and.has_nkxc_gga).or.(optxc==-1.and.has_nkxc_gga).or.&
1775 : & (optxc/=-1.and.nkxc/=nkxc_cur)) then
1776 1 : if (usepaw==1.and.dtset%xclevel==2.and.usexcnhat>0.and.dtset%pawnhatxc>0) then
1777 0 : nhatgrdim=1
1778 0 : ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3))
1779 0 : ider=1;cplex=1;ipert=0;idir=0
1780 : call pawmknhat(dum,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
1781 : & nfft,ngfft,nhatgrdim,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,&
1782 : & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qq,rprimd,ucvol,dtset%usewvl,xred,&
1783 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1784 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1785 0 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1786 : else
1787 1 : ABI_MALLOC(nhatgr,(0,0,0))
1788 : end if
1789 : else
1790 6481 : ABI_MALLOC(nhatgr,(0,0,0))
1791 : end if
1792 :
1793 6482 : ABI_MALLOC(dummy,(0))
1794 : !First case: Kxc has already been computed
1795 : !-----------------------------------------
1796 6482 : if (nkxc==nkxc_cur.or.optxc==-1) then
1797 :
1798 : ! Compute VH(n^res)(r)
1799 6481 : if (dtset%icoulomb == 0) then
1800 : call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,&
1801 6481 : &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres)
1802 : else
1803 0 : comm=mpi_enreg%comm_cell
1804 0 : nproc=xmpi_comm_size(comm)
1805 0 : me=xmpi_comm_rank(comm)
1806 : call psolver_hartree(energy, (/ rprimd(1,1) / dtset%ngfft(1), &
1807 : & rprimd(2,2) / dtset%ngfft(2), rprimd(3,3) / dtset%ngfft(3) /), dtset%icoulomb, &
1808 : & me, comm, dtset%nfft, dtset%ngfft(1:3), nproc, dtset%nscforder, dtset%nspden, &
1809 0 : & nresid(:,1), vhres, dtset%usewvl)
1810 : end if
1811 :
1812 : ! Compute Kxc(r).n^res(r)
1813 6481 : if (optxc/=-1) then
1814 :
1815 : ! Collinear magnetism or non-polarized
1816 6481 : if (dtset%nspden/=4) then
1817 : !Note: imposing usexcnhat=1 avoid nhat to be substracted
1818 : call dfpt_mkvxc(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,&
1819 6392 : & nkxc,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy)
1820 : else
1821 : !FR call routine for Non-collinear magnetism
1822 356 : ABI_MALLOC(rhor0,(nfft,dtset%nspden))
1823 3635833 : rhor0(:,:)=rhor(:,:)-nresid(:,:)
1824 : !Note: imposing usexcnhat=1 avoid nhat to be substracted
1825 : call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,&
1826 89 : & nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d)
1827 89 : ABI_FREE(rhor0)
1828 : end if
1829 :
1830 : else
1831 1 : vresid=zero
1832 : end if
1833 :
1834 : end if
1835 :
1836 : !2nd case: Kxc has to be computed
1837 : !--------------------------------
1838 6482 : if (nkxc/=nkxc_cur.and.optxc/=-1) then
1839 :
1840 : ! Has to use the "initial" density to compute Kxc
1841 4 : ABI_MALLOC(rhor0,(nfft,dtset%nspden))
1842 6753 : rhor0(:,:)=rhor(:,:)-nresid(:,:)
1843 :
1844 : ! Compute VH(n^res) and XC kernel (Kxc) together
1845 4 : ABI_MALLOC(kxc_cur,(nfft,nkxc_cur))
1846 :
1847 1 : option=2;if (dtset%xclevel==2.and.optxc==0) option=12
1848 :
1849 : call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,&
1850 1 : &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres)
1851 1 : call xcdata_init(xcdata,dtset=dtset)
1852 :
1853 : ! To be adjusted for the call to rhotoxc
1854 1 : nk3xc=1
1855 : call rhotoxc(energy,dum,kxc_cur,mpi_enreg,nfft,ngfft,&
1856 : & nhat,usepaw,nhatgr,nhatgrdim,nkxc_cur,nk3xc,non_magnetic_xc,n3xccc,option,&
1857 : & rhor0,rprimd,usexcnhat,vresid,vxcavg,xccc3d,xcdata,vhartr=vhres,&
1858 1 : & xcctau3d=xcctau3d) !vresid=work space
1859 1 : if (dtset%nspden/=4) then
1860 1 : ABI_FREE(rhor0)
1861 : end if
1862 :
1863 : ! Compute Kxc(r).n^res(r)
1864 :
1865 1 : if (dtset%nspden/=4) then
1866 : ! Collinear magnetism or non-polarized
1867 : call dfpt_mkvxc(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,&
1868 1 : & nkxc_cur,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy)
1869 : else
1870 : ! Non-collinear magnetism
1871 0 : ABI_MALLOC(rhor0,(nfft,dtset%nspden))
1872 0 : rhor0(:,:)=rhor(:,:)-nresid(:,:)
1873 : call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,&
1874 0 : & nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d)
1875 0 : ABI_FREE(rhor0)
1876 : end if
1877 :
1878 2 : ABI_FREE(kxc_cur)
1879 : end if
1880 :
1881 : !if (nhatgrdim>0) then
1882 6482 : ABI_FREE(nhatgr)
1883 : !end if
1884 :
1885 : !Assemble potential residual: V^res(r)=VH(n^res)(r) + Kxc(r).n^res(r)
1886 : !--------------------------------------------------------------------
1887 13573 : do ispden=1,dtset%nspden/optnc
1888 353303490 : vresid(:,ispden)=vresid(:,ispden)+vhres(:)
1889 : end do
1890 :
1891 6482 : if (dtset%icoulomb==0) then
1892 6482 : ABI_FREE(nresg)
1893 : end if
1894 6482 : ABI_FREE(vhres)
1895 6482 : ABI_FREE(dummy)
1896 :
1897 12964 : end subroutine nres2vres
1898 : !!***
1899 :
1900 : !----------------------------------------------------------------------
1901 :
1902 : !!****f* ABINIT/stress_mGGA
1903 : !!
1904 : !! NAME
1905 : !! stress_mGGA
1906 : !!
1907 : !! FUNCTION
1908 : !! Compute non-local metaGGA contribution to stress tensor (i.e. contribution from Div.(Vtau.Grad))
1909 : !! of (one block of) band Psi_n
1910 : !! This contribution writes: -1/m_e * Sum_n[ f_n . Re{<Grad_alpha(Psi_n)|V_tau|Grad_beta(Psi_n)>} ]
1911 : !!
1912 : !! INPUTS
1913 : !! cwavef(2,npw_k*my_nspinor*ndat)=planewave coefficients of wavefunction.
1914 : !! effmass_free=effective mass for electrons (1. in common case)
1915 : !! gbound_k(2*mgfft+4)=sphere boundary info
1916 : !! gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1917 : !! istwf_k=input parameter that describes the storage of wfs
1918 : !! kg_k(3,npw_k)=G vec coordinates wrt recip lattice transl.
1919 : !! kpt(3)=current k point
1920 : !! mgfft=maximum single fft dimension
1921 : !! mpi_enreg=information about MPI parallelization
1922 : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
1923 : !! ndat=number of wave functions used to accumulate stress tensor (block of bands)
1924 : !! ngfft(18)=contain all needed information about 3D FFT
1925 : !! npw_k=number of planewaves in basis for given k point.
1926 : !! nvtau=number of spin components of vxctaulocal
1927 : !! n4,n5,n6=for dimensionning of vxctaulocal
1928 : !! occ(ndat)=occupancies of bands at various k points
1929 : !! ucvol=unit cell volume in bohr**3
1930 : !! vxctaulocal(n4,n5,n6,nvtau,4)= local potential corresponding to the derivative of XC energy
1931 : !! with respect to kinetic energy density, in real space, on the augmented fft grid.
1932 : !! This array contains also the gradient of vxctaulocal (gvxctaulocal) in vxctaulocal(:,:,:,:,2:4).
1933 : !! wtk=weights associated with current k-point
1934 : !! [gpu_option]=GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU) [default=0]
1935 : !!
1936 : !! OUTPUT
1937 : !!
1938 : !! SIDE EFFECTS
1939 : !! mggastr(6)=metaGGA contribution to stress tensor updated
1940 : !!
1941 : !! SOURCE
1942 :
1943 197 : subroutine stress_mGGA(mggastr,cwavef,effmass_free,gbound_k,gprimd,istwf_k,kg_k,kpt,mgfft, &
1944 197 : & mpi_enreg,my_nspinor,ndat,ngfft,npw_k,nvtau,n4,n5,n6,occ,ucvol,vxctaulocal,wtk, &
1945 : & gpu_option) ! optional
1946 :
1947 : !Arguments ------------------------------------
1948 : !scalars
1949 : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvtau,n4,n5,n6
1950 : integer,intent(in),optional :: gpu_option
1951 : real(dp),intent(in) :: effmass_free,ucvol,wtk
1952 : type(MPI_type),intent(in) :: mpi_enreg
1953 : !arrays
1954 : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
1955 : real(dp),intent(inout),target :: cwavef(2,npw_k*my_nspinor*ndat)
1956 : real(dp),intent(in) :: gprimd(3,3),kpt(3),occ(ndat)
1957 : real(dp),intent(inout) :: vxctaulocal(n4,n5,n6,nvtau,4)
1958 : real(dp),intent(inout) :: mggastr(6)
1959 :
1960 : !Local variables-------------------------------
1961 : !scalars
1962 : integer,parameter :: opt_fourwf=2,tim_fourwf=1
1963 : integer :: gpu_option_,ia,ib,idat,idir,ierr,ipw,ispinor,mu,nspinortot
1964 : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
1965 : real(dp) :: gp2pi1,gp2pi2,gp2pi3,kg_k_cart,kpt_cart,renorm_factor,weight_dum=1
1966 : !arrays
1967 : integer,parameter :: voigt1(6)=[1,2,3,3,3,2],voigt2(6)=[1,2,3,2,1,1]
1968 394 : real(dp) :: dotr(ndat),doti(ndat),my_mggastr(6)
1969 197 : real(dp),allocatable,target :: gcwavef(:,:,:),vtau_gcwavef(:,:,:)
1970 197 : real(dp),allocatable :: weight_array(:),work(:,:,:,:)
1971 197 : real(dp),pointer :: gcwavef_ndat(:,:,:,:),vtau_gcwavef_ndat(:,:,:,:)
1972 197 : real(dp),pointer :: my_cwavef(:,:)
1973 :
1974 : ! *********************************************************************
1975 :
1976 197 : if (nvtau/=1) then
1977 0 : ABI_BUG("mGGA potential Vtau not compatible with non-collinear magnetism!")
1978 : end if
1979 :
1980 : !Some inits
1981 197 : gpu_option_=0;if (present(gpu_option)) gpu_option_=gpu_option
1982 197 : renorm_factor=-one/effmass_free/ucvol
1983 197 : my_mggastr(:)=zero
1984 591 : ABI_MALLOC(weight_array,(ndat))
1985 509 : weight_array(1:ndat)=wtk*occ(1:ndat)
1986 :
1987 : !Parallelization over spinors
1988 197 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
1989 197 : if (mpi_enreg%paral_spinor==0) then
1990 197 : nspinor1TreatedByThisProc=.true.
1991 197 : nspinor2TreatedByThisProc=(nspinortot==2)
1992 : else
1993 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
1994 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
1995 : end if
1996 :
1997 : !Manage memory
1998 197 : if (nspinortot==1) then
1999 132 : my_cwavef => cwavef
2000 : else ! nspinortot==2
2001 195 : ABI_MALLOC(my_cwavef,(2,npw_k*ndat))
2002 : end if
2003 985 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
2004 788 : ABI_MALLOC(gcwavef,(2,npw_k*ndat,3))
2005 394 : ABI_MALLOC(vtau_gcwavef,(2,npw_k*ndat,3))
2006 985 : call c_f_pointer(c_loc(gcwavef(1,1,1)),gcwavef_ndat,[2,npw_k,ndat,3])
2007 985 : call c_f_pointer(c_loc(vtau_gcwavef(1,1,1)),vtau_gcwavef_ndat,[2,npw_k,ndat,3])
2008 : #ifdef HAVE_OPENMP_OFFLOAD
2009 : !$OMP TARGET ENTER DATA MAP(alloc:gcwavef,vtau_gcwavef,dotr,doti) IF(gpu_option_==ABI_GPU_OPENMP)
2010 : !$OMP TARGET ENTER DATA MAP(alloc:my_cwavef) IF(gpu_option_==ABI_GPU_OPENMP .and. nspinortot==2)
2011 : !$OMP TARGET ENTER DATA MAP(to:cwavef) IF(gpu_option_==ABI_GPU_OPENMP)
2012 : #endif
2013 :
2014 : !Loop over spinors (if any)
2015 459 : do ispinor=1,my_nspinor
2016 :
2017 : ! Select spinor component of WF
2018 262 : if (nspinortot==2) then
2019 130 : if (ispinor==1.and.nspinor1TreatedByThisProc) then
2020 65 : if(gpu_option_==ABI_GPU_OPENMP) then
2021 : #ifdef HAVE_OPENMP_OFFLOAD
2022 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:my_cwavef,cwavef) PRIVATE(idat)
2023 : do idat=1,ndat
2024 : !$OMP PARALLEL DO PRIVATE(ipw)
2025 : do ipw=1,npw_k
2026 : my_cwavef(1,ipw+(idat-1)*npw_k)=cwavef(1,ipw+(idat-1)*my_nspinor*npw_k)
2027 : my_cwavef(2,ipw+(idat-1)*npw_k)=cwavef(2,ipw+(idat-1)*my_nspinor*npw_k)
2028 : end do
2029 : end do
2030 : #endif
2031 : else
2032 169 : do idat=1,ndat
2033 67249 : do ipw=1,npw_k
2034 402584 : my_cwavef(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k)
2035 : end do
2036 : end do
2037 : end if ! gpu_option_
2038 65 : else if (ispinor==2.and.nspinor2TreatedByThisProc) then
2039 65 : if(gpu_option_==ABI_GPU_OPENMP) then
2040 : #ifdef HAVE_OPENMP_OFFLOAD
2041 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:my_cwavef,cwavef) PRIVATE(idat)
2042 : do idat=1,ndat
2043 : !$OMP PARALLEL DO PRIVATE(ipw)
2044 : do ipw=1,npw_k
2045 : my_cwavef(1,ipw+(idat-1)*npw_k)=cwavef(1,ipw+(idat-1)*my_nspinor*npw_k+npw_k)
2046 : my_cwavef(2,ipw+(idat-1)*npw_k)=cwavef(2,ipw+(idat-1)*my_nspinor*npw_k+npw_k)
2047 : end do
2048 : end do
2049 : #endif
2050 : else
2051 169 : do idat=1,ndat
2052 67249 : do ipw=1,npw_k
2053 402584 : my_cwavef(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k+npw_k)
2054 : end do
2055 : end do
2056 : end if ! gpu_option_
2057 : else
2058 : cycle
2059 : end if
2060 : end if
2061 :
2062 : ! Loop over cartesian directions
2063 1048 : do idir=1,3
2064 786 : gp2pi1=gprimd(idir,1)*two_pi
2065 786 : gp2pi2=gprimd(idir,2)*two_pi
2066 786 : gp2pi3=gprimd(idir,3)*two_pi
2067 786 : kpt_cart=gp2pi1*kpt(1)+gp2pi2*kpt(2)+gp2pi3*kpt(3)
2068 :
2069 : ! Compute grad of WF (multiplication by 2pi i (G+k)_idir in reciprocal space)
2070 786 : if(gpu_option_==ABI_GPU_OPENMP) then
2071 : #ifdef HAVE_OPENMP_OFFLOAD
2072 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:my_cwavef,gcwavef,kg_k) PRIVATE(idat)
2073 : do idat=1,ndat
2074 : !$OMP PARALLEL DO PRIVATE(ipw,kg_k_cart)
2075 : do ipw=1,npw_k
2076 : kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
2077 : gcwavef(1,ipw+(idat-1)*npw_k,idir)= my_cwavef(2,ipw+(idat-1)*npw_k)*kg_k_cart
2078 : gcwavef(2,ipw+(idat-1)*npw_k,idir)=-my_cwavef(1,ipw+(idat-1)*npw_k)*kg_k_cart
2079 : end do
2080 : end do
2081 : #endif
2082 : else
2083 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(idat,ipw,kg_k_cart)
2084 2034 : do idat=1,ndat
2085 909816 : do ipw=1,npw_k
2086 907782 : kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
2087 907782 : gcwavef(1,ipw+(idat-1)*npw_k,idir)= my_cwavef(2,ipw+(idat-1)*npw_k)*kg_k_cart
2088 909030 : gcwavef(2,ipw+(idat-1)*npw_k,idir)=-my_cwavef(1,ipw+(idat-1)*npw_k)*kg_k_cart
2089 : end do
2090 : end do
2091 : end if
2092 :
2093 : ! Compute vxctaulocal*(grad of WF) in reciprocal space
2094 : call fourwf(1,vxctaulocal(:,:,:,:,1),gcwavef(:,:,idir),vtau_gcwavef(:,:,idir), &
2095 : & work,gbound_k,gbound_k,istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft, &
2096 : & npw_k,npw_k,n4,n5,n6,opt_fourwf,tim_fourwf,weight_dum,weight_dum, &
2097 1048 : & weight_array_r=weight_array,gpu_option=gpu_option_)
2098 :
2099 : end do ! idir
2100 :
2101 : ! Accumulate stress tensor components (Re{<Grad_alpha(Psi_n)|V_tau.Grad_beta(Psi_n)>})
2102 2031 : do mu=1,6
2103 1572 : ia=voigt1(mu) ; ib=voigt2(mu)
2104 : call dotprod_g_batch_full(dotr,doti,istwf_k,npw_k,ndat,1, &
2105 : & gcwavef_ndat(:,:,:,ia),vtau_gcwavef_ndat(:,:,:,ib), &
2106 1572 : & mpi_enreg%me_g0,mpi_enreg%comm_fft,gpu_option=gpu_option_)
2107 : #ifdef HAVE_OPENMP_OFFLOAD
2108 : !$OMP TARGET UPDATE FROM(dotr,doti) IF(gpu_option_==ABI_GPU_OPENMP)
2109 : #endif
2110 4330 : my_mggastr(mu)=my_mggastr(mu) + renorm_factor*sum(weight_array(1:ndat)*dotr(1:ndat))
2111 : end do
2112 :
2113 : end do ! ispinor
2114 :
2115 : #ifdef HAVE_OPENMP_OFFLOAD
2116 : !$OMP TARGET EXIT DATA MAP(delete:cwavef,gcwavef,vtau_gcwavef,dotr,doti) IF(gpu_option_==ABI_GPU_OPENMP)
2117 : !$OMP TARGET EXIT DATA MAP(delete:my_cwavef) IF(gpu_option_==ABI_GPU_OPENMP .and. nspinortot==2)
2118 : #endif
2119 : !Release memory
2120 197 : if (nspinortot==2) then
2121 65 : ABI_FREE(my_cwavef)
2122 : end if
2123 197 : ABI_FREE(gcwavef)
2124 197 : ABI_FREE(vtau_gcwavef)
2125 197 : ABI_FREE(weight_array)
2126 197 : ABI_FREE(work)
2127 :
2128 : !Take into account MPI parallelism (bands, spinors)
2129 197 : call xmpi_sum(my_mggastr,mpi_enreg%comm_bandspinor ,ierr)
2130 :
2131 : !Final accumulation of stresses
2132 1379 : mggastr(1:6) = mggastr(1:6) + my_mggastr(1:6)
2133 :
2134 197 : end subroutine stress_mGGA
2135 : !!***
2136 :
2137 : end module m_forstr
2138 : !!***
|
