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