Branch data 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, initylmg_k
61 : : use m_xchybrid, only : xchybrid_ncpp_cc
62 : : use m_kg, only : mkkpg
63 : : use m_hamiltonian, only : 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 : 6364 : subroutine forstr(atindx1,cg,cprj,diffor,dtefield,dtset,eigen,electronpositron,energies,favg,fcart,fock,&
268 : 6364 : & forold,gred,grchempottn,grcondft,gresid,grewtn,grhf,grvdw,grxc,gsqcut,extfpmd,indsym,&
269 : 6364 : & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,&
270 : 6364 : & nfftf,ngfftf,ngrvdw,nhat,nkxc,npwarr,&
271 : 6364 : & ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,&
272 : 6364 : & pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,psps,rhog,rhor,rprimd,stress_needed,&
273 : 6364 : & strscondft,strsxc,strten,symrec,synlgr,ucvol,usecprj,usevxctau,vhartr,vpsp,&
274 : 6364 : & 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 : 6364 : 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 : 6364 : real(dp),allocatable :: grnl(:),vlocal(:,:),vxc_hf(:,:),xcart(:,:),ylmbz(:,:),ylmgrbz(:,:,:)
336 : 6364 : real(dp), contiguous, pointer :: resid(:,:)
337 : : ! *************************************************************************
338 : :
339 : 6364 : call timab(910,1,tsec)
340 : : ABI_NVTX_START_RANGE(NVTX_FORSTR)
341 : 6364 : call timab(911,1,tsec)
342 : :
343 : : !Do nothing if nothing is required
344 [ + + + - ]: 6364 : 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 [ + - ]: 6364 : if (dtset%usewvl==0) then
348 [ + + + - : 6364 : 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 [ + + + - : 6364 : 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 [ + + + + : 19092 : ABI_MALLOC(grnl,(3*dtset%natom*optfor))
+ - + - +
+ - + -
+ ]
362 [ + + ]: 40711 : 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 [ + + + - ]: 6364 : 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 [ + - ]: 6321 : else if (dtset%usewvl==0) then
370 : 6321 : occopt_=0 ! This means that occ are now fixed
371 [ + + + - ]: 6321 : 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 : 6321 : usexg = 0
421 [ + + ]: 6321 : 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 : 6321 : & 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 : 6364 : call timab(911,2,tsec)
437 : 6364 : call timab(912,1,tsec)
438 : :
439 : : !PAW: add gradients due to Dij derivatives to non-local term
440 [ + + ]: 6364 : if (psps%usepaw==1) then
441 [ - + + - : 5760 : ABI_MALLOC(vlocal,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
442 : :
443 : : !$OMP PARALLEL DO COLLAPSE(2)
444 [ + + ]: 3137 : do ispden=1,min(dtset%nspden,2)
445 [ + + ]: 59710415 : do ifft=1,nfftf
446 : 59708975 : vlocal(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
447 : : end do
448 : : end do
449 [ + + ]: 1440 : if (dtset%nspden==4) then
450 : : !$OMP PARALLEL DO COLLAPSE(2)
451 [ + + ]: 174 : do ispden=3,4
452 [ + + ]: 2319660 : do ifft=1,nfftf
453 : 2319602 : vlocal(ifft,ispden)=vxc(ifft,ispden)
454 : : end do
455 : : end do
456 : : end if
457 : 1440 : 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 : 1440 : optgr=optfor;optgr2=0;optstr=stress_needed;optstr2=0
462 [ - + ]: 1440 : 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 : 1440 : & 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 [ - + ]: 1440 : ABI_FREE(vlocal)
471 : :
472 : : end if
473 : 6364 : call timab(912,2,tsec)
474 : 6364 : call timab(913,1,tsec)
475 : :
476 : : !==========================================================================
477 : : !Here compute forces (if required)
478 : : !==========================================================================
479 [ + + ]: 6364 : if (optfor==1) then
480 : : apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. &
481 [ + + + - : 3279 : & 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 [ - + + - : 2860 : ABI_MALLOC(resid,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
487 [ - + ]: 715 : option=0; if (dtset%densfor_pred<0) option=1
488 [ + + + + ]: 715 : 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 : 715 : & rhor,rprimd,psps%usepaw,resid,xccc3d,xred,vxc,xcctau3d)
492 : : else
493 : 2564 : 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 : 3279 : & electronpositron=electronpositron)
502 : :
503 [ + + ]: 3279 : if (apply_residual) then
504 [ - + ]: 715 : ABI_FREE(resid)
505 : : end if
506 : : end if
507 : :
508 : 6364 : call timab(913,2,tsec)
509 : 6364 : call timab(914,1,tsec)
510 : :
511 : : !==========================================================================
512 : : !Here compute stress tensor (if required)
513 : : !==========================================================================
514 : :
515 [ + + + - ]: 6364 : 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 [ + + + - ]: 6333 : 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 : 6333 : & electronpositron=electronpositron)
537 : : end if
538 : :
539 : : !Memory deallocation
540 [ - + ]: 6364 : ABI_FREE(grnl)
541 [ + + ]: 6364 : if (allocated(vxc_hf)) then
542 : 30 : ABI_FREE(vxc_hf)
543 : : end if
544 : :
545 : :
546 : 6364 : call timab(914,2,tsec)
547 : : ABI_NVTX_END_RANGE()
548 : 6364 : call timab(910,2,tsec)
549 : :
550 [ - + - + : 12728 : 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 : 6321 : subroutine forstrnps(cg,cprj,ecut,ecutsm,effmass_free,eigen,electronpositron,fock,&
635 : 6321 : & grnl,istwfk,kg,kinstr,npsstr,kpt,mband,mcg,mcprj,mgfft,mggastr,mkmem,mpi_enreg,mpsang,&
636 : 6321 : & mpw,my_natom,natom,nband,nfft,nfftf,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,nsym,&
637 : 6321 : & ntypat,nucdipmom,occ,optfor,paw_ij,pawfgr,pawtab,ph1d,psps,qgbt,rprimd,&
638 : 6321 : & stress_needed,symrec,typat,use_gbt,usecprj,usefock,usevxctau,vxctau,usexg,&
639 : 6321 : & 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 : 6321 : 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 : 6321 : type(xgTransposer_t) :: xgTransposer
687 : : !arrays
688 : 6321 : integer,allocatable,target :: kg_k(:,:)
689 : 6321 : 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 : 6321 : real(dp),allocatable,target :: cwavef(:,:)
696 : : #endif
697 : 6321 : real(dp),allocatable :: enlout(:),ffnl_sav(:,:,:,:),ffnl_str(:,:,:,:),ffnl_kphq_gather(:,:,:,:)
698 : 6321 : real(dp),allocatable :: ghc_dum(:,:),gprimd(:,:),kinpw_kphq_gather(:)
699 : 6321 : real(dp),allocatable :: kpg_k(:,:),kpg_kphq(:,:),kpg_k_sav(:,:),kpg_kphq_gather(:,:)
700 : 6321 : real(dp),allocatable :: kstr1(:),kstr2(:),kstr3(:),kstr4(:),kstr5(:),kstr6(:)
701 : 6321 : real(dp),allocatable :: kstr1_kphq(:),kstr2_kphq(:),kstr3_kphq(:),kstr4_kphq(:),kstr5_kphq(:),kstr6_kphq(:)
702 : 6321 : real(dp),allocatable :: lambda(:),occblock(:),ph3d(:,:,:),ph3d_kphq(:,:,:),ph3d_sav(:,:,:),ph3d_kphq_gather(:,:,:)
703 : 6321 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
704 : 6321 : real(dp),allocatable :: weight(:),ylm_k(:,:),ylm_kphq(:,:),ylmgr_k(:,:,:),ylmgr_kphq(:,:,:)
705 : 6321 : real(dp),allocatable,target :: ffnl(:,:,:,:),ffnl_kphq(:,:,:,:)
706 : 6321 : real(dp),pointer :: vxctau_ptr(:,:,:)
707 : : type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
708 : 6321 : type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
709 : 6321 : type(pawcprj_type),pointer :: cwaveprj_idat(:,:)
710 : : type(xgBlock_t) :: xgx0,xgeigen,xgforces,xgstress
711 : : type(xg_t) :: cprj_xgx0,cprj_work
712 : 6321 : real(dp),allocatable :: enlout_2d(:,:),enlout_2d_stress(:,:)
713 : 6321 : real(dp),allocatable :: cwavef_spin(:,:),enlout_spin(:)
714 : 6321 : real(dp),pointer :: cwavef_tr(:,:)
715 : : !*************************************************************************
716 : :
717 : : ABI_NVTX_START_RANGE(NVTX_FORSTRNPS)
718 : 6321 : call timab(920,1,tsec) ; call timab(921,-1,tsec)
719 : :
720 : : !Init mpicomm and me
721 [ + + ]: 6321 : if(mpi_enreg%paral_kgb==1)then
722 : 920 : spaceComm=mpi_enreg%comm_kpt
723 : 920 : me_distrb=mpi_enreg%me_kpt
724 : : else
725 : : !* In case of HF calculation
726 [ - + ]: 5401 : if (mpi_enreg%paral_hf==1) then
727 : 0 : spaceComm=mpi_enreg%comm_kpt
728 : 0 : me_distrb=mpi_enreg%me_kpt
729 : : else
730 : 5401 : spaceComm=mpi_enreg%comm_cell
731 : 5401 : me_distrb=mpi_enreg%me_cell
732 : : end if
733 : : end if
734 : :
735 : : !Some constants
736 : : ipositron=abs(electronpositron_calctype(electronpositron))
737 : 6321 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
738 : : !Smearing of plane wave kinetic energy
739 [ + + ]: 6321 : ecutsm_inv=zero;if( ecutsm>1.0d-20) ecutsm_inv=1/ecutsm
740 : : !htpisq is (1/2) (2 Pi) **2:
741 : 6321 : htpisq=0.5_dp*(two_pi)**2
742 : :
743 : : !Check that fock is present if want to use fock option
744 : 6321 : compute_gbound=.false.
745 [ + + - + ]: 12520 : usefock_loc = (usefock==1 .and. associated(fock))
746 : : !Arrays initializations
747 [ + + ]: 37473 : grnl(:)=zero
748 [ + + ]: 6321 : if (usefock_loc) then
749 : 122 : fockcommon =>fock%fock_common
750 : 122 : fockcommon%optfor=.false.
751 : 122 : fockcommon%optstr=.false.
752 : 122 : use_ACE_old=fockcommon%use_ACE
753 : 122 : fockcommon%use_ACE=0
754 : : if (fockcommon%optfor) compute_gbound=.true.
755 : : end if
756 [ + + ]: 6321 : if (stress_needed==1) then
757 : 6290 : kinstr(:)=zero;mggastr(:)=0;npsstr(:)=zero
758 [ + + ]: 6290 : if (usevxctau>0) compute_gbound=.true.
759 [ + + ]: 6290 : if (usefock_loc) then
760 : 121 : fockcommon%optstr=.TRUE.
761 [ + + ]: 847 : fockcommon%stress=zero
762 : 121 : compute_gbound=.true.
763 : : end if
764 : : end if
765 : :
766 : 6321 : usecprj_local=usecprj
767 : :
768 [ + + + + ]: 6321 : if ((usefock_loc).and.(psps%usepaw==1)) then
769 : 26 : usecprj_local=1
770 [ + + ]: 26 : if(optfor==1)then
771 : 10 : fockcommon%optfor=.true.
772 [ + - ]: 10 : if (.not.allocated(fockcommon%forces_ikpt)) then
773 [ - + + - : 40 : ABI_MALLOC(fockcommon%forces_ikpt,(3,natom,mband))
- + + - +
- + - - +
- + - + ]
774 : : end if
775 [ + - ]: 10 : if (.not.allocated(fockcommon%forces)) then
776 [ - + + - : 30 : ABI_MALLOC(fockcommon%forces,(3,natom))
+ - + - -
+ - + -
+ ]
777 : : end if
778 [ + + + + ]: 90 : fockcommon%forces=zero
779 : 10 : compute_gbound=.true.
780 : : end if
781 : : end if
782 : :
783 : : !Initialize Hamiltonian (k-independent terms)
784 : : call gs_hamk%init(psps,pawtab,nspinor,nsppol,nspden,natom,&
785 : : & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj_local,&
786 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
787 : : & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,fock=fock,&
788 : 6321 : & nucdipmom=nucdipmom,gpu_option=gpu_option,nfft_blocks=gpu_nfft_blocks)
789 [ + + + + : 252840 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
+ + ]
790 : :
791 [ + + ]: 6321 : if (usevxctau>0) then
792 [ - + + - : 259 : ABI_MALLOC(vxctaulocal,(gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,gs_hamk%nvloc,4))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + - + ]
793 : : end if
794 : :
795 [ + + ]: 6321 : if (usexg==1) then
796 : 166 : call xg_nonlop_update_weight(xg_nonlop,gs_hamk%ucvol)
797 [ + + ]: 166 : if (xg_nonlop%paw) call xg_nonlop_make_Dij(xg_nonlop,paw_ij,nsppol,gs_hamk%atindx)
798 : : end if
799 : :
800 : : !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted)
801 [ + + + + ]: 6321 : if (psps%usepaw==1.and.usecprj_local==1) then
802 : 106 : call pawcprj_reorder(cprj,gs_hamk%atindx)
803 : : end if
804 : :
805 : : !Common data for "nonlop" routine
806 : 6321 : signs=1 ; idir=0 ; ishift=0 ; tim_nonlop=4 ; tim_nonlop_prep=12
807 [ + + ]: 6321 : choice=2*optfor;if (stress_needed==1) choice=10*choice+3
808 [ + + + + ]: 6321 : if (optfor==1.and.stress_needed==1) ishift=6
809 : 6321 : nnlout=max(1,6*stress_needed+3*natom*optfor)
810 [ + + ]: 6321 : if (psps%usepaw==0) then
811 : 4881 : paw_opt=0 ; cpopt=-1
812 : : else
813 : 1440 : paw_opt=2 ; cpopt=-1+3*usecprj_local
814 : : end if
815 : :
816 : 6321 : call timab(921,2,tsec)
817 : :
818 : : !LOOP OVER SPINS
819 : 6321 : bdtot_index=0;ibg=0;icg=0
820 [ + + ]: 13250 : do isppol=1,nsppol
821 : :
822 : : ! Continue to initialize the Hamiltonian (PAW DIJ coefficients)
823 : 6929 : call gs_hamk%load_spin(isppol,with_nonlocal=.true.)
824 [ + + ]: 6929 : if (usefock_loc) fockcommon%isppol=isppol
825 : :
826 : : ! If any, set up local potential vtau on the coarse FFT mesh
827 [ + + ]: 6929 : if (usevxctau>0) then
828 : 39 : vxctau_ptr => vxctau ! This is to bepass annoying inout attribute
829 : : call gspot_transgrid_and_pack(isppol,psps%usepaw,mpi_enreg%paral_kgb,nfft,ngfft,nfftf, &
830 [ - + ]: 39 : nspden,gs_hamk%nvloc,4,pawfgr,mpi_enreg,vxctau_ptr,vxctaulocal)
831 : 39 : call gs_hamk%load_spin(isppol,vxctaulocal=vxctaulocal)
832 : : end if
833 : :
834 [ + + ]: 6929 : if (usexg==1) then
835 [ + + ]: 170 : if (xg_nonlop%paw) call xg_nonlop_set_Dij_spin(xg_nonlop,isppol)
836 : : end if
837 : :
838 : : ! Loop over k points
839 : 6929 : ikg=0
840 [ + + ]: 73776 : do ikpt=1,nkpt
841 [ + + ]: 60526 : if (usefock_loc) fockcommon%ikpt=ikpt
842 : 60526 : nband_k=nband(ikpt+(isppol-1)*nkpt)
843 : 60526 : istwf_k=istwfk(ikpt)
844 : 60526 : npw_k=npwarr(ikpt)
845 [ + + ]: 242104 : kpoint(:)=kpt(:,ikpt)
846 [ + + ]: 60526 : if (use_gbt/=0) then
847 [ + + ]: 448 : kpoint(:)=kpt(:,ikpt)-half*qgbt(:)
848 [ + + ]: 448 : kphq(:)=kpt(:,ikpt)+half*qgbt(:)
849 : : end if
850 : :
851 [ + + ]: 60526 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
852 : 15916 : bdtot_index=bdtot_index+nband_k
853 : 15916 : cycle
854 : : end if
855 : :
856 : 44610 : call timab(922,1,tsec)
857 : :
858 : 44610 : my_ikpt=mpi_enreg%my_kpttab(ikpt)
859 : : ! Parallelism over FFT and/or bands: define sizes and tabs
860 [ + + ]: 44610 : if (mpi_enreg%paral_kgb==1) then
861 : 3764 : nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
862 : 3764 : bandpp=mpi_enreg%bandpp
863 : 3764 : my_bandfft_kpt => bandfft_kpt(my_ikpt)
864 : : else
865 : 40846 : bandpp=mpi_enreg%bandpp
866 : 40846 : nblockbd=nband_k/bandpp
867 : : end if
868 : 44610 : blocksize=nband_k/nblockbd
869 : 44610 : mband_cprj=mband/mpi_enreg%nproc_band
870 : 44610 : nband_cprj_k=nband_k/mpi_enreg%nproc_band
871 : :
872 [ + - ]: 44610 : if(gpu_option == ABI_GPU_KOKKOS) then
873 : : #if defined HAVE_GPU && defined HAVE_YAKL
874 : : ABI_MALLOC_MANAGED(cwavef,(/2,npw_k*my_nspinor*blocksize/))
875 : : #endif
876 : : else
877 [ - + + - : 133830 : ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
+ - + - -
+ - + -
+ ]
878 : : end if
879 : :
880 [ + + + + ]: 44610 : if (psps%usepaw==1.and.usecprj_local==1) then
881 [ - + + - : 10589 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp))
- + + - +
- + - - +
- + - + +
+ + + ]
882 : 1517 : call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
883 : : else
884 [ - + - + ]: 43093 : ABI_MALLOC(cwaveprj,(0,0))
885 : : end if
886 : :
887 [ + + ]: 44610 : if (stress_needed==1) then
888 [ - + + - : 133143 : ABI_MALLOC(kstr1,(npw_k))
+ - + - -
+ - + -
+ ]
889 [ - + + - : 88762 : ABI_MALLOC(kstr2,(npw_k))
- + - + -
+ ]
890 [ - + + - : 88762 : ABI_MALLOC(kstr3,(npw_k))
- + - + -
+ ]
891 [ - + + - : 88762 : ABI_MALLOC(kstr4,(npw_k))
- + - + -
+ ]
892 [ - + + - : 88762 : ABI_MALLOC(kstr5,(npw_k))
- + - + -
+ ]
893 [ - + + - : 88762 : ABI_MALLOC(kstr6,(npw_k))
- + - + -
+ ]
894 [ + + ]: 44381 : if (use_gbt /= 0) then
895 [ - + + - : 224 : ABI_MALLOC(kstr1_kphq,(npw_k))
- + - + -
+ ]
896 [ - + + - : 224 : ABI_MALLOC(kstr2_kphq,(npw_k))
- + - + -
+ ]
897 [ - + + - : 224 : ABI_MALLOC(kstr3_kphq,(npw_k))
- + - + -
+ ]
898 [ - + + - : 224 : ABI_MALLOC(kstr4_kphq,(npw_k))
- + - + -
+ ]
899 [ - + + - : 224 : ABI_MALLOC(kstr5_kphq,(npw_k))
- + - + -
+ ]
900 [ - + + - : 224 : ABI_MALLOC(kstr6_kphq,(npw_k))
- + - + -
+ ]
901 : : end if
902 : : end if
903 : :
904 [ - + + - : 133830 : ABI_MALLOC(kg_k,(3,mpw))
+ - + - -
+ - + -
+ ]
905 : : !$OMP PARALLEL DO
906 [ + + ]: 10162915 : do ipw=1,npw_k
907 [ + + ]: 40517830 : kg_k(:,ipw)=kg(:,ipw+ikg)
908 : : end do
909 : :
910 [ - + + - : 178440 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
+ + + - +
- + + - +
- + - + ]
911 [ + + - + : 44610 : if (stress_needed==1 .or. (use_gbt/=0 .and. psps%useylm==1)) then
- - ]
912 [ - + + - : 177524 : ABI_MALLOC(ylmgr_k,(npw_k,3,mpsang*mpsang*psps%useylm))
+ + + - +
- + + - +
- + - + ]
913 : : else
914 [ - + - + ]: 229 : ABI_MALLOC(ylmgr_k,(0,0,0))
915 : : end if
916 [ + + ]: 44610 : if (psps%useylm==1) then
917 [ + + ]: 9511 : if (use_gbt/=0) then
918 : 8 : call initylmg_k(npw_k,mpsang,stress_needed,rprimd,gs_hamk%gprimd,kpoint,kg_k,ylm_k,ylmgr_k)
919 : : else
920 : : !$OMP PARALLEL DO COLLAPSE(2)
921 [ + + ]: 62835 : do ilm=1,mpsang*mpsang
922 [ + + ]: 13194080 : do ipw=1,npw_k
923 : 13184577 : ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
924 : : end do
925 : : end do
926 [ + + ]: 9503 : if (stress_needed==1) then
927 : : !$OMP PARALLEL DO COLLAPSE(2)
928 [ + + ]: 61540 : do ilm=1,mpsang*mpsang
929 [ + + ]: 218284 : do ii=1,3
930 [ + + ]: 39434310 : do ipw=1,npw_k
931 : 39382062 : ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm)
932 : : end do
933 : : end do
934 : : end do
935 : : end if
936 : : end if
937 : : end if
938 [ + + + + ]: 44610 : if (use_gbt/=0 .and. psps%useylm==1) then
939 [ - + + - : 32 : ABI_MALLOC(ylm_kphq,(npw_k,mpsang*mpsang))
- + + - +
- + - - +
- + - + ]
940 [ - + + - : 32 : ABI_MALLOC(ylmgr_kphq,(npw_k,3,mpsang*mpsang))
- + + - +
- + - - +
- + - + ]
941 : 8 : call initylmg_k(npw_k,mpsang,stress_needed,rprimd,gs_hamk%gprimd,kphq,kg_k,ylm_kphq,ylmgr_kphq)
942 : : else
943 [ - + - + ]: 44602 : ABI_MALLOC(ylm_kphq,(0,0))
944 [ - + - + ]: 44602 : ABI_MALLOC(ylmgr_kphq,(0,0,0))
945 : : end if
946 : :
947 : :
948 : : ! Prepare kinetic contribution to stress tensor (Warning : the symmetry
949 : : ! has not been broken, like in mkkin.f or kpg3.f . It should be, in order to be coherent).
950 [ + + ]: 44610 : if (stress_needed==1) then
951 [ - + ]: 44381 : ABI_MALLOC(gprimd,(3,3))
952 [ - + + + : 621334 : gprimd=gs_hamk%gprimd
+ + ]
953 : : !$OMP PARALLEL DO PRIVATE(fact_kin,ipw,kgr,kgr_kphq,kgc,kgc_kphq,kin,xx,fsm,dfsm) &
954 : : !$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)
955 [ + + ]: 10141954 : do ipw=1,npw_k
956 : : ! Compute Cartesian coordinates of (k+G)
957 [ + + ]: 40390292 : kgr = kpoint + kg_k(:,ipw)
958 [ + + + + ]: 131268449 : kgc = MATMUL(gprimd, kgr)
959 [ + + ]: 40390292 : kin = htpisq * DOT_PRODUCT(kgc, kgc)
960 : : ! kgc1=gprimd(1,1)*(kpoint(1)+kg_k(1,ipw))+&
961 : : ! & gprimd(1,2)*(kpoint(2)+kg_k(2,ipw))+&
962 : : ! & gprimd(1,3)*(kpoint(3)+kg_k(3,ipw))
963 : : ! kgc2=gprimd(2,1)*(kpoint(1)+kg_k(1,ipw))+&
964 : : ! & gprimd(2,2)*(kpoint(2)+kg_k(2,ipw))+&
965 : : ! & gprimd(2,3)*(kpoint(3)+kg_k(3,ipw))
966 : : ! kgc3=gprimd(3,1)*(kpoint(1)+kg_k(1,ipw))+&
967 : : ! & gprimd(3,2)*(kpoint(2)+kg_k(2,ipw))+&
968 : : ! & gprimd(3,3)*(kpoint(3)+kg_k(3,ipw))
969 : : ! kin=htpisq* ( kgc1**2 + kgc2**2 + kgc3**2 )
970 : 10097573 : fact_kin=1.0_dp
971 [ + + ]: 10097573 : if (kin>ecut-ecutsm) then
972 [ + + ]: 495314 : if (kin>ecut) then
973 : : fact_kin=0.0_dp
974 : : else
975 : : ! See the routine mkkin.f, for the smearing procedure
976 : 204208 : xx=(ecut-kin)*ecutsm_inv
977 : : ! This kinetic cutoff smoothing function and its xx derivatives
978 : : ! were produced with Mathematica and the fortran code has been
979 : : ! numerically checked against Mathematica.
980 : 204208 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
981 : 204208 : dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
982 : : ! d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
983 : : ! & (-144+45*xx))))))*fsm**3
984 : 204208 : fact_kin=fsm+kin*(-ecutsm_inv)*dfsm
985 : : end if
986 : : end if
987 : 10097573 : kstr1(ipw) = fact_kin * kgc(1) * kgc(1)
988 : 10097573 : kstr2(ipw) = fact_kin * kgc(2) * kgc(2)
989 : 10097573 : kstr3(ipw) = fact_kin * kgc(3) * kgc(3)
990 : 10097573 : kstr4(ipw) = fact_kin * kgc(3) * kgc(2)
991 : 10097573 : kstr5(ipw) = fact_kin * kgc(3) * kgc(1)
992 : 10141954 : kstr6(ipw) = fact_kin * kgc(2) * kgc(1)
993 : : ! kstr1(ipw)=fact_kin*kgc1*kgc1
994 : : ! kstr2(ipw)=fact_kin*kgc2*kgc2
995 : : ! kstr3(ipw)=fact_kin*kgc3*kgc3
996 : : ! kstr4(ipw)=fact_kin*kgc3*kgc2
997 : : ! kstr5(ipw)=fact_kin*kgc3*kgc1
998 : : ! kstr6(ipw)=fact_kin*kgc2*kgc1
999 : : end do ! ipw
1000 [ + + ]: 44381 : if (use_gbt/=0) then
1001 [ + + ]: 49796 : do ipw=1,npw_k
1002 [ + + ]: 198736 : kgr_kphq = kphq+kg_k(:,ipw)
1003 [ + + + + ]: 645892 : kgc_kphq = MATMUL(gprimd, kgr_kphq)
1004 [ + + ]: 198736 : kin_kphq = htpisq * DOT_PRODUCT(kgc_kphq, kgc_kphq)
1005 : 49684 : fact_kin=1.0_dp
1006 [ + + ]: 49684 : if (kin_kphq>ecut-ecutsm) then
1007 [ + + ]: 1243 : if (kin_kphq>ecut) then
1008 : : fact_kin=0.0_dp
1009 : : else
1010 : 291 : xx=(ecut-kin_kphq)*ecutsm_inv
1011 : 291 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
1012 : 291 : dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
1013 : 291 : fact_kin=fsm+kin_kphq*(-ecutsm_inv)*dfsm
1014 : : end if
1015 : : end if
1016 : 49684 : kstr1_kphq(ipw)=fact_kin*kgc_kphq(1)*kgc_kphq(1)
1017 : 49684 : kstr2_kphq(ipw)=fact_kin*kgc_kphq(2)*kgc_kphq(2)
1018 : 49684 : kstr3_kphq(ipw)=fact_kin*kgc_kphq(3)*kgc_kphq(3)
1019 : 49684 : kstr4_kphq(ipw)=fact_kin*kgc_kphq(3)*kgc_kphq(2)
1020 : 49684 : kstr5_kphq(ipw)=fact_kin*kgc_kphq(3)*kgc_kphq(1)
1021 : 49796 : kstr6_kphq(ipw)=fact_kin*kgc_kphq(2)*kgc_kphq(1)
1022 : : end do ! ipw
1023 : : end if ! GBT
1024 [ - + ]: 44381 : ABI_FREE(gprimd)
1025 : : end if
1026 : :
1027 : :
1028 : : ! Compute (k+G) vectors (only if useylm=1)
1029 [ + + ]: 44610 : if (usexg/=1) then
1030 : 43450 : nkpg=3*nloalg(3)
1031 [ - + + - : 173800 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ + + - +
- + + - +
- + - + ]
1032 [ + + ]: 43450 : if (nkpg>0) call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1033 [ + + ]: 43450 : if (use_gbt/=0) then
1034 [ - + + - : 336 : ABI_MALLOC(kpg_kphq,(npw_k,nkpg))
+ - - - -
+ - + -
+ ]
1035 [ - + ]: 112 : if (nkpg>0) call mkkpg(kg_k,kpg_kphq,kphq,nkpg,npw_k)
1036 : : end if
1037 : : else ! cprj_in_memory = 1
1038 : 1160 : nkpg=3
1039 [ - + + - : 4640 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ - + - +
- - + - +
- + ]
1040 : 1160 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1041 [ - + ]: 1160 : if (use_gbt/=0) then
1042 [ # # # # : 0 : ABI_MALLOC(kpg_kphq,(npw_k,nkpg))
# # # # #
# # # ]
1043 : 0 : call mkkpg(kg_k,kpg_kphq,kphq,nkpg,npw_k)
1044 : : end if
1045 : : end if
1046 : :
1047 : : ! Compute nonlocal form factors ffnl at all (k+G)
1048 : 44610 : ider=0;idir=0;dimffnl=1
1049 [ + + ]: 44610 : if (stress_needed==1) then
1050 : 44381 : ider=1;dimffnl=2+2*psps%useylm
1051 : : end if
1052 [ - + + - : 267660 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1053 : : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1054 : : & ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1055 : 44610 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
1056 [ + + ]: 44610 : if (use_gbt/=0) then
1057 [ - + + - : 672 : ABI_MALLOC(ffnl_kphq,(npw_k,dimffnl,psps%lmnmax,ntypat)) ! 0
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
1058 : : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_kphq,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1059 : : & ider,idir,psps%indlmn,kg_k,kpg_kphq,kphq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1060 : 112 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_kphq,ylmgr_kphq)
1061 : : end if
1062 : : #ifdef HAVE_OPENMP_OFFLOAD
1063 : : !$OMP TARGET ENTER DATA MAP(to:ffnl) IF(gpu_option==ABI_GPU_OPENMP)
1064 : : #endif
1065 [ + + + + : 44610 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
+ + ]
1066 : 54 : ider_str=1; dimffnl_str=7;idir_str=-7
1067 [ - + + - : 324 : ABI_MALLOC(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,ntypat))
+ - - + +
- - + + -
+ - + - -
+ - + -
+ ]
1068 : : call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1069 : : & ider_str,idir_str,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1070 : 54 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
1071 : : #ifdef HAVE_OPENMP_OFFLOAD
1072 : : !$OMP TARGET ENTER DATA MAP(to:ffnl_str) IF(gpu_option==ABI_GPU_OPENMP)
1073 : : #endif
1074 : : end if
1075 : :
1076 : : ! Load k-dependent part in the Hamiltonian datastructure
1077 : : ! - Compute 3D phase factors
1078 : : ! - Prepare various tabs in case of band-FFT parallelism
1079 : : ! - Load k-dependent quantities in the Hamiltonian
1080 [ - + + - : 178440 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
- + + - +
- + - - +
- + - + ]
1081 : : call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
1082 : 44610 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.)
1083 [ + + ]: 44610 : if (use_gbt/=0) then
1084 [ - + + - : 448 : ABI_MALLOC(ph3d_kphq,(2,npw_k,gs_hamk%matblk))
- + + - +
- + - - +
- + - + ]
1085 : : call gs_hamk%load_kprime(kpt_kp=kphq,&
1086 : 112 : & kpg_kp=kpg_kphq,ffnl_kp=ffnl_kphq,ph3d_kp=ph3d_kphq,compute_gbound=compute_gbound,compute_ph3d=.true.)
1087 : : end if
1088 : :
1089 : : ! Load band-FFT tabs (transposed k-dependent arrays)
1090 [ + + + + ]: 44610 : if (mpi_enreg%paral_kgb==1.and.usexg/=1) then
1091 : 2950 : call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1092 : 2950 : call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg)
1093 : : call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
1094 : : & gbound_k =my_bandfft_kpt%gbound, &
1095 : : & kg_k =my_bandfft_kpt%kg_k_gather, &
1096 : : & kpg_k =my_bandfft_kpt%kpg_k_gather, &
1097 : : & ffnl_k =my_bandfft_kpt%ffnl_gather, &
1098 : 2950 : & ph3d_k =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound)
1099 [ - + ]: 2950 : if (use_gbt/=0) then
1100 : : call gs_hamk%load_kprime(npw_fft_kp=my_bandfft_kpt%ndatarecv,&
1101 : : & gbound_kp=my_bandfft_kpt%gbound,kg_kp=my_bandfft_kpt%kg_k_gather,&
1102 : : & kpg_kp=my_bandfft_kpt%kpg_kp_gather,ffnl_kp=my_bandfft_kpt%ffnl_kp_gather,&
1103 : 0 : & ph3d_kp=my_bandfft_kpt%ph3d_kp_gather)
1104 : : end if
1105 : : end if
1106 : :
1107 : : ! If OpenMP GPU, load "hamiltonian" on GPU device
1108 [ - + ]: 44610 : if (gpu_option == ABI_GPU_OPENMP) then
1109 [ # # ]: 0 : if(mpi_enreg%paral_kgb==0) then
1110 : 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k)
1111 [ # # ]: 0 : else if(gs_hamk%istwf_k==1) then
1112 : 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)
1113 : : else
1114 : 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)
1115 : : end if
1116 : : end if
1117 : :
1118 : : ! Setup gemm_nonlop
1119 [ + + ]: 44610 : if (gemm_nonlop_use_gemm) then
1120 : : call set_gemm_nonlop_ikpt(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%istwf_k,gs_hamk%indlmn,&
1121 : 108 : & gs_hamk%ntypat,gs_hamk%nattyp,gs_hamk%gpu_option)
1122 : 108 : call reset_gemm_nonlop()
1123 : :
1124 [ - + ]: 108 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1125 : 0 : gemm_nonlop_block_size = gpu_nl_splitsize
1126 : : call get_gemm_nonlop_ompgpu_blocksize(my_ikpt,gs_hamk,mpi_enreg%bandpp,nband_k,&
1127 : : & nspinor,1,mpi_enreg%paral_kgb,mpi_enreg%nproc_band,&
1128 : : & optfor,stress_needed,-1,gs_hamk%gpu_option,(gpu_nl_distrib/=0),&
1129 : : & gemm_nonlop_block_size,nblk_gemm_nonlop,gs_hamk%nfft_blocks,&
1130 : 0 : & warn_on_fail=.true.,disable_output=(ikpt>1))
1131 [ # # # # ]: 0 : gemm_nonlop_is_distributed = (gpu_nl_distrib/=0 .and. nblk_gemm_nonlop > 0)
1132 [ # # ]: 0 : if(nblk_gemm_nonlop==-1) then
1133 : 0 : gs_hamk%gpu_option=ABI_GPU_DISABLED
1134 : 0 : ABI_WARNING("GPU has been disabled for forces and stress computation due to memory constraints.")
1135 : : end if
1136 : : end if
1137 : : end if
1138 : :
1139 [ + + ]: 44610 : if (usexg==1) then
1140 : 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.)
1141 : : end if
1142 : :
1143 : : ! Loop over (blocks of) bands; accumulate forces and/or stresses
1144 : : ! The following is now wrong. In sequential, nblockbd=nband_k/bandpp
1145 : : ! blocksize= bandpp (JB 2016/04/16)
1146 : : ! Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
1147 [ - + + - : 133830 : ABI_MALLOC(lambda,(blocksize))
+ - + - -
+ - + -
+ ]
1148 [ - + + - : 89220 : ABI_MALLOC(occblock,(blocksize))
- + - + ]
1149 [ - + + - : 89220 : ABI_MALLOC(weight,(blocksize))
- + - + ]
1150 [ - + + - : 133830 : ABI_MALLOC(enlout,(nnlout*blocksize))
+ - + - -
+ - + -
+ ]
1151 [ + + ]: 44610 : if (usexg==1) then
1152 [ + - - + : 4640 : ABI_MALLOC(enlout_2d,(3*natom,blocksize*optfor))
+ - + + +
- + - + +
- + - + -
+ ]
1153 [ + + + - : 3480 : ABI_MALLOC(enlout_2d_stress,(6,blocksize*stress_needed))
+ - + + -
+ - + -
+ ]
1154 [ + + + + : 48644 : if (optfor==1) enlout_2d=zero
+ + ]
1155 [ + + + + : 67604 : if (stress_needed==1) enlout_2d_stress=zero
+ + ]
1156 : : end if
1157 [ + + + + : 1663484 : occblock=zero;weight=zero;enlout(:)=zero
+ + ]
1158 [ + + ]: 44610 : if (usefock_loc) then
1159 [ + + ]: 536 : if (fockcommon%optstr) then
1160 [ - + + - : 1599 : ABI_MALLOC(fockcommon%stress_ikpt,(6,nband_k))
+ - + - -
+ - + -
+ ]
1161 [ + + + + ]: 19041 : fockcommon%stress_ikpt=zero
1162 : : end if
1163 : : end if
1164 [ + + ]: 536 : if ((usefock_loc).and.(psps%usepaw==1)) then
1165 [ + + ]: 57 : if (fockcommon%optfor) then
1166 [ + + + + : 1528 : fockcommon%forces_ikpt=zero
+ + ]
1167 : : end if
1168 : : end if
1169 : :
1170 : 44610 : call timab(922,2,tsec)
1171 : :
1172 [ + + ]: 44610 : if (usexg==1) then
1173 : 1160 : ncols_cprj = bandpp*my_nspinor
1174 : 1160 : call xg_init(cprj_xgx0,xg_nonlop%space_cprj,xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
1175 : 1160 : call xg_init(cprj_work,xg_nonlop%space_cprj,xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
1176 : : end if
1177 : :
1178 [ + + ]: 278536 : do iblock=1,nblockbd
1179 : :
1180 : 233926 : iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
1181 : 233926 : iband_cprj=(iblock-1)*bandpp+1
1182 [ + - ]: 233926 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
1183 : :
1184 : : ! Select occupied bandsddk
1185 [ + + ]: 513606 : occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
1186 [ + - - + : 747532 : if( abs(maxval(occblock))>=tol8 ) then
- - + + +
+ ]
1187 : 159546 : call timab(923,1,tsec)
1188 [ + + ]: 356798 : weight(:)=wtk(ikpt)*occblock(:)
1189 : :
1190 : : ! gs_hamk%ffnl_k is changed in fock_getghc, so that it is necessary to restore it when stresses are to be calculated.
1191 [ + + + + : 159546 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
+ + ]
1192 : 202 : call gs_hamk%load_k(ffnl_k=ffnl)
1193 : : end if
1194 : :
1195 : : ! Load contribution from n,k
1196 : : cwavef(:,1:npw_k*my_nspinor*blocksize)=&
1197 [ + + + + ]: 211007016 : & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
1198 [ + + + + ]: 159546 : if (psps%usepaw==1.and.usecprj_local==1) then
1199 : : call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,&
1200 : : & mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,&
1201 : 8441 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1202 : : end if
1203 : :
1204 : 159546 : call timab(923,2,tsec)
1205 : 159546 : call timab(924,-1,tsec)
1206 : :
1207 [ + + ]: 356798 : lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
1208 : : ABI_NVTX_START_RANGE(NVTX_FORSTR_NONLOP)
1209 [ + + + + ]: 159546 : if (mpi_enreg%paral_kgb/=1.and.usexg/=1) then
1210 [ + + ]: 150664 : if (use_gbt == 0) then
1211 : : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
1212 : 148411 : & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
1213 : : else
1214 : 2253 : gs_hamk%nspinor=1
1215 [ - + + - : 6759 : ABI_MALLOC(cwavef_spin, (2, npw_k*blocksize))
+ - + - -
+ - + -
+ ]
1216 [ - + + - : 6759 : ABI_MALLOC(enlout_spin,(nnlout*blocksize))
+ - + - -
+ - + -
+ ]
1217 [ + + ]: 25068 : enlout_spin(:) = zero
1218 : 2253 : gs_hamk%ispin_gbt = 1
1219 : 2253 : call cg_copy_spin(1,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1220 : : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
1221 : 2253 : & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef_spin,cwavef_spin,select_k=K_H_K)
1222 : :
1223 : 2253 : gs_hamk%ispin_gbt = 2
1224 : 2253 : call cg_copy_spin(2,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1225 : : call nonlop(choice,cpopt,cwaveprj,enlout_spin,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
1226 : 2253 : & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef_spin,cwavef_spin,select_k=KPRIME_H_KPRIME)
1227 : :
1228 [ + + ]: 25068 : enlout(1:nnlout*blocksize) = enlout(1:nnlout*blocksize) + enlout_spin(1:nnlout*blocksize)
1229 : 2253 : gs_hamk%nspinor = 2
1230 [ - + ]: 2253 : ABI_FREE(cwavef_spin)
1231 [ - + ]: 2253 : ABI_FREE(enlout_spin)
1232 : : end if ! GBT
1233 [ + + + + ]: 150664 : if((stress_needed==1).and.(usevxctau==1)) then
1234 : 66 : cwavef_tr=>cwavef
1235 : : endif
1236 [ + + ]: 8882 : else if (usexg/=1) then ! paral_kgb = 1
1237 : : ! Here we MUST pass gpu_option=ABI_GPU_DISABLED, as cwavef is a host buffer.
1238 : : ! For GBT, transpose and evaluate the two diagonal spin blocks separately;
1239 : : ! the default prep_nonlop operation is <k'|Vnl|k> and is not the GBT
1240 : : ! force/stress operator.
1241 [ + - ]: 7534 : if (use_gbt == 0) then
1242 [ + - + + : 7534 : if((stress_needed==1).and.(usevxctau==1).and.mpi_enreg%nproc_band>1) then
+ + ]
1243 [ - + + - : 192 : ABI_MALLOC(cwavef_tr,(2,my_bandfft_kpt%ndatarecv*my_nspinor*mpi_enreg%bandpp))
+ - + - -
+ - + ]
1244 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,&
1245 : : & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,cwavef,&
1246 : 64 : & already_transposed=.False.,gpu_option=ABI_GPU_DISABLED,cwavef_tr=cwavef_tr)
1247 : : else
1248 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,&
1249 : : & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,cwavef,&
1250 : 7470 : & already_transposed=.False.,gpu_option=ABI_GPU_DISABLED)
1251 [ + - + + ]: 7470 : if((stress_needed==1).and.(usevxctau==1)) cwavef_tr=>cwavef
1252 : : end if
1253 : : else
1254 [ # # ]: 0 : ABI_CHECK(usevxctau == 0, 'GBT band-parallel force/stress with meta-GGA is not implemented')
1255 : 0 : gs_hamk%nspinor=1
1256 [ # # # # : 0 : ABI_MALLOC(cwavef_spin,(2,npw_k*blocksize))
# # # # #
# # # #
# ]
1257 [ # # # # : 0 : ABI_MALLOC(enlout_spin,(nnlout*blocksize))
# # # # #
# # # #
# ]
1258 [ # # ]: 0 : enlout_spin=zero
1259 : :
1260 : 0 : gs_hamk%ispin_gbt=1
1261 : 0 : call cg_copy_spin(1,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1262 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,&
1263 : : & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef_spin,cwavef_spin,&
1264 : 0 : & already_transposed=.False.,gpu_option=ABI_GPU_DISABLED,select_k=K_H_K)
1265 : :
1266 : 0 : gs_hamk%ispin_gbt=2
1267 : 0 : call cg_copy_spin(2,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1268 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout_spin,gs_hamk,idir,lambda,blocksize,&
1269 : : & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef_spin,cwavef_spin,&
1270 : 0 : & already_transposed=.False.,gpu_option=ABI_GPU_DISABLED,select_k=KPRIME_H_KPRIME)
1271 : :
1272 [ # # # # ]: 0 : enlout=enlout+enlout_spin
1273 : 0 : gs_hamk%nspinor=2
1274 [ # # ]: 0 : ABI_FREE(cwavef_spin)
1275 [ # # ]: 0 : ABI_FREE(enlout_spin)
1276 : : end if
1277 : : else ! usexg==1
1278 [ + + ]: 1348 : if ( istwf_k > 1 ) then ! Real only
1279 : 529 : space = SPACE_CR
1280 : : else ! complex
1281 : 819 : space = SPACE_C
1282 : : end if
1283 : 1348 : me_g0 = -1
1284 : 1348 : me_g0_fft = -1
1285 [ + + ]: 1348 : if (space==SPACE_CR) then
1286 : 529 : me_g0 = 0
1287 : 529 : me_g0_fft = 0
1288 [ + + ]: 529 : if (istwf_k == 2) then
1289 [ + + ]: 81 : if (mpi_enreg%me_g0 == 1) me_g0 = 1
1290 [ + - ]: 81 : if (mpi_enreg%me_g0_fft == 1) me_g0_fft = 1
1291 : : end if
1292 : : end if
1293 : :
1294 : : call xgBlock_map(xgx0,cwavef,space,npw_k*my_nspinor,blocksize,comm=mpi_enreg%comm_band,me_g0=me_g0,&
1295 : 1348 : & gpu_option=gpu_option)
1296 : :
1297 [ + + + + ]: 1348 : if ((stress_needed==1).and.(usevxctau==1)) then
1298 [ + + ]: 66 : if(mpi_enreg%nproc_band>1) then
1299 : : call xgTransposer_constructor(xgTransposer,xgx0,xgx0_tr,my_nspinor,&
1300 : 64 : STATE_LINALG,TRANS_ALL2ALL,mpi_enreg%comm_spinorfft,mpi_enreg%comm_band,0,0,me_g0_fft)
1301 : 64 : call xgTransposer_transpose(xgTransposer,STATE_COLSROWS)
1302 : 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)
1303 : : else
1304 : 2 : cwavef_tr=>cwavef
1305 : : endif
1306 : : endif
1307 : :
1308 : :
1309 : 1348 : call xgBlock_map_1d(xgeigen,lambda,SPACE_R,blocksize)
1310 : :
1311 [ + + - + ]: 1348 : if (psps%usepaw==1.and.usecprj_local==1) then
1312 : 0 : call xg_cprj_copy(cwaveprj,cprj_xgx0%self,xg_nonlop,CPRJ_TO_XG)
1313 : : else
1314 : 1348 : call xg_nonlop_getcprj(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self)
1315 : : end if
1316 : :
1317 [ + + ]: 1348 : if (optfor==1) call xgBlock_map(xgforces,enlout_2d,SPACE_R,3*natom,blocksize)
1318 [ + + ]: 1348 : if (stress_needed==1) call xgBlock_map(xgstress,enlout_2d_stress,SPACE_R,6,blocksize)
1319 : :
1320 [ + + + + ]: 1348 : if (optfor==1.and.stress_needed==0) then
1321 : : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,&
1322 : 208 : forces=xgforces)
1323 : : end if
1324 [ + + + - ]: 1348 : if (optfor==0.and.stress_needed==1) then
1325 : : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,&
1326 : 800 : stress=xgstress,gprimd=gs_hamk%gprimd)
1327 : : end if
1328 [ + + + + ]: 1348 : if (optfor==1.and.stress_needed==1) then
1329 : : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,&
1330 : 340 : forces=xgforces,stress=xgstress,gprimd=gs_hamk%gprimd)
1331 : : end if
1332 : :
1333 : : end if ! end usexg==1
1334 : : ABI_NVTX_END_RANGE()
1335 : :
1336 [ + + + + : 159546 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
+ + ]
1337 : 202 : call gs_hamk%load_k(ffnl_k=ffnl_str)
1338 : : end if
1339 : :
1340 : : ! Accumulate non-local contributions from n,k
1341 [ + + ]: 159546 : if (usexg/=1) then
1342 [ + + ]: 158198 : if (optfor==1) then
1343 [ + + ]: 242677 : do iblocksize=1,blocksize
1344 : 129599 : ibs=nnlout*(iblocksize-1)
1345 [ + + ]: 2100904 : grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout(ibs+1+ishift:ibs+3*natom+ishift)
1346 : : end do
1347 : : end if
1348 [ + + ]: 158198 : if (stress_needed==1) then
1349 [ + + ]: 342358 : do iblocksize=1,blocksize
1350 : 184264 : ibs=nnlout*(iblocksize-1)
1351 [ + + ]: 1447942 : npsstr(1:6)=npsstr(1:6) + weight(iblocksize)*enlout(ibs+1:ibs+6)
1352 : : end do
1353 : : end if
1354 : : else
1355 [ + + ]: 1348 : if (optfor==1) then
1356 [ + + ]: 7128 : do iblocksize=1,blocksize
1357 [ + + ]: 49488 : grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout_2d(1:3*natom,iblocksize)
1358 : : end do
1359 : : end if
1360 [ + + ]: 1348 : if (stress_needed==1) then
1361 [ + + ]: 11208 : do iblocksize=1,blocksize
1362 [ + + ]: 71616 : npsstr(1:6) = npsstr(1:6) + weight(iblocksize)*enlout_2d_stress(1:6,iblocksize)
1363 : : end do
1364 : : end if
1365 : : end if
1366 : :
1367 : 159546 : call timab(924,2,tsec)
1368 : :
1369 : : #if defined HAVE_GPU && defined HAVE_YAKL
1370 : : if(gpu_option==ABI_GPU_KOKKOS) then
1371 : : ! the following is done on CPU, so prefetch wave functions from device to host (for efficiency)
1372 : : call gpu_data_prefetch_async(C_LOC(cwavef), INT(2, c_size_t)*npw_k*my_nspinor*blocksize, CPU_DEVICE_ID)
1373 : : call gpu_device_synchronize()
1374 : : end if
1375 : : #endif
1376 : :
1377 : : ! Accumulate stress tensor kinetic contributions
1378 [ + + ]: 159546 : if (stress_needed==1) then
1379 : 159234 : call timab(925,1,tsec)
1380 [ + + ]: 159234 : if (use_gbt == 0) then
1381 [ + + ]: 349060 : do iblocksize=1,blocksize
1382 : : call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1383 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1384 : 192079 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1385 : 192079 : kinstr(1)=kinstr(1)+weight(iblocksize)*ar
1386 : : call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1387 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1388 : 192079 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1389 : 192079 : kinstr(2)=kinstr(2)+weight(iblocksize)*ar
1390 : : call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1391 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1392 : 192079 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1393 : 192079 : kinstr(3)=kinstr(3)+weight(iblocksize)*ar
1394 : : call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1395 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1396 : 192079 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1397 : 192079 : kinstr(4)=kinstr(4)+weight(iblocksize)*ar
1398 : : call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1399 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1400 : 192079 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1401 : 192079 : kinstr(5)=kinstr(5)+weight(iblocksize)*ar
1402 : : call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1403 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1404 : 192079 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1405 : 349060 : kinstr(6)=kinstr(6)+weight(iblocksize)*ar
1406 : : end do
1407 : : else
1408 [ + + ]: 4506 : do iblocksize=1,blocksize
1409 : : call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,1,&
1410 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1411 : 2253 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1412 : : call meanvalue_g(ar2,kstr1_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1413 : : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1414 : 2253 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1415 : 2253 : kinstr(1)=kinstr(1)+weight(iblocksize)*(ar+ar2) ! (1,1)
1416 : : call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,1,&
1417 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1418 : 2253 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1419 : : call meanvalue_g(ar2,kstr2_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1420 : : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1421 : 2253 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1422 : 2253 : kinstr(2)=kinstr(2)+weight(iblocksize)*(ar+ar2) ! (2,2)
1423 : : call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,1,&
1424 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1425 : 2253 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1426 : : call meanvalue_g(ar2,kstr3_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1427 : : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1428 : 2253 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1429 : 2253 : kinstr(3)=kinstr(3)+weight(iblocksize)*(ar+ar2) ! (3,3)
1430 : : call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,1,&
1431 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1432 : 2253 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1433 : : call meanvalue_g(ar2,kstr4_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1434 : : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1435 : 2253 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1436 : 2253 : kinstr(4)=kinstr(4)+weight(iblocksize)*(ar+ar2) ! (3,2)
1437 : : call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,1,&
1438 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1439 : 2253 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1440 : : call meanvalue_g(ar2,kstr5_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1441 : : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1442 : 2253 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1443 : 2253 : kinstr(5)=kinstr(5)+weight(iblocksize)*(ar+ar2) ! (3,1)
1444 : : call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,1,&
1445 : : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),&
1446 : 2253 : & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:(iblocksize-1)*npw_k*my_nspinor+npw_k),0)
1447 : : call meanvalue_g(ar2,kstr6_kphq,0,istwf_k,mpi_enreg,npw_k,1,&
1448 : : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),&
1449 : 2253 : & cwavef(:,1 + (iblocksize-1)*npw_k*my_nspinor+npw_k:),0)
1450 : 4506 : kinstr(6)=kinstr(6)+weight(iblocksize)*(ar+ar2) ! (2,1)`
1451 : : end do
1452 : : end if ! GBT
1453 : 159234 : call timab(925,2,tsec)
1454 : : end if
1455 : :
1456 : : ! Accumulate stress tensor in case meta-GGA using v_tau
1457 [ + + + + ]: 159546 : if ((stress_needed==1).and.(usevxctau==1)) then
1458 [ + + ]: 197 : if(mpi_enreg%paral_kgb==1) then
1459 : 131 : gbound_k => my_bandfft_kpt%gbound
1460 : 131 : kg_k_ => my_bandfft_kpt%kg_k_gather
1461 : 131 : npw_k_=my_bandfft_kpt%ndatarecv
1462 : : else
1463 : 66 : gbound_k => gs_hamk%gbound_k
1464 : 66 : kg_k_ => kg_k
1465 : 66 : npw_k_=npw_k
1466 : : endif
1467 : : call stress_mGGA(mggastr,cwavef_tr,effmass_free,gbound_k,gs_hamk%gprimd,istwf_k, &
1468 : : & kg_k_,kpoint,mgfft,mpi_enreg,my_nspinor,mpi_enreg%bandpp,ngfft,npw_k_,gs_hamk%nvloc, &
1469 : : & 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, &
1470 [ - + - + : 197 : & wtk(ikpt),gpu_option=gpu_option)
- + ]
1471 [ + + + + ]: 197 : if(mpi_enreg%paral_kgb==1.and.mpi_enreg%nproc_band>1) then
1472 [ + + ]: 128 : if(usexg==1) then
1473 : 64 : call xgTransposer_free(xgTransposer)
1474 : : else
1475 [ - + ]: 64 : ABI_FREE(cwavef_tr)
1476 : : endif
1477 : : else
1478 : 69 : nullify(cwavef_tr)
1479 : : endif
1480 : : end if
1481 : :
1482 : : ! Accumulate stress tensor and forces for the Fock part
1483 [ + + ]: 159546 : if (usefock_loc) then
1484 [ + + + - ]: 2052 : if(fockcommon%optstr.or.fockcommon%optfor) then
1485 : 2052 : call timab(926,1,tsec)
1486 [ - + ]: 2052 : if (mpi_enreg%paral_kgb==1) then
1487 : 0 : msg='forsrtnps: Paral_kgb is not yet implemented for fock stresses'
1488 : 0 : ABI_BUG(msg)
1489 : : end if
1490 : 2052 : ndat=mpi_enreg%bandpp
1491 [ + + ]: 2052 : if (gs_hamk%usepaw==0) cwaveprj_idat => cwaveprj
1492 [ - + - + ]: 2052 : ABI_MALLOC(ghc_dum,(0,0))
1493 : 2052 : fockcommon%ieigen=(iblock-1)*blocksize+1
1494 : 2052 : fockcommon%iband=(iblock-1)*blocksize+1
1495 [ + + ]: 2052 : if (gs_hamk%usepaw==1) then
1496 : 214 : cwaveprj_idat => cwaveprj(:,1+blocksize*(my_nspinor-1):blocksize*my_nspinor)
1497 : : end if
1498 : : call fock_getghc(cwavef(:,1+blocksize*npw_k*(my_nspinor-1):blocksize*npw_k*my_nspinor),cwaveprj_idat,&
1499 : 2052 : & ghc_dum,gs_hamk,mpi_enreg,blocksize)
1500 [ + + ]: 4104 : do iblocksize=1,blocksize
1501 [ + + ]: 2052 : if (fockcommon%optstr) then
1502 [ + + ]: 14280 : fockcommon%stress(:)=fockcommon%stress(:)+weight(iblocksize)*fockcommon%stress_ikpt(:,fockcommon%ieigen+iblocksize-1)
1503 : : end if
1504 [ + + ]: 4104 : if (fockcommon%optfor) then
1505 [ + + + + ]: 1224 : fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen+iblocksize-1)
1506 : : end if
1507 : : end do
1508 [ - + ]: 2052 : ABI_FREE(ghc_dum)
1509 : 2052 : call timab(926,2,tsec)
1510 : : end if
1511 : : end if ! usefock_loc
1512 : : end if
1513 [ - + ]: 278536 : if ( gpu_option == ABI_GPU_OPENMP) then
1514 : 0 : call ompgpu_free_hamilt_buffers()
1515 : : end if
1516 : :
1517 : : end do ! End of loop on block of bands
1518 : :
1519 [ + + ]: 44610 : if (usexg==1) then
1520 : 1160 : call xg_free(cprj_xgx0)
1521 : 1160 : call xg_free(cprj_work)
1522 : : end if
1523 : :
1524 : 44610 : call timab(927,1,tsec)
1525 : :
1526 : : ! Restore the bandfft tabs
1527 [ + + ]: 44610 : if (mpi_enreg%paral_kgb==1) then
1528 : 3764 : call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1529 : : end if
1530 : :
1531 : : ! Increment indexes
1532 : 44610 : bdtot_index=bdtot_index+nband_k
1533 [ + - ]: 44610 : if (mkmem/=0) then
1534 : 44610 : ibg=ibg+my_nspinor*nband_cprj_k
1535 : 44610 : icg=icg+npw_k*my_nspinor*nband_k
1536 : 44610 : ikg=ikg+npw_k
1537 : : end if
1538 : :
1539 [ + + ]: 44610 : if (usefock_loc) then
1540 [ + + ]: 536 : if (fockcommon%optstr) then
1541 [ - + ]: 533 : ABI_FREE(fockcommon%stress_ikpt)
1542 : : end if
1543 : : end if
1544 : :
1545 [ + + ]: 44610 : if (psps%usepaw==1) then
1546 : 9147 : call pawcprj_free(cwaveprj)
1547 : : end if
1548 [ + - + + : 47539 : ABI_FREE(cwaveprj)
- + - + -
+ ]
1549 : :
1550 [ + - ]: 44610 : if(gpu_option == ABI_GPU_KOKKOS) then
1551 : : #if defined HAVE_GPU && defined HAVE_YAKL
1552 : : ABI_FREE_MANAGED(cwavef)
1553 : : #endif
1554 : : else
1555 [ - + ]: 44610 : ABI_FREE(cwavef)
1556 : : end if
1557 : :
1558 [ - + ]: 44610 : ABI_FREE(lambda)
1559 [ - + ]: 44610 : ABI_FREE(occblock)
1560 [ - + ]: 44610 : ABI_FREE(weight)
1561 [ - + ]: 44610 : ABI_FREE(enlout)
1562 [ + + ]: 44610 : if (usexg==1) then
1563 [ - + ]: 1160 : ABI_FREE(enlout_2d)
1564 [ - + ]: 1160 : ABI_FREE(enlout_2d_stress)
1565 : : end if
1566 : : #ifdef HAVE_OPENMP_OFFLOAD
1567 : : !$OMP TARGET EXIT DATA MAP(delete:ffnl) IF(gpu_option==ABI_GPU_OPENMP)
1568 : : #endif
1569 [ - + ]: 44610 : ABI_FREE(ffnl)
1570 [ + + ]: 44610 : ABI_SFREE(ffnl_kphq)
1571 [ - + ]: 44610 : ABI_FREE(kg_k)
1572 [ - + ]: 44610 : ABI_FREE(kpg_k)
1573 [ + + ]: 44610 : ABI_SFREE(kpg_kphq)
1574 [ - + ]: 44610 : ABI_FREE(ylm_k)
1575 [ - + ]: 44610 : ABI_FREE(ylmgr_k)
1576 [ - + ]: 44610 : ABI_FREE(ylm_kphq)
1577 [ - + ]: 44610 : ABI_FREE(ylmgr_kphq)
1578 [ - + ]: 44610 : ABI_FREE(ph3d)
1579 [ + + ]: 44610 : ABI_SFREE(ph3d_kphq)
1580 : : ABI_SFREE(ffnl_kphq_gather)
1581 : : ABI_SFREE(kinpw_kphq_gather)
1582 : : ABI_SFREE(kpg_kphq_gather)
1583 : : ABI_SFREE(ph3d_kphq_gather)
1584 [ + + ]: 44610 : if (stress_needed==1) then
1585 [ - + ]: 44381 : ABI_FREE(kstr1)
1586 [ - + ]: 44381 : ABI_FREE(kstr2)
1587 [ - + ]: 44381 : ABI_FREE(kstr3)
1588 [ - + ]: 44381 : ABI_FREE(kstr4)
1589 [ - + ]: 44381 : ABI_FREE(kstr5)
1590 [ - + ]: 44381 : ABI_FREE(kstr6)
1591 [ + + ]: 44381 : if (use_gbt /= 0) then
1592 [ - + ]: 112 : ABI_FREE(kstr1_kphq)
1593 [ - + ]: 112 : ABI_FREE(kstr2_kphq)
1594 [ - + ]: 112 : ABI_FREE(kstr3_kphq)
1595 [ - + ]: 112 : ABI_FREE(kstr4_kphq)
1596 [ - + ]: 112 : ABI_FREE(kstr5_kphq)
1597 [ - + ]: 112 : ABI_FREE(kstr6_kphq)
1598 : : end if
1599 : : end if
1600 [ + + + + : 44610 : if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
+ + ]
1601 : : #ifdef HAVE_OPENMP_OFFLOAD
1602 : : !$OMP TARGET EXIT DATA MAP(delete:ffnl_str) IF(gpu_option==ABI_GPU_OPENMP)
1603 : : #endif
1604 [ - + ]: 54 : ABI_FREE(ffnl_str)
1605 : : end if
1606 : :
1607 : 51539 : call timab(927,2,tsec)
1608 : :
1609 : : end do ! End k point loop
1610 : :
1611 : : end do ! End loop over spins
1612 : :
1613 : 6321 : call timab(928,1,tsec)
1614 : :
1615 [ + + ]: 6321 : if (usexg==1) then
1616 [ + + ]: 166 : if (xg_nonlop%paw) call xg_nonlop_destroy_Dij(xg_nonlop)
1617 : : end if
1618 : :
1619 : : !Stress is equal to dE/d_strain * (1/ucvol)
1620 [ + + ]: 44247 : npsstr(:)=npsstr(:)/gs_hamk%ucvol
1621 : :
1622 : : !Parallel case: accumulate (n,k) contributions
1623 : : if (xmpi_paral==1) then
1624 : : ! Forces
1625 [ + + ]: 6321 : if (optfor==1) then
1626 : 3236 : call timab(65,1,tsec)
1627 : 3236 : call xmpi_sum(grnl,spaceComm,ierr)
1628 : 3236 : call timab(65,2,tsec)
1629 [ + + + + ]: 3236 : if ((usefock_loc).and.(psps%usepaw==1)) then
1630 : 10 : call xmpi_sum(fockcommon%forces,spaceComm,ierr)
1631 : : end if
1632 : : end if
1633 : : ! Stresses
1634 [ + + ]: 6321 : if (stress_needed==1) then
1635 : 6290 : call timab(65,1,tsec)
1636 : 6290 : call xmpi_sum(kinstr,spaceComm,ierr)
1637 : 6290 : call xmpi_sum(npsstr,spaceComm,ierr)
1638 [ + + ]: 6290 : if (usevxctau>0) then
1639 : 37 : call xmpi_sum(mggastr,spaceComm,ierr)
1640 : : end if
1641 [ + + ]: 6290 : if (usefock_loc) then
1642 [ + - ]: 121 : if (fockcommon%optstr) then
1643 : 121 : call xmpi_sum(fockcommon%stress,spaceComm,ierr)
1644 : : end if
1645 : : end if
1646 : 6290 : call timab(65,2,tsec)
1647 : : end if
1648 : : end if
1649 : :
1650 : : !Do final normalizations and symmetrizations of stress tensor contributions
1651 [ + + ]: 6321 : if (stress_needed==1) then
1652 : 6290 : renorm_factor=-(two_pi**2)/effmass_free/gs_hamk%ucvol
1653 [ + + ]: 44030 : kinstr(:)=kinstr(:)*renorm_factor
1654 [ + + ]: 6290 : if (nsym>1) then
1655 : 4860 : call stresssym(gs_hamk%gprimd,nsym,kinstr,symrec)
1656 : 4860 : call stresssym(gs_hamk%gprimd,nsym,npsstr,symrec)
1657 [ + + ]: 4860 : if (usevxctau>0) then
1658 : 37 : call stresssym(gs_hamk%gprimd,nsym,mggastr,symrec)
1659 : : end if
1660 [ + + ]: 4860 : if (usefock_loc) then
1661 [ + - ]: 121 : if (fockcommon%optstr) then
1662 : 121 : call stresssym(gs_hamk%gprimd,nsym,fockcommon%stress,symrec)
1663 : : end if
1664 : : end if
1665 : : end if
1666 : : end if
1667 : :
1668 : : !Need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted)
1669 [ + + + + ]: 6321 : if (psps%usepaw==1.and.usecprj_local==1) then
1670 : 106 : call pawcprj_reorder(cprj,gs_hamk%atindx1)
1671 : : end if
1672 : :
1673 : : !Deallocate temporary space
1674 [ + + ]: 6321 : if (usevxctau>0) then
1675 [ - + ]: 37 : ABI_FREE(vxctaulocal)
1676 : : end if
1677 : 6321 : call gs_hamk%free()
1678 [ + + ]: 6321 : if (usefock_loc) then
1679 : 122 : fockcommon%use_ACE=use_ACE_old
1680 : : end if
1681 : :
1682 : 6321 : call timab(928,2,tsec) ; call timab(920,-2,tsec)
1683 : : ABI_NVTX_END_RANGE()
1684 : :
1685 [ - + - + : 12642 : end subroutine forstrnps
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - + -
+ ]
1686 : : !!***
1687 : :
1688 : : !----------------------------------------------------------------------
1689 : :
1690 : : !!****f* ABINIT/nres2vres
1691 : : !!
1692 : : !! NAME
1693 : : !! nres2vres
1694 : : !!
1695 : : !! FUNCTION
1696 : : !! Convert a density residual into a potential residual
1697 : : !! using a first order formula:
1698 : : !! V^res(r)=dV/dn.n^res(r)
1699 : : !! =V_hartree(n^res)(r) + Kxc.n^res(r)
1700 : : !! INPUTS
1701 : : !! dtset <type(dataset_type)>=all input variables in this dataset
1702 : : !! | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver
1703 : : !! | natom= number of atoms in cell
1704 : : !! | nspden=number of spin-density components
1705 : : !! | ntypat=number of atom types
1706 : : !! | typat(natom)=type (integer) for each atom
1707 : : !! gsqcut=cutoff value on G**2 for sphere inside fft box
1708 : : !! izero=if 1, unbalanced components of Vhartree(g) are set to zero
1709 : : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
1710 : : !! mpi_enreg=information about MPI parallelization
1711 : : !! my_natom=number of atoms treated by current processor
1712 : : !! nfft=(effective) number of FFT grid points (for this processor)
1713 : : !! ngfft(18)=contain all needed information about 3D FFT
1714 : : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
1715 : : !! nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
1716 : : !! nresid(nfft,nspden)= the input density residual
1717 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
1718 : : !! optnc=option for non-collinear magnetism (nspden=4):
1719 : : !! 1: the whole 2x2 Vres matrix is computed
1720 : : !! 2: only Vres^{11} and Vres^{22} are computed
1721 : : !! optxc=0 if LDA part of XC kernel has only to be taken into account (even for GGA)
1722 : : !! 1 if XC kernel has to be fully taken into
1723 : : !! -1 if XC kernel does not have to be taken into account
1724 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
1725 : : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
1726 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
1727 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1728 : : !! rhor(nfft,nspden)=electron density in real space
1729 : : !! (used only if Kxc was not computed before)
1730 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
1731 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
1732 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
1733 : : !! xcctau3d(n3xccc)=3D core electron kinetic energy density for XC core correction (bohr^-3)
1734 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
1735 : : !!
1736 : : !! === optional inputs ===
1737 : : !! vxc(cplex*nfft,nspden)=XC GS potential
1738 : : !!
1739 : : !! OUTPUT
1740 : : !! vresid(nfft,nspden)= the output potential residual
1741 : : !!
1742 : : !! SOURCE
1743 : :
1744 : 6729 : subroutine nres2vres(dtset,gsqcut,izero,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
1745 : 6729 : & nkxc,nresid,n3xccc,optnc,optxc,pawang,pawfgrtab,pawrhoij,pawtab,&
1746 : 6729 : & rhor,rprimd,usepaw,vresid,xccc3d,xred,&
1747 : 715 : & vxc, xcctau3d)
1748 : :
1749 : : !Arguments ------------------------------------
1750 : : !scalars
1751 : : integer,intent(in) :: izero,my_natom,n3xccc,nfft,nkxc,optnc,optxc,usepaw
1752 : : real(dp),intent(in) :: gsqcut
1753 : : type(MPI_type),intent(in) :: mpi_enreg
1754 : : type(dataset_type),intent(in) :: dtset
1755 : : type(pawang_type),intent(in) :: pawang
1756 : : !arrays
1757 : : integer,intent(in) :: ngfft(18)
1758 : : real(dp),intent(in) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden)
1759 : : real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3),xccc3d(n3xccc),xred(3,dtset%natom)
1760 : : real(dp),intent(in),optional :: xcctau3d(n3xccc*dtset%usekden)
1761 : : real(dp),intent(inout) :: nhat(nfft,dtset%nspden*usepaw)
1762 : : real(dp),intent(out) :: vresid(nfft,dtset%nspden)
1763 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*usepaw)
1764 : : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*usepaw)
1765 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw)
1766 : : real(dp),intent(in) :: vxc(nfft,dtset%nspden) !FR TODO:cplex?
1767 : :
1768 : : !Local variables-------------------------------
1769 : : !scalars
1770 : : integer :: cplex,ider,idir,ipert,ispden,nhatgrdim,nkxc_cur,option,me,nproc,comm,usexcnhat
1771 : : logical :: has_nkxc_gga,non_magnetic_xc
1772 : : real(dp) :: dum,energy,m_norm_min,ucvol,vxcavg
1773 : : character(len=500) :: message
1774 : : type(xcdata_type) :: xcdata
1775 : : !arrays
1776 : : integer :: nk3xc
1777 : : real(dp) :: gmet(3,3),gprimd(3,3),qq(3),rmet(3,3)
1778 : 6729 : real(dp),allocatable :: dummy(:),kxc_cur(:,:),nhatgr(:,:,:)
1779 : 6729 : real(dp),allocatable :: nresg(:,:),rhor0(:,:),vhres(:)
1780 : : ! *************************************************************************
1781 : :
1782 : : !Compatibility tests:
1783 : 6729 : has_nkxc_gga=(nkxc==7.or.nkxc==19)
1784 : :
1785 [ - + ]: 6729 : if(optxc<-1.or.optxc>1)then
1786 : 0 : write(message,'(a,i0)')' Wrong value for optxc ',optxc
1787 : 0 : ABI_BUG(message)
1788 : : end if
1789 : :
1790 [ + - + + : 6729 : if((optnc/=1.and.optnc/=2).or.(dtset%nspden/=4.and.optnc/=1))then
- + ]
1791 : 0 : write(message,'(a,i0)')' Wrong value for optnc ',optnc
1792 : 0 : ABI_BUG(message)
1793 : : end if
1794 : :
1795 [ - + - - ]: 6729 : if(dtset%icoulomb==1.and.optxc/=-1)then
1796 : 0 : write(message,'(a)')' This routine is not compatible with icoulomb==1 and optxc/=-1 !'
1797 : 0 : ABI_BUG(message)
1798 : : end if
1799 : :
1800 [ + + + + : 6729 : if(dtset%nspden==4.and.dtset%xclevel==2.and.optxc==1.and.(.not.has_nkxc_gga))then
- + - - ]
1801 : 0 : ABI_ERROR(' Wrong values for optxc and nkxc !')
1802 : : end if
1803 : :
1804 : 6729 : qq=zero
1805 : 6729 : nkxc_cur=0
1806 : 6729 : m_norm_min=EPSILON(0.0_dp)**2
1807 [ + + + + ]: 17374 : usexcnhat=0;if (usepaw==1) usexcnhat=maxval(pawtab(1:dtset%ntypat)%usexcnhat)
1808 [ + + + + ]: 6729 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
1809 [ + + + - ]: 6729 : if (dtset%xclevel==1.or.optxc==0) nkxc_cur= 2*min(dtset%nspden,2)-1 ! LDA: nkxc=1,3
1810 [ + + - + ]: 6729 : if (dtset%xclevel==2.and.optxc==1)nkxc_cur=12*min(dtset%nspden,2)-5 ! GGA: nkxc=7,19
1811 [ - + + - : 20187 : ABI_MALLOC(vhres,(nfft))
+ - + - -
+ - + ]
1812 : :
1813 : : !Compute different geometric tensor, as well as ucvol, from rprimd
1814 : 6729 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1815 : :
1816 : : !Compute density residual in reciprocal space
1817 [ + - ]: 6729 : if (dtset%icoulomb==0) then
1818 [ - + + - : 20187 : ABI_MALLOC(nresg,(2,nfft))
+ - + - -
+ - + ]
1819 [ - + + - : 13458 : ABI_MALLOC(dummy,(nfft))
- + + - ]
1820 [ + + ]: 335772490 : dummy(:)=nresid(:,1)
1821 : 6729 : call fourdp(1,nresg,dummy,-1,mpi_enreg,nfft,1,ngfft,0)
1822 [ - + ]: 6729 : ABI_FREE(dummy)
1823 : : end if
1824 : :
1825 : : !For GGA, has to recompute gradients of nhat
1826 : 6729 : nhatgrdim=0
1827 [ + + - + : 6729 : if ((nkxc==nkxc_cur.and.has_nkxc_gga).or.(optxc==-1.and.has_nkxc_gga).or.&
+ - + - +
+ ]
1828 : : & (optxc/=-1.and.nkxc/=nkxc_cur)) then
1829 [ - + - - : 1 : if (usepaw==1.and.dtset%xclevel==2.and.usexcnhat>0.and.dtset%pawnhatxc>0) then
- - - - ]
1830 : 0 : nhatgrdim=1
1831 [ # # # # : 0 : ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3))
# # # # #
# # # # #
# # # # ]
1832 : 0 : ider=1;cplex=1;ipert=0;idir=0
1833 : : call pawmknhat(dum,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
1834 : : & nfft,ngfft,nhatgrdim,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,&
1835 : : & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qq,rprimd,ucvol,dtset%usewvl,xred,&
1836 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1837 : : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1838 : 0 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1839 : : else
1840 [ - + ]: 1 : ABI_MALLOC(nhatgr,(0,0,0))
1841 : : end if
1842 : : else
1843 [ - + ]: 6728 : ABI_MALLOC(nhatgr,(0,0,0))
1844 : : end if
1845 : :
1846 [ - + ]: 6729 : ABI_MALLOC(dummy,(0))
1847 : : !First case: Kxc has already been computed
1848 : : !-----------------------------------------
1849 [ + + - + ]: 6729 : if (nkxc==nkxc_cur.or.optxc==-1) then
1850 : :
1851 : : ! Compute VH(n^res)(r)
1852 [ + - ]: 6728 : if (dtset%icoulomb == 0) then
1853 : : call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,&
1854 : 6728 : &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres)
1855 : : else
1856 : 0 : comm=mpi_enreg%comm_cell
1857 : 0 : nproc=xmpi_comm_size(comm)
1858 : 0 : me=xmpi_comm_rank(comm)
1859 : : call psolver_hartree(energy, (/ rprimd(1,1) / dtset%ngfft(1), &
1860 : : & rprimd(2,2) / dtset%ngfft(2), rprimd(3,3) / dtset%ngfft(3) /), dtset%icoulomb, &
1861 : : & me, comm, dtset%nfft, dtset%ngfft(1:3), nproc, dtset%nscforder, dtset%nspden, &
1862 [ # # ]: 0 : & nresid(:,1), vhres, dtset%usewvl)
1863 : : end if
1864 : :
1865 : : ! Compute Kxc(r).n^res(r)
1866 [ - + ]: 6728 : if (optxc/=-1) then
1867 : :
1868 : : ! Collinear magnetism or non-polarized
1869 [ + + ]: 6728 : if (dtset%nspden/=4) then
1870 : : !Note: imposing usexcnhat=1 avoid nhat to be substracted
1871 : : call dfpt_mkvxc(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,&
1872 : 6639 : & nkxc,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy)
1873 : : else
1874 : : !FR call routine for Non-collinear magnetism
1875 [ - + + - : 356 : ABI_MALLOC(rhor0,(nfft,dtset%nspden))
- + + - +
- + - - +
+ - ]
1876 [ + + + + ]: 3635833 : rhor0(:,:)=rhor(:,:)-nresid(:,:)
1877 : : !Note: imposing usexcnhat=1 avoid nhat to be substracted
1878 : : call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,&
1879 : 89 : & nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d)
1880 [ - + ]: 89 : ABI_FREE(rhor0)
1881 : : end if
1882 : :
1883 : : else
1884 [ - - - - ]: 1 : vresid=zero
1885 : : end if
1886 : :
1887 : : end if
1888 : :
1889 : : !2nd case: Kxc has to be computed
1890 : : !--------------------------------
1891 [ + + + - ]: 6729 : if (nkxc/=nkxc_cur.and.optxc/=-1) then
1892 : :
1893 : : ! Has to use the "initial" density to compute Kxc
1894 [ - + + - : 4 : ABI_MALLOC(rhor0,(nfft,dtset%nspden))
- + + - +
- + - - +
+ - ]
1895 [ + + + + ]: 6753 : rhor0(:,:)=rhor(:,:)-nresid(:,:)
1896 : :
1897 : : ! Compute VH(n^res) and XC kernel (Kxc) together
1898 [ - + + - : 4 : ABI_MALLOC(kxc_cur,(nfft,nkxc_cur))
- + + - +
- + - - +
- + ]
1899 : :
1900 [ + - + - ]: 1 : option=2;if (dtset%xclevel==2.and.optxc==0) option=12
1901 : :
1902 : : call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,&
1903 : 1 : &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres)
1904 : 1 : call xcdata_init(xcdata,dtset=dtset)
1905 : :
1906 : : ! To be adjusted for the call to rhotoxc
1907 : 1 : nk3xc=1
1908 : : call rhotoxc(energy,dum,kxc_cur,mpi_enreg,nfft,ngfft,&
1909 : : & nhat,usepaw,nhatgr,nhatgrdim,nkxc_cur,nk3xc,non_magnetic_xc,n3xccc,option,&
1910 : : & rhor0,rprimd,usexcnhat,vresid,vxcavg,xccc3d,xcdata,vhartr=vhres,&
1911 [ + - - + ]: 1 : & xcctau3d=xcctau3d) !vresid=work space
1912 [ + - ]: 1 : if (dtset%nspden/=4) then
1913 [ - + ]: 1 : ABI_FREE(rhor0)
1914 : : end if
1915 : :
1916 : : ! Compute Kxc(r).n^res(r)
1917 : :
1918 [ + - ]: 1 : if (dtset%nspden/=4) then
1919 : : ! Collinear magnetism or non-polarized
1920 : : call dfpt_mkvxc(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,&
1921 : 1 : & nkxc_cur,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy)
1922 : : else
1923 : : ! Non-collinear magnetism
1924 [ # # # # : 0 : ABI_MALLOC(rhor0,(nfft,dtset%nspden))
# # # # #
# # # #
# ]
1925 [ # # # # ]: 0 : rhor0(:,:)=rhor(:,:)-nresid(:,:)
1926 : : call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,&
1927 : 0 : & nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d)
1928 [ # # ]: 0 : ABI_FREE(rhor0)
1929 : : end if
1930 : :
1931 [ - + ]: 2 : ABI_FREE(kxc_cur)
1932 : : end if
1933 : :
1934 : : !if (nhatgrdim>0) then
1935 [ - + ]: 6729 : ABI_FREE(nhatgr)
1936 : : !end if
1937 : :
1938 : : !Assemble potential residual: V^res(r)=VH(n^res)(r) + Kxc(r).n^res(r)
1939 : : !--------------------------------------------------------------------
1940 [ + + ]: 14068 : do ispden=1,dtset%nspden/optnc
1941 [ + + ]: 356939049 : vresid(:,ispden)=vresid(:,ispden)+vhres(:)
1942 : : end do
1943 : :
1944 [ + - ]: 6729 : if (dtset%icoulomb==0) then
1945 [ - + ]: 6729 : ABI_FREE(nresg)
1946 : : end if
1947 [ - + ]: 6729 : ABI_FREE(vhres)
1948 [ - + ]: 6729 : ABI_FREE(dummy)
1949 : :
1950 [ + + - + : 13458 : end subroutine nres2vres
- + ]
1951 : : !!***
1952 : :
1953 : : !----------------------------------------------------------------------
1954 : :
1955 : : !!****f* ABINIT/stress_mGGA
1956 : : !!
1957 : : !! NAME
1958 : : !! stress_mGGA
1959 : : !!
1960 : : !! FUNCTION
1961 : : !! Compute non-local metaGGA contribution to stress tensor (i.e. contribution from Div.(Vtau.Grad))
1962 : : !! of (one block of) band Psi_n
1963 : : !! This contribution writes: -1/m_e * Sum_n[ f_n . Re{<Grad_alpha(Psi_n)|V_tau|Grad_beta(Psi_n)>} ]
1964 : : !!
1965 : : !! INPUTS
1966 : : !! cwavef(2,npw_k*my_nspinor*ndat)=planewave coefficients of wavefunction.
1967 : : !! effmass_free=effective mass for electrons (1. in common case)
1968 : : !! gbound_k(2*mgfft+4)=sphere boundary info
1969 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1970 : : !! istwf_k=input parameter that describes the storage of wfs
1971 : : !! kg_k(3,npw_k)=G vec coordinates wrt recip lattice transl.
1972 : : !! kpt(3)=current k point
1973 : : !! mgfft=maximum single fft dimension
1974 : : !! mpi_enreg=information about MPI parallelization
1975 : : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
1976 : : !! ndat=number of wave functions used to accumulate stress tensor (block of bands)
1977 : : !! ngfft(18)=contain all needed information about 3D FFT
1978 : : !! npw_k=number of planewaves in basis for given k point.
1979 : : !! nvtau=number of spin components of vxctaulocal
1980 : : !! n4,n5,n6=for dimensionning of vxctaulocal
1981 : : !! occ(ndat)=occupancies of bands at various k points
1982 : : !! ucvol=unit cell volume in bohr**3
1983 : : !! vxctaulocal(n4,n5,n6,nvtau,4)= local potential corresponding to the derivative of XC energy
1984 : : !! with respect to kinetic energy density, in real space, on the augmented fft grid.
1985 : : !! This array contains also the gradient of vxctaulocal (gvxctaulocal) in vxctaulocal(:,:,:,:,2:4).
1986 : : !! wtk=weights associated with current k-point
1987 : : !! [gpu_option]=GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU) [default=0]
1988 : : !!
1989 : : !! OUTPUT
1990 : : !!
1991 : : !! SIDE EFFECTS
1992 : : !! mggastr(6)=metaGGA contribution to stress tensor updated
1993 : : !!
1994 : : !! SOURCE
1995 : :
1996 : 197 : subroutine stress_mGGA(mggastr,cwavef,effmass_free,gbound_k,gprimd,istwf_k,kg_k,kpt,mgfft, &
1997 : 197 : & mpi_enreg,my_nspinor,ndat,ngfft,npw_k,nvtau,n4,n5,n6,occ,ucvol,vxctaulocal,wtk, &
1998 : : & gpu_option) ! optional
1999 : :
2000 : : !Arguments ------------------------------------
2001 : : !scalars
2002 : : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvtau,n4,n5,n6
2003 : : integer,intent(in),optional :: gpu_option
2004 : : real(dp),intent(in) :: effmass_free,ucvol,wtk
2005 : : type(MPI_type),intent(in) :: mpi_enreg
2006 : : !arrays
2007 : : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
2008 : : real(dp),intent(inout),target :: cwavef(2,npw_k*my_nspinor*ndat)
2009 : : real(dp),intent(in) :: gprimd(3,3),kpt(3),occ(ndat)
2010 : : real(dp),intent(inout) :: vxctaulocal(n4,n5,n6,nvtau,4)
2011 : : real(dp),intent(inout) :: mggastr(6)
2012 : :
2013 : : !Local variables-------------------------------
2014 : : !scalars
2015 : : integer,parameter :: opt_fourwf=2,tim_fourwf=1
2016 : : integer :: gpu_option_,ia,ib,idat,idir,ierr,ipw,ispinor,mu,nspinortot
2017 : : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
2018 : : real(dp) :: gp2pi1,gp2pi2,gp2pi3,kg_k_cart,kpt_cart,renorm_factor,weight_dum=1
2019 : : !arrays
2020 : : integer,parameter :: voigt1(6)=[1,2,3,3,3,2],voigt2(6)=[1,2,3,2,1,1]
2021 : 394 : real(dp) :: dotr(ndat),doti(ndat),my_mggastr(6)
2022 : 197 : real(dp),allocatable,target :: gcwavef(:,:,:),vtau_gcwavef(:,:,:)
2023 : 197 : real(dp),allocatable :: weight_array(:),work(:,:,:,:)
2024 : 197 : real(dp),pointer :: gcwavef_ndat(:,:,:,:),vtau_gcwavef_ndat(:,:,:,:)
2025 : 197 : real(dp),pointer :: my_cwavef(:,:)
2026 : :
2027 : : ! *********************************************************************
2028 : :
2029 : 197 : if (nvtau/=1) then
2030 : 0 : ABI_BUG("mGGA potential Vtau not compatible with non-collinear magnetism!")
2031 : : end if
2032 : :
2033 : : !Some inits
2034 [ + - ]: 197 : gpu_option_=0;if (present(gpu_option)) gpu_option_=gpu_option
2035 : 197 : renorm_factor=-one/effmass_free/ucvol
2036 : 197 : my_mggastr(:)=zero
2037 [ - + + - : 591 : ABI_MALLOC(weight_array,(ndat))
+ - + - -
+ - + ]
2038 [ + + ]: 509 : weight_array(1:ndat)=wtk*occ(1:ndat)
2039 : :
2040 : : !Parallelization over spinors
2041 : 197 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
2042 [ + - ]: 197 : if (mpi_enreg%paral_spinor==0) then
2043 : 197 : nspinor1TreatedByThisProc=.true.
2044 : 197 : nspinor2TreatedByThisProc=(nspinortot==2)
2045 : : else
2046 : 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
2047 : 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
2048 : : end if
2049 : :
2050 : : !Manage memory
2051 [ + + ]: 197 : if (nspinortot==1) then
2052 : 132 : my_cwavef => cwavef
2053 : : else ! nspinortot==2
2054 [ - + + - : 195 : ABI_MALLOC(my_cwavef,(2,npw_k*ndat))
+ - + - -
+ - + ]
2055 : : end if
2056 [ - + + - : 985 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
2057 [ - + + - : 788 : ABI_MALLOC(gcwavef,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
2058 [ - + + - : 394 : ABI_MALLOC(vtau_gcwavef,(2,npw_k*ndat,3))
- + - + ]
2059 [ + + ]: 985 : call c_f_pointer(c_loc(gcwavef(1,1,1)),gcwavef_ndat,[2,npw_k,ndat,3])
2060 [ + + ]: 985 : call c_f_pointer(c_loc(vtau_gcwavef(1,1,1)),vtau_gcwavef_ndat,[2,npw_k,ndat,3])
2061 : : #ifdef HAVE_OPENMP_OFFLOAD
2062 : : !$OMP TARGET ENTER DATA MAP(alloc:gcwavef,vtau_gcwavef,dotr,doti) IF(gpu_option_==ABI_GPU_OPENMP)
2063 : : !$OMP TARGET ENTER DATA MAP(alloc:my_cwavef) IF(gpu_option_==ABI_GPU_OPENMP .and. nspinortot==2)
2064 : : !$OMP TARGET ENTER DATA MAP(to:cwavef) IF(gpu_option_==ABI_GPU_OPENMP)
2065 : : #endif
2066 : :
2067 : : !Loop over spinors (if any)
2068 [ + + ]: 459 : do ispinor=1,my_nspinor
2069 : :
2070 : : ! Select spinor component of WF
2071 [ + + ]: 262 : if (nspinortot==2) then
2072 [ + + ]: 130 : if (ispinor==1.and.nspinor1TreatedByThisProc) then
2073 [ + - ]: 65 : if(gpu_option_==ABI_GPU_OPENMP) then
2074 : : #ifdef HAVE_OPENMP_OFFLOAD
2075 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:my_cwavef,cwavef) PRIVATE(idat)
2076 : : do idat=1,ndat
2077 : : !$OMP PARALLEL DO PRIVATE(ipw)
2078 : : do ipw=1,npw_k
2079 : : my_cwavef(1,ipw+(idat-1)*npw_k)=cwavef(1,ipw+(idat-1)*my_nspinor*npw_k)
2080 : : my_cwavef(2,ipw+(idat-1)*npw_k)=cwavef(2,ipw+(idat-1)*my_nspinor*npw_k)
2081 : : end do
2082 : : end do
2083 : : #endif
2084 : : else
2085 [ + + ]: 169 : do idat=1,ndat
2086 [ + + ]: 67249 : do ipw=1,npw_k
2087 [ + + + + ]: 402584 : my_cwavef(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k)
2088 : : end do
2089 : : end do
2090 : : end if ! gpu_option_
2091 [ + - ]: 65 : else if (ispinor==2.and.nspinor2TreatedByThisProc) then
2092 [ + - ]: 65 : if(gpu_option_==ABI_GPU_OPENMP) then
2093 : : #ifdef HAVE_OPENMP_OFFLOAD
2094 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:my_cwavef,cwavef) PRIVATE(idat)
2095 : : do idat=1,ndat
2096 : : !$OMP PARALLEL DO PRIVATE(ipw)
2097 : : do ipw=1,npw_k
2098 : : my_cwavef(1,ipw+(idat-1)*npw_k)=cwavef(1,ipw+(idat-1)*my_nspinor*npw_k+npw_k)
2099 : : my_cwavef(2,ipw+(idat-1)*npw_k)=cwavef(2,ipw+(idat-1)*my_nspinor*npw_k+npw_k)
2100 : : end do
2101 : : end do
2102 : : #endif
2103 : : else
2104 [ + + ]: 169 : do idat=1,ndat
2105 [ + + ]: 67249 : do ipw=1,npw_k
2106 [ + + + + ]: 402584 : my_cwavef(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k+npw_k)
2107 : : end do
2108 : : end do
2109 : : end if ! gpu_option_
2110 : : else
2111 : : cycle
2112 : : end if
2113 : : end if
2114 : :
2115 : : ! Loop over cartesian directions
2116 [ + + ]: 1048 : do idir=1,3
2117 : 786 : gp2pi1=gprimd(idir,1)*two_pi
2118 : 786 : gp2pi2=gprimd(idir,2)*two_pi
2119 : 786 : gp2pi3=gprimd(idir,3)*two_pi
2120 : 786 : kpt_cart=gp2pi1*kpt(1)+gp2pi2*kpt(2)+gp2pi3*kpt(3)
2121 : :
2122 : : ! Compute grad of WF (multiplication by 2pi i (G+k)_idir in reciprocal space)
2123 [ + - ]: 786 : if(gpu_option_==ABI_GPU_OPENMP) then
2124 : : #ifdef HAVE_OPENMP_OFFLOAD
2125 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:my_cwavef,gcwavef,kg_k) PRIVATE(idat)
2126 : : do idat=1,ndat
2127 : : !$OMP PARALLEL DO PRIVATE(ipw,kg_k_cart)
2128 : : do ipw=1,npw_k
2129 : : kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
2130 : : gcwavef(1,ipw+(idat-1)*npw_k,idir)= my_cwavef(2,ipw+(idat-1)*npw_k)*kg_k_cart
2131 : : gcwavef(2,ipw+(idat-1)*npw_k,idir)=-my_cwavef(1,ipw+(idat-1)*npw_k)*kg_k_cart
2132 : : end do
2133 : : end do
2134 : : #endif
2135 : : else
2136 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(idat,ipw,kg_k_cart)
2137 [ + + ]: 2034 : do idat=1,ndat
2138 [ + + ]: 909816 : do ipw=1,npw_k
2139 : 907782 : kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
2140 : 907782 : gcwavef(1,ipw+(idat-1)*npw_k,idir)= my_cwavef(2,ipw+(idat-1)*npw_k)*kg_k_cart
2141 : 909030 : gcwavef(2,ipw+(idat-1)*npw_k,idir)=-my_cwavef(1,ipw+(idat-1)*npw_k)*kg_k_cart
2142 : : end do
2143 : : end do
2144 : : end if
2145 : :
2146 : : ! Compute vxctaulocal*(grad of WF) in reciprocal space
2147 : : call fourwf(1,vxctaulocal(:,:,:,:,1),gcwavef(:,:,idir),vtau_gcwavef(:,:,idir), &
2148 : : & work,gbound_k,gbound_k,istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft, &
2149 : : & npw_k,npw_k,n4,n5,n6,opt_fourwf,tim_fourwf,weight_dum,weight_dum, &
2150 : 1048 : & weight_array_r=weight_array,gpu_option=gpu_option_)
2151 : :
2152 : : end do ! idir
2153 : :
2154 : : ! Accumulate stress tensor components (Re{<Grad_alpha(Psi_n)|V_tau.Grad_beta(Psi_n)>})
2155 [ + + ]: 2031 : do mu=1,6
2156 : 1572 : ia=voigt1(mu) ; ib=voigt2(mu)
2157 : : call dotprod_g_batch_full(dotr,doti,istwf_k,npw_k,ndat,1, &
2158 : : & gcwavef_ndat(:,:,:,ia),vtau_gcwavef_ndat(:,:,:,ib), &
2159 [ - + - + ]: 1572 : & mpi_enreg%me_g0,mpi_enreg%comm_fft,gpu_option=gpu_option_)
2160 : : #ifdef HAVE_OPENMP_OFFLOAD
2161 : : !$OMP TARGET UPDATE FROM(dotr,doti) IF(gpu_option_==ABI_GPU_OPENMP)
2162 : : #endif
2163 [ + + ]: 4330 : my_mggastr(mu)=my_mggastr(mu) + renorm_factor*sum(weight_array(1:ndat)*dotr(1:ndat))
2164 : : end do
2165 : :
2166 : : end do ! ispinor
2167 : :
2168 : : #ifdef HAVE_OPENMP_OFFLOAD
2169 : : !$OMP TARGET EXIT DATA MAP(delete:cwavef,gcwavef,vtau_gcwavef,dotr,doti) IF(gpu_option_==ABI_GPU_OPENMP)
2170 : : !$OMP TARGET EXIT DATA MAP(delete:my_cwavef) IF(gpu_option_==ABI_GPU_OPENMP .and. nspinortot==2)
2171 : : #endif
2172 : : !Release memory
2173 [ + + ]: 197 : if (nspinortot==2) then
2174 [ - + ]: 65 : ABI_FREE(my_cwavef)
2175 : : end if
2176 [ - + ]: 197 : ABI_FREE(gcwavef)
2177 [ - + ]: 197 : ABI_FREE(vtau_gcwavef)
2178 [ - + ]: 197 : ABI_FREE(weight_array)
2179 [ - + ]: 197 : ABI_FREE(work)
2180 : :
2181 : : !Take into account MPI parallelism (bands, spinors)
2182 : 197 : call xmpi_sum(my_mggastr,mpi_enreg%comm_bandspinor ,ierr)
2183 : :
2184 : : !Final accumulation of stresses
2185 [ + + ]: 1379 : mggastr(1:6) = mggastr(1:6) + my_mggastr(1:6)
2186 : :
2187 [ - + ]: 197 : end subroutine stress_mGGA
2188 : : !!***
2189 : :
2190 : : end module m_forstr
2191 : : !!***
|
