Line data Source code
1 : !!****m* ABINIT/m_dfpt_vtorho
2 : !! NAME
3 : !! m_dfpt_vtorho
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, AR, DRH, MB, XW, MT)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : ! nvtx related macro definition
22 : #include "nvtx_macros.h"
23 :
24 : module m_dfpt_vtorho
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_xmpi
29 : use m_errors
30 : use m_efield
31 : use m_hamiltonian
32 : use m_wfk
33 : use m_cgtools
34 : use m_dtset
35 : use m_dtfil
36 : use m_ompgpu_utils
37 :
38 :
39 : use defs_datatypes, only : pseudopotential_type
40 : use defs_abitypes, only : MPI_type
41 : use m_time, only : timab
42 : use m_occ, only : occeig
43 : use m_hdr, only : hdr_skip, hdr_io
44 : use m_pawang, only : pawang_type
45 : use m_pawtab, only : pawtab_type
46 : use m_paw_ij, only : paw_ij_type
47 : use m_pawfgrtab,only : pawfgrtab_type
48 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, &
49 : & pawrhoij_init_unpacked, pawrhoij_free_unpacked, &
50 : & pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim
51 : use m_pawcprj, only : pawcprj_type
52 : use m_pawfgr, only : pawfgr_type
53 : use m_paw_mkrho,only : pawmkrho
54 : use m_fft, only : fftpac
55 : use m_spacepar, only : symrhg
56 : use m_getgh1c, only : rf_transgrid_and_pack, getgh1c_setup
57 : use m_dfpt_vtowfk, only : dfpt_vtowfk
58 : use m_dfpt_fef, only : dfptff_gradberry, dfptff_gbefd
59 : use m_mpinfo, only : proc_distrb_cycle,proc_distrb_nband
60 : use m_fourier_interpol, only : transgrid
61 : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, gemm_nonlop_use_gemm
62 :
63 : #if defined(HAVE_GPU_MARKERS)
64 : use m_nvtx_data
65 : #endif
66 :
67 : implicit none
68 :
69 : private
70 : !!***
71 :
72 : public :: dfpt_vtorho
73 : !!***
74 :
75 : contains
76 : !!***
77 :
78 : !!****f* ABINIT/dfpt_vtorho
79 : !! NAME
80 : !! dfpt_vtorho
81 : !!
82 : !! FUNCTION
83 : !! This routine compute the new 1-density from a fixed 1-potential (vtrial1)
84 : !! but might also simply compute eigenvectors and eigenvalues.
85 : !! The main part of it is a wf update over all k points
86 : !!
87 : !! INPUTS
88 : !! cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions
89 : !! cgq(2,mpw1*nspinor*mband_mem*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
90 : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
91 : !! cplex: if 1, real space 1-order functions on FFT grid are REAL; if 2, COMPLEX
92 : !! TODO distribute cprj over bands mband_mem
93 : !! cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
94 : !! projected with non-local projectors: cprj=<p_i|Cnk>
95 : !! cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q
96 : !! projected with non-local projectors: cprjq=<p_i|Cnk+q>
97 : !! dbl_nnsclo=if 1, will double the value of dtset%nnsclo
98 : !! doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
99 : !! docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
100 : !! dtefield = variables related to response Berry-phase calculation
101 : !! dtfil <type(datafiles_type)>=variables related to files
102 : !! dtset <type(dataset_type)>=all input variables for this dataset
103 : !! eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
104 : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
105 : !! fermie1=derivative of fermi energy wrt (strain) perturbation
106 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
107 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
108 : !! idir=direction of the perturbation
109 : !! indsy1(4,nsym1,natom)=indirect indexing array for atom labels
110 : !! ipert=type of the perturbation
111 : !! irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
112 : !! istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
113 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
114 : !! kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
115 : !! kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
116 : !! mband=maximum number of bands
117 : !! mband_mem=maximum number of bands on this cpu
118 : !! mkmem =number of k points treated by this node (GS data).
119 : !! mkqmem =number of k+q points treated by this node (GS data)
120 : !! mk1mem =number of k points treated by this node (RF data)
121 : !! mpw=maximum dimensioned size of npw or wfs at k
122 : !! mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
123 : !! my_natom=number of atoms treated by current processor
124 : !! natom=number of atoms in cell.
125 : !! nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
126 : !! ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
127 : !! nfftf= -PAW ONLY- number of FFT grid points for the fine grid
128 : !! (nfftf=nfft for norm-conserving potential runs - see comment in respfn.F90)
129 : !! nkpt_rbz=number of k points in the IBZ for this perturbation
130 : !! mpi_enreg=information about MPI parallelization
131 : !! npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
132 : !! npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
133 : !! nspden=number of spin-density components
134 : !! nsppol=1 for unpolarized, 2 for spin-polarized
135 : !! nsym1=number of symmetry elements in space group consistent with
136 : !! perturbation
137 : !! ntypat=number of types of atoms in unit cell.
138 : !! occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
139 : !! at each k+q point of the reduced Brillouin zone.
140 : !! occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k
141 : !! (usually 2)
142 : !! omega= frequency of the perturbation
143 : !! optres=0: the new value of the density is computed in place of the input value
144 : !! 1: only the density residual is computed ; the input density is kept
145 : !! paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
146 : !! paw_ij1(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels
147 : !! pawang <type(pawang_type)>=paw angular mesh and related data
148 : !! pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
149 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
150 : !! pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
151 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
152 : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
153 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
154 : !! phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
155 : !! ph1d(2,3*(2*dtset%mgfft+1)*natom)=one-dimensional structure factor information
156 : !! prtvol=control print volume and debugging output
157 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
158 : !! pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices
159 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
160 : !! inverse of the overlap matrix
161 : !! rmet(3,3)=real space metric (bohr**2)
162 : !! rprimd(3,3)=dimensional real space primitive translations
163 : !! symaf1(nsym1)=(anti)ferromagnetic part of symmetry operations
164 : !! symrc1(3,3,nsym1)=symmetry operations in reciprocal space
165 : !! symrl1(3,3,nsym1)=symmetry operations in real space
166 : !! tnons1(3,nsym1)=non-symmorphic translations
167 : !! ucvol=unit cell volume in bohr**3.
168 : !! usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
169 : !! useylmgr1= 1 if ylmgr1 array is allocated
170 : !! usevxctau=1 if if XC functional depends on kinetic energy density
171 : !! ddk<wfk_t)=struct info DDK file
172 : !! vectornd(with_vectornd*nfftf,nspden,3)=nuclear dipole moment vector potential
173 : !! vtrial(nfftf,nspden)=GS Vtrial(r).
174 : !! vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r).
175 : !! vxctau(nfftf,nspden,4*usevxctau)=derivative of e_xc wrt kin energy density, for mGGA
176 : !! with_vectornd = 1 if vectornd allocated
177 : !! wtk_rbz(nkpt_rbz)=weight assigned to each k point.
178 : !! xred(3,natom)=reduced dimensionless atomic coordinates
179 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
180 : !! ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+g point
181 : !! ylmgr1(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics for each G and k+g point
182 : !!
183 : !!
184 : !! OUTPUT
185 : !! cg1(2,mpw*nspinor*mband_mem*mk1mem*nsppol)=updated wavefunctions, orthogonalized to the occupied states
186 : !! cg1_active(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF
187 : !! wavefunctions at k,q. They are orthogonalized to the active.
188 : !! eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues
189 : !! (hartree)
190 : !! edocc=correction to 2nd-order total energy coming from changes of occupation
191 : !! eeig0=0th-order eigenenergies part of 2nd-order total energy
192 : !! ek0=0th-order kinetic energy part of 2nd-order total energy.
193 : !! ek1=1st-order kinetic energy part of 2nd-order total energy
194 : !! (not for phonons)
195 : !! eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
196 : !! end0=0th-order nuclear dipole energy part of 2nd-order total energy.
197 : !! end1=1st-order nuclear dipole energy part of 2nd-order total energy
198 : !! enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
199 : !! enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
200 : !! evxctau0=0th-order energy from vxctau
201 : !! evxctau1=1st-order energy from vxctau
202 : !! gh1c_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
203 : !! gh0c1_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}|\Psi^{(1)}>
204 : !! The wavefunction is orthogonal to the active space (for metals). It is not
205 : !! coherent with cg1.
206 : !! resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points.
207 : !! residm=maximum value from resid array (except for nbdbuf highest bands)
208 : !! rhog1(2,nfftf)=RF electron density in reciprocal space
209 : !! ==== if optres==1
210 : !! nres2=square of the norm of the residual
211 : !! nvresid1(cplex*nfftf,nspden)=1st-order density residual
212 : !! ==== if psps%usepaw==1
213 : !! cprj1(natom,nspinor*mband_mem*mk1mem*nsppol*usecprj)=
214 : !! 1st-order wave functions at k,q projected with non-local projectors:
215 : !! cprj1=<p_i|C1nk,q> where p_i is a non-local projector
216 : !! nhat1(cplex*nfftf,nspden*psps%usepaw)=1st-order compensation charge density
217 : !!
218 : !! SIDE EFFECTS
219 : !! pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
220 : !! rhor1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3.
221 : !!
222 : !! SOURCE
223 :
224 18776 : subroutine dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dbl_nnsclo,&
225 18776 : & dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,qphon,&
226 18776 : & edocc,eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,end0,end1,enl0,enl1,&
227 18776 : & evxctau0,evxctau1,fermie1,gh0c1_set,gh1c_set,gmet,gprimd,idir,indsy1,&
228 18776 : & ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,mband_mem,&
229 : & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,&
230 18776 : & natom,nband_rbz,ncpgr,nfftf,nhat1,nkpt_rbz,npwarr,npwar1,nres2,nspden,&
231 18776 : & nsppol,nsym1,ntypat,nvresid1,occkq,occ_rbz,optres,&
232 18776 : & paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,pawrhoij1,pawtab,&
233 18776 : & phnons1,ph1d,prtvol,psps,pwindall,qmat,resid,residm,rhog1,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,ucvol,&
234 18776 : & usecprj,useylmgr1,usevxctau,ddk_f,vectornd,vtrial,vtrial1,vxctau,with_vectornd,wtk_rbz,xred,ylm,ylm1,ylmgr1,cg1_out,eta,omega)
235 :
236 : !Arguments -------------------------------
237 : !scalars
238 : integer,intent(in) :: cplex,dbl_nnsclo,dim_eig2rf,idir,ipert,mband,mk1mem,mkmem
239 : integer,intent(in) :: mband_mem
240 : integer,intent(in) :: mkqmem,mpw,mpw1,my_natom,natom,ncpgr,nfftf,nkpt_rbz,nspden
241 : integer,intent(in) :: nsppol,nsym1,ntypat,optres,prtvol,usecprj,useylmgr1,usevxctau,with_vectornd
242 : integer,optional,intent(in) :: cg1_out
243 : real(dp),intent(in) :: fermie1,ucvol
244 : real(dp),optional,intent(in) :: eta,omega
245 : real(dp),intent(out) :: edocc,eeig0,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,nres2,residm
246 : type(MPI_type),intent(in) :: mpi_enreg
247 : type(datafiles_type),intent(in) :: dtfil
248 : type(dataset_type),intent(in) :: dtset
249 : type(efield_type),intent(in) :: dtefield
250 : type(pawang_type),intent(in) :: pawang,pawang1
251 : type(pawfgr_type),intent(in) :: pawfgr
252 : type(pseudopotential_type),intent(in) :: psps
253 :
254 : !arrays
255 : integer,intent(in) :: indsy1(4,nsym1,natom)
256 : ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise
257 : integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))
258 : integer,intent(in) :: istwfk_rbz(nkpt_rbz)
259 : integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
260 : integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),npwar1(nkpt_rbz,2)
261 : integer,intent(in) :: npwarr(nkpt_rbz,2),symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1)
262 : real(dp),intent(in) :: qphon(3)
263 : real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem*mkmem*nsppol)
264 : real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol)
265 : real(dp),intent(inout):: cg1_active(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)
266 : real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)
267 : real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)
268 : real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol)
269 : real(dp),intent(in) :: doccde_rbz(mband*nkpt_rbz*nsppol)
270 : real(dp),intent(in) :: docckqde(mband*nkpt_rbz*nsppol)
271 : real(dp),intent(in) :: eigen0(mband*nkpt_rbz*nsppol)
272 : real(dp),intent(out) :: eigen1(2*mband*mband*nkpt_rbz*nsppol)
273 : real(dp),intent(in) :: eigenq(mband*nkpt_rbz*nsppol),gmet(3,3),gprimd(3,3)
274 : real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),occ_rbz(mband*nkpt_rbz*nsppol)
275 : real(dp),intent(in) :: occkq(mband*nkpt_rbz*nsppol),ph1d(2,3*(2*dtset%mgfft+1)*natom)
276 : real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))
277 : real(dp), intent(out) :: nhat1(cplex*nfftf,dtset%nspden*psps%usepaw)
278 : real(dp),intent(out) :: resid(mband*nkpt_rbz*nsppol),rhog1(2,nfftf)
279 : real(dp),intent(inout) :: nvresid1(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden)
280 : real(dp),intent(inout) :: vectornd(with_vectornd*nfftf,dtset%nspden,3)
281 : real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
282 : real(dp),intent(in) :: tnons1(3,nsym1)
283 : real(dp),intent(in),target :: vtrial(nfftf,nspden)
284 : real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden)
285 : real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*usevxctau)
286 : real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom)
287 : real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
288 : real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
289 : real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
290 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
291 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt_rbz,2,3)
292 : type(pawcprj_type),intent(in) :: cprj (natom,dtset%nspinor*mband_mem*mkmem *nsppol*usecprj)
293 : type(pawcprj_type),intent(in) :: cprjq(natom,dtset%nspinor*mband_mem*mkqmem*nsppol*usecprj)
294 : type(pawcprj_type),intent(inout) :: cprj1(natom,dtset%nspinor*mband_mem*mk1mem*nsppol*usecprj)
295 : type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw),paw_ij1(my_natom*psps%usepaw)
296 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
297 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw)
298 : type(pawrhoij_type),target,intent(inout) :: pawrhoij1(my_natom*psps%usepaw)
299 : type(pawtab_type), intent(in) :: pawtab(ntypat*psps%usepaw)
300 : type(wfk_t),intent(inout) :: ddk_f(4)
301 :
302 : !Local variables-------------------------------
303 : !scalars
304 : integer,parameter :: level=13
305 : integer :: bd2tot_index,bdtot_index,buffer_size,counter,cplex_rhoij
306 : integer :: iband,nlines_done,ibdkpt,ibg,ibg1,ibgq,icg,icg1,icgq,ierr
307 : integer :: ii,ikg,ikg1,ikpt,ilm,index1,ispden,iscf_mod,isppol,istwf_k
308 : integer :: mbd2kpsp,mbdkpsp,mcgq,mcgq_disk,mcprjq
309 : integer :: mcprjq_disk,me,n1,n2,n3,n4,n5,n6,nband_k,nband_kq,nkpg,nkpg1
310 : integer :: nband_eff
311 : integer :: nnsclo_now,npw1_k,npw_k,nspden_rhoij,qphase_rhoij,spaceworld,test_dot
312 : integer :: nband_me
313 : logical :: has_vectornd,paral_atom,qne0,with_vxctau
314 : real(dp) :: arg,eta_,omega_,wtk_k
315 18776 : type(gs_hamiltonian_type) :: gs_hamkq
316 18776 : type(rf_hamiltonian_type) :: rf_hamkq,rf_hamk_dir2
317 : !arrays
318 18776 : integer,allocatable :: kg1_k(:,:),kg_k(:,:)
319 18776 : integer, pointer :: my_atmtab(:)
320 : real(dp) :: kpoint(3),kpq(3)
321 : real(dp) :: tsec(2)
322 18776 : real(dp),allocatable :: buffer1(:)
323 18776 : real(dp),allocatable :: ddkinpw(:),dkinpw(:),dkinpw2(:)
324 18776 : real(dp),allocatable :: doccde_k(:),doccde_kq(:)
325 18776 : real(dp),allocatable :: edocc_k(:),eeig0_k(:),eig0_k(:),eig0_kq(:),eig1_k(:)
326 18776 : real(dp),allocatable :: ek0_k(:),ek1_k(:),eloc0_k(:),end0_k(:),end1_k(:),enl0_k(:),enl1_k(:)
327 18776 : real(dp),allocatable :: evxctau0_k(:),evxctau1_k(:)
328 18776 : real(dp),allocatable :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:),ffnlk(:,:,:,:)
329 18776 : real(dp),allocatable :: grad_berry(:,:,:),kinpw1(:),kpg1_k(:,:)
330 18776 : real(dp),allocatable :: kpg_k(:,:),occ_k(:),occ_kq(:)
331 18776 : real(dp),allocatable :: ph3d(:,:,:),ph3d1(:,:,:),resid_k(:)
332 18776 : real(dp),allocatable :: rho1wfg(:,:),rho1wfr(:,:),rhoaug1(:,:,:,:),rocceig(:,:)
333 18776 : real(dp),allocatable :: vectornd_pac(:,:,:,:,:),vectornd_pac_idir(:,:,:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:)
334 18776 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
335 18776 : real(dp),allocatable :: ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:)
336 18776 : type(pawrhoij_type),pointer :: pawrhoij1_unsym(:)
337 : ! *********************************************************************
338 :
339 : DBG_ENTER('COLL')
340 :
341 : ABI_NVTX_START_RANGE(NVTX_DFPT_VTORHO)
342 :
343 : !Keep track of total time spent in this routine
344 18776 : call timab(121,1,tsec)
345 18776 : call timab(124,1,tsec)
346 :
347 : !Retrieve parallelism data
348 18776 : spaceworld=mpi_enreg%comm_cell
349 18776 : me=mpi_enreg%me_kpt
350 18776 : paral_atom=(my_natom/=natom)
351 18776 : my_atmtab=>mpi_enreg%my_atmtab
352 :
353 : if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
354 18776 : & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
355 388 : ABI_MALLOC(grad_berry,(2,mpw1,dtefield%mband_occ))
356 : else
357 18679 : ABI_MALLOC(grad_berry,(0,0,0))
358 : end if
359 :
360 : !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
361 18776 : if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
362 0 : ABI_BUG('wrong values for nfft, nfftf!')
363 : end if
364 :
365 : !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f
366 18776 : iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3
367 :
368 18776 : edocc=zero ; eeig0=zero ; ek0=zero ; ek1=zero
369 18776 : eloc0=zero ; end0=zero ; end1=zero ; enl0=zero ; enl1=zero
370 18776 : evxctau0=zero; evxctau1=zero
371 18776 : bdtot_index=0
372 18776 : bd2tot_index=0
373 18776 : ibg=0;icg=0
374 18776 : ibgq=0;icgq=0
375 18776 : ibg1=0;icg1=0
376 18776 : mbdkpsp=mband*nkpt_rbz*nsppol
377 18776 : mbd2kpsp=2*mband**2*nkpt_rbz*nsppol
378 :
379 18776 : n1=dtset%ngfft(1); n2=dtset%ngfft(2); n3=dtset%ngfft(3)
380 18776 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6)
381 18776 : qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=tol14)
382 :
383 18776 : omega_=zero; if(present(omega))omega_=omega
384 18776 : eta_=zero; if(present(eta))eta_=eta
385 :
386 : !Initialize PW 1st-order density if needed
387 : !Also store old rho1 in case of density mixing
388 18776 : if (iscf_mod>0) then
389 15118 : if (optres==1) nvresid1=rhor1
390 15118 : if (psps%usepaw==0) then
391 160367912 : rhor1(:,:)=zero
392 : else
393 10588 : ABI_MALLOC(rho1wfr,(cplex*dtset%nfft,dtset%nspden))
394 7941 : ABI_MALLOC(rho1wfg,(2,dtset%nfft))
395 25152725 : rho1wfr(:,:)=zero
396 : end if
397 : end if
398 :
399 : !Set max number of non-self-consistent loops nnsclo_now for use in dfpt_vtowfk
400 18776 : if(iscf_mod<=0 .and. iscf_mod/=-3)then
401 124 : nnsclo_now=dtset%nstep
402 : else
403 18652 : if(dtset%nnsclo>0)then
404 0 : nnsclo_now=dtset%nnsclo
405 : else
406 18652 : nnsclo_now=1
407 : end if
408 18652 : if(dbl_nnsclo==1) nnsclo_now=nnsclo_now*2
409 : end if
410 :
411 : !Prepare GS k+q wf
412 18776 : mcgq=mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol;mcgq_disk=0
413 :
414 : !Prepare RF PAW files
415 18776 : if (psps%usepaw==1) then
416 3496 : mcprjq=dtset%nspinor*mband_mem*mkqmem*nsppol*usecprj;mcprjq_disk=0
417 : else
418 15280 : mcprjq=0;mcprjq_disk=0
419 : end if
420 :
421 : !Initialisation of the wfdot file in case of electric field (or 2nd order Sternheimer equation)
422 18776 : test_dot=0
423 : if (ipert==natom+2.and.sum((qphon(1:3))**2 )<=tol7.and.&
424 : & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
425 : & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17).or.&
426 : & (ipert==natom+10.or.ipert==natom+11)) then
427 18776 : test_dot=1
428 : end if
429 :
430 : !==== Initialize most of the Hamiltonian (and derivative) ====
431 : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
432 : !2) Perform the setup needed for the non-local factors:
433 : !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
434 : !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
435 :
436 : call gs_hamkq%init(psps,pawtab,dtset%nspinor,nsppol,nspden,natom,&
437 : & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
438 : & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
439 18776 : & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option)
440 :
441 : call rf_hamkq%init(cplex,gs_hamkq,ipert,has_e1kbsc=.true.,paw_ij1=paw_ij1,&
442 18776 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
443 18776 : if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then
444 : call rf_hamk_dir2%init(cplex,gs_hamkq,ipert,has_e1kbsc=.true.,paw_ij1=paw_ij1,&
445 130 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
446 : end if
447 :
448 : !PAW:allocate memory for non-symetrized 1st-order occupancies matrix (pawrhoij1)
449 18776 : pawrhoij1_unsym => pawrhoij1
450 18776 : if (psps%usepaw==1.and.iscf_mod>0) then
451 2647 : if (paral_atom) then
452 9170 : ABI_MALLOC(pawrhoij1_unsym,(natom))
453 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
454 1254 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
455 : call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
456 : & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,&
457 1254 : & use_rhoijp=0,use_rhoij_=1)
458 : else
459 1393 : pawrhoij1_unsym => pawrhoij1
460 1393 : call pawrhoij_init_unpacked(pawrhoij1_unsym)
461 : end if
462 : end if
463 :
464 112656 : ABI_MALLOC(rhoaug1,(cplex*n4,n5,n6,gs_hamkq%nvloc))
465 112656 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
466 93880 : ABI_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc))
467 :
468 18776 : with_vxctau = ( usevxctau > 0 )
469 18776 : if(with_vxctau) then
470 294 : ABI_MALLOC(vxctaulocal,(n4,n5,n6,gs_hamkq%nvloc,4))
471 : end if
472 :
473 18776 : has_vectornd = (with_vectornd .EQ. 1)
474 18776 : if(has_vectornd) then
475 2562 : ABI_MALLOC(vectornd_pac,(n4,n5,n6,gs_hamkq%nvloc,3))
476 1830 : ABI_MALLOC(vectornd_pac_idir,(n4,n5,n6,gs_hamkq%nvloc))
477 : end if
478 :
479 18776 : nlines_done = 0
480 5410696 : resid = zero
481 :
482 : !LOOP OVER SPINS
483 39139 : do isppol=1,nsppol
484 :
485 : ! Rewind kpgsph data file if needed:
486 20363 : ikg=0;ikg1=0
487 :
488 : ! Set up local potential vlocal1 with proper dimensioning, from vtrial1
489 : ! Same thing for vlocal from vtrial Also take into account the spin.
490 :
491 : call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,&
492 20363 : gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1)
493 :
494 : ! Continue to initialize the Hamiltonian
495 20363 : call gs_hamkq%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.)
496 20363 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.)
497 20363 : if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then
498 170 : call rf_hamk_dir2%load_spin(isppol,with_nonlocal=.true.)
499 170 : if (ipert==natom+11) then ! load vlocal1
500 78 : call rf_hamk_dir2%load_spin(isppol,vlocal1=vlocal1)
501 : end if
502 : end if
503 :
504 20363 : if (ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11)) then !SPr deb, in case of magnetic field perturbation, no non-local
505 169 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1)
506 : end if
507 :
508 : ! Nullify contribution to 1st-order density from this k-point
509 246495323 : rhoaug1(:,:,:,:)=zero
510 :
511 : ! if vectornd is present, set it up for addition to gs_hamkq and rf_hamkq.
512 : ! Note that it must be done for the three Cartesian directions. Also, the following
513 : ! code assumes explicitly and implicitly that nvloc = 1. This should eventually be generalized.
514 20363 : if(has_vectornd) then
515 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
516 390 : & dtset%nspden, gs_hamkq%nvloc, 3, pawfgr, mpi_enreg, vectornd, vectornd_pac)
517 390 : call gs_hamkq%load_spin(isppol, vectornd=vectornd_pac)
518 12273369 : vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,idir)
519 390 : call rf_hamkq%load_spin(isppol, vectornd=vectornd_pac_idir)
520 : end if
521 :
522 : !! add vxctau for mGGA to GS hamiltonian and RF hamiltonian
523 20363 : if (with_vxctau) then
524 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
525 42 : dtset%nspden, gs_hamkq%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
526 42 : call gs_hamkq%load_spin(isppol, vxctaulocal=vxctaulocal)
527 42 : call rf_hamkq%load_spin(isppol, vxctaulocal=vxctaulocal)
528 : end if
529 :
530 20363 : call timab(125,1,tsec)
531 :
532 : !======================================================================
533 : !============== BIG FAT K POINT LOOP ================================
534 : !======================================================================
535 :
536 1045820 : do ikpt=1,nkpt_rbz
537 1025457 : counter=100*ikpt+isppol
538 :
539 1025457 : nband_k = nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
540 : ! enables variable nband less than the block size in serial case
541 1025457 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
542 1025457 : istwf_k = istwfk_rbz(ikpt)
543 1025457 : npw_k = npwarr(ikpt,1)
544 1025457 : npw1_k = npwar1(ikpt,1)
545 1025457 : wtk_k = wtk_rbz(ikpt)
546 :
547 1025457 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
548 11134860 : eigen1(1+bd2tot_index : 2*nband_k**2+bd2tot_index) = zero
549 1134428 : resid(1+bdtot_index : nband_k+bdtot_index) = zero
550 : bdtot_index=bdtot_index+nband_k
551 : bd2tot_index=bd2tot_index+2*nband_k**2
552 :
553 : cycle ! Skip the rest of the k-point loop
554 : end if
555 :
556 3379348 : kpoint(:)=kpt_rbz(:,ikpt)
557 844837 : kpq(:)=kpoint(:)
558 844837 : if (ipert<natom+3.or.ipert==natom+5.or.ipert==natom+6.or.&
559 3104800 : &(ipert>natom+11.and.ipert<=2*natom+11)) kpq(:)=kpq(:)+qphon(1:3)
560 2534511 : ABI_MALLOC(kg_k,(3,npw_k))
561 2534511 : ABI_MALLOC(kg1_k,(3,npw1_k))
562 3379348 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
563 3379348 : ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
564 4224185 : ABI_MALLOC(ylmgr1_k,(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
565 2534511 : ABI_MALLOC(doccde_k,(nband_k))
566 1689674 : ABI_MALLOC(doccde_kq,(nband_k))
567 1689674 : ABI_MALLOC(eig0_k,(nband_k))
568 1689674 : ABI_MALLOC(eig0_kq,(nband_k))
569 2534511 : ABI_MALLOC(eig1_k,(2*nband_k**2))
570 1689674 : ABI_MALLOC(edocc_k,(nband_k))
571 1689674 : ABI_MALLOC(eeig0_k,(nband_k))
572 1689674 : ABI_MALLOC(ek0_k,(nband_k))
573 1689674 : ABI_MALLOC(ek1_k,(nband_k))
574 1689674 : ABI_MALLOC(eloc0_k,(nband_k))
575 1689674 : ABI_MALLOC(end0_k,(nband_k))
576 1689674 : ABI_MALLOC(end1_k,(nband_k))
577 1689674 : ABI_MALLOC(enl0_k,(nband_k))
578 1689674 : ABI_MALLOC(enl1_k,(nband_k))
579 1689674 : ABI_MALLOC(evxctau0_k,(nband_k))
580 1689674 : ABI_MALLOC(evxctau1_k,(nband_k))
581 1689674 : ABI_MALLOC(occ_k,(nband_k))
582 1689674 : ABI_MALLOC(occ_kq,(nband_k))
583 1689674 : ABI_MALLOC(resid_k,(nband_k))
584 3379348 : ABI_MALLOC(rocceig,(nband_k,nband_k))
585 :
586 69535129 : eig1_k(:)=zero
587 5280849 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
588 5280849 : eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
589 5280849 : edocc_k(:)=zero
590 14152873 : eeig0_k(:)=zero ; ek0_k(:)=zero ; ek1_k(:)=zero
591 14152873 : eloc0_k(:)=zero ; end0_k(:)=zero ; end1_k(:)=zero
592 9716861 : enl0_k(:)=zero ; enl1_k(:)=zero
593 9716861 : evxctau0_k(:)=zero; evxctau1_k(:)=zero
594 5280849 : occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
595 5280849 : occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index)
596 5280849 : doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index)
597 5280849 : doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index)
598 5280849 : resid_k(:)=zero
599 :
600 : ! For each pair of active bands (m,n), generates the ratios
601 : ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
602 : ! and decide to which band to attribute it.
603 844837 : call occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig)
604 :
605 : ! These arrays are not needed anymore.
606 844837 : ABI_FREE(doccde_k)
607 844837 : ABI_FREE(doccde_kq)
608 844837 : ABI_FREE(occ_kq)
609 :
610 470509413 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
611 844837 : if (psps%useylm==1) then
612 2612756 : do ilm=1,psps%mpsang*psps%mpsang
613 316306573 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
614 : end do
615 : end if
616 :
617 : ! Get (k+q+G) wave vectors and associated spherical harmonics
618 470542733 : kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
619 844837 : if (psps%useylm==1) then
620 2612756 : do ilm=1,psps%mpsang*psps%mpsang
621 316306573 : ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
622 : end do
623 268010 : if (useylmgr1==1) then
624 782386 : do ilm=1,psps%mpsang*psps%mpsang
625 3661018 : do ii=1,3+6*((ipert-natom)/10)
626 394594596 : ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
627 : end do
628 : end do
629 : end if
630 : end if
631 :
632 : ! Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
633 : call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In
634 : kpoint,kpq,idir,ipert,natom,rmet,gprimd,gmet,istwf_k,& ! In
635 : npw_k,npw1_k,useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,& ! In
636 : dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out
637 : ddkinpw=ddkinpw,dkinpw2=dkinpw2,rf_hamk_dir2=rf_hamk_dir2,& ! Optional
638 844837 : ffnl1_test=ffnl1_test) ! Optional
639 :
640 844837 : if (gs_hamkq%gpu_option == ABI_GPU_OPENMP) then
641 0 : call ompgpu_load_hamilt_buffers(kg_k,kg1_k,ffnlk,ph3d)
642 : end if
643 :
644 : ! Compute the gradient of the Berry-phase term
645 : if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
646 : & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
647 776 : if (ipert<=natom) then
648 : ! phonon perturbation
649 : call dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,mband,mband_mem,mpw,mpw1,mkmem,mk1mem,&
650 : & mpi_enreg,nkpt_rbz,&
651 488 : & npwarr,npwar1,dtset%nspinor,nsppol,qmat,pwindall)
652 : else
653 : ! electric field perturbation
654 : call dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir,ikpt,isppol,mband,mband_mem,mpw,mpw1,mkmem,mk1mem,&
655 : & mpi_enreg,nkpt_rbz,&
656 : & npwarr,npwar1,dtset%nspinor,&
657 288 : & nsppol,qmat,pwindall,rprimd)
658 : end if
659 : end if
660 :
661 : ! Setup gemm_nonlop
662 844837 : if (gemm_nonlop_use_gemm) then
663 : call set_gemm_nonlop_ikpt(ikpt,gs_hamkq%npw_fft_k,gs_hamkq%istwf_k,gs_hamkq%indlmn,&
664 0 : & gs_hamkq%ntypat,gs_hamkq%nattyp,gs_hamkq%gpu_option)
665 : end if ! gemm_nonlop_use_gemm
666 :
667 : ! Free some memory before calling dfpt_vtowfk
668 844837 : ABI_FREE(ylm_k)
669 844837 : ABI_FREE(ylm1_k)
670 844837 : ABI_FREE(ylmgr1_k)
671 :
672 : ! Compute the eigenvalues, wavefunction, residuals,
673 : ! contributions to kinetic energy, nonlocal energy, forces,
674 : ! and update of 1st-order density to this k-point and this spin polarization.
675 844837 : nband_kq = nband_k !Note that the calculation only works for same number of bands on all K points.
676 : ! Note that dfpt_vtowfk is called with kpoint, while kpt is used inside vtowfk3
677 : call dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dim_eig2rf,dtfil,&
678 : & dtset,edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,ek0_k,ek1_k,eloc0_k,end0_k,end1_k,enl0_k,enl1_k,&
679 : & evxctau0_k,evxctau1_k,fermie1,&
680 : & ffnl1,ffnl1_test,gh0c1_set,gh1c_set,grad_berry,gs_hamkq,ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,isppol,&
681 : & mband,mband_mem,mcgq,mcprjq,mkmem,mk1mem,mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,nnsclo_now,&
682 : & npw_k,npw1_k,dtset%nspinor,nsppol,n4,n5,n6,occ_k,pawrhoij1_unsym,prtvol,psps,resid_k,&
683 844837 : & rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,ddk_f,wtk_k,nlines_done,cg1_out,eta=eta_,omega=omega_)
684 :
685 844837 : if ( gs_hamkq%gpu_option == ABI_GPU_OPENMP) then
686 0 : call ompgpu_free_hamilt_buffers()
687 : end if
688 :
689 : ! Free temporary storage
690 844837 : ABI_FREE(kinpw1)
691 844837 : ABI_FREE(kg_k)
692 844837 : ABI_FREE(kg1_k)
693 844837 : ABI_FREE(kpg_k)
694 844837 : ABI_FREE(kpg1_k)
695 844837 : ABI_FREE(dkinpw)
696 844837 : if (ipert==natom+10) then
697 11160 : ABI_FREE(ddkinpw)
698 11160 : if (idir>3) then
699 7056 : ABI_FREE(dkinpw2)
700 : end if
701 : end if
702 844837 : ABI_FREE(ffnlk)
703 844837 : ABI_FREE(ffnl1)
704 844837 : ABI_SFREE(ffnl1_test)
705 844837 : ABI_FREE(eig0_k)
706 844837 : ABI_FREE(eig0_kq)
707 844837 : ABI_FREE(rocceig)
708 844837 : ABI_FREE(ph3d)
709 844837 : ABI_SFREE(ph3d1)
710 :
711 : ! Save eigenvalues (hartree), residuals (hartree**2)
712 69535129 : eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:)
713 5280849 : resid (1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
714 :
715 : ! Accumulate sum over k points for nonlocal and kinetic energies,
716 : ! also accumulate gradients of Enonlocal:
717 844837 : if (iscf_mod>0 .or. iscf_mod==-3 .or. iscf_mod==-2)then
718 5280849 : do iband=1,nband_k
719 4436012 : edocc=edocc+wtk_k*occ_k(iband)*edocc_k(iband)
720 4436012 : eeig0=eeig0+wtk_k*occ_k(iband)*eeig0_k(iband)
721 4436012 : ek0=ek0+wtk_k*occ_k(iband)*ek0_k(iband)
722 4436012 : ek1=ek1+wtk_k*occ_k(iband)*ek1_k(iband)
723 4436012 : eloc0=eloc0+wtk_k*occ_k(iband)*eloc0_k(iband)
724 4436012 : end0=end0+wtk_k*occ_k(iband)*end0_k(iband)
725 4436012 : end1=end1+wtk_k*occ_k(iband)*end1_k(iband)
726 4436012 : enl0=enl0+wtk_k*occ_k(iband)*enl0_k(iband)
727 4436012 : enl1=enl1+wtk_k*occ_k(iband)*enl1_k(iband)
728 4436012 : evxctau0=evxctau0+wtk_k*occ_k(iband)*evxctau0_k(iband)
729 5280849 : evxctau1=evxctau1+wtk_k*occ_k(iband)*evxctau1_k(iband)
730 : end do
731 : end if
732 :
733 844837 : ABI_FREE(eig1_k)
734 844837 : ABI_FREE(occ_k)
735 844837 : ABI_FREE(resid_k)
736 844837 : ABI_FREE(edocc_k)
737 844837 : ABI_FREE(eeig0_k)
738 844837 : ABI_FREE(ek0_k)
739 844837 : ABI_FREE(ek1_k)
740 844837 : ABI_FREE(eloc0_k)
741 844837 : ABI_FREE(end0_k)
742 844837 : ABI_FREE(end1_k)
743 844837 : ABI_FREE(enl0_k)
744 844837 : ABI_FREE(enl1_k)
745 844837 : ABI_FREE(evxctau0_k)
746 844837 : ABI_FREE(evxctau1_k)
747 :
748 : ! Keep track of total number of bands (all k points so far, even for k points not treated by me)
749 844837 : bdtot_index=bdtot_index+nband_k
750 844837 : bd2tot_index=bd2tot_index+2*nband_k**2
751 :
752 : ! Shift array memory
753 844837 : if (mkmem/=0) then
754 844837 : ibg=ibg+dtset%nspinor*nband_me
755 844837 : icg=icg+npw_k*dtset%nspinor*nband_me
756 844837 : ikg=ikg+npw_k
757 : end if
758 844837 : if (mkqmem/=0) then
759 844837 : ibgq=ibgq+dtset%nspinor*nband_me
760 844837 : icgq=icgq+npw1_k*dtset%nspinor*nband_me
761 : end if
762 865200 : if (mk1mem/=0) then
763 844837 : ibg1=ibg1+dtset%nspinor*nband_me
764 844837 : icg1=icg1+npw1_k*dtset%nspinor*nband_me
765 844837 : ikg1=ikg1+npw1_k
766 : end if
767 :
768 : end do !ikpt loop
769 :
770 : !======================================================================
771 : !================== END BIG K POINT LOOP ============================
772 : !======================================================================
773 :
774 20363 : call timab(125,2,tsec)
775 :
776 : ! Transfer density on augmented fft grid to normal fft grid in real space. Also take into account the spin.
777 : ! FR EB for the non-collinear part see vtorho.F90
778 39139 : if(iscf_mod>0) then
779 16513 : if (psps%usepaw==0) then
780 13465 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,1),1)
781 13465 : if(nspden==4)then
782 1012 : do ispden=2,4
783 1012 : call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,ispden),1)
784 : end do
785 : end if
786 : else
787 3048 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,1),1)
788 3048 : if(nspden==4)then
789 0 : do ispden=2,4
790 0 : call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,ispden),1)
791 : end do
792 : end if
793 : end if
794 : end if
795 :
796 : end do ! End loop over spins
797 :
798 : !More memory cleaning
799 18776 : call gs_hamkq%free()
800 18776 : call rf_hamkq%free()
801 18776 : if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) call rf_hamk_dir2%free()
802 :
803 18776 : ABI_FREE(rhoaug1)
804 18776 : ABI_FREE(vlocal)
805 18776 : ABI_FREE(vlocal1)
806 18776 : ABI_SFREE(vxctaulocal)
807 18776 : ABI_SFREE(vectornd_pac)
808 18776 : ABI_SFREE(vectornd_pac_idir)
809 :
810 18776 : call timab(124,2,tsec)
811 :
812 : !=== MPI communications ==================
813 : if(xmpi_paral==1)then
814 18776 : call timab(129,1,tsec)
815 :
816 : ! MG: For the record, buffer1 can be pretty big if we have dense k-meshes e.g. metals
817 : ! this is what you get with nband 26 and ngkpt 42**3:
818 : !
819 : ! [0] <var=buffer1, A@m_dfpt_vtorho.F90:773, addr=0x150d5a59e010, size_mb=401.130>
820 : !
821 : ! and this can lead to OOM if we have 2Gb per core also because xmpi_sum allocates another array of the same size!
822 : ! TODO: Avoid packing rhor1 in buffer
823 :
824 : ! Compute buffer size
825 18776 : buffer_size=11
826 18776 : ABI_MALLOC(buffer1,(buffer_size))
827 :
828 : ! Pack edocc,eeig0,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,eigen1,resid
829 18776 : index1=0
830 18776 : buffer1(index1+1)=edocc;buffer1(index1+2)=eeig0
831 18776 : buffer1(index1+3)=ek0 ;buffer1(index1+4)=ek1
832 18776 : buffer1(index1+5)=eloc0;buffer1(index1+6)=enl0
833 18776 : buffer1(index1+7)=enl1
834 18776 : buffer1(index1+8)=end0;buffer1(index1+9)=end1
835 18776 : buffer1(index1+10)=evxctau0;buffer1(index1+11)=evxctau1
836 18776 : index1=index1+11
837 : if(index1<buffer_size)buffer1(index1+1:buffer_size)=zero
838 :
839 : ! Build sum of everything
840 18776 : call timab(48,1,tsec)
841 18776 : call xmpi_sum(buffer1,buffer_size,spaceworld,ierr)
842 18776 : call timab(48,2,tsec)
843 :
844 : ! Unpack the final result
845 18776 : index1=0
846 18776 : edocc=buffer1(index1+1);eeig0=buffer1(index1+2)
847 18776 : ek0=buffer1(index1+3) ;ek1=buffer1(index1+4)
848 18776 : eloc0=buffer1(index1+5);enl0=buffer1(index1+6)
849 18776 : enl1=buffer1(index1+7)
850 18776 : end0=buffer1(index1+8);end1=buffer1(index1+9)
851 18776 : evxctau0=buffer1(index1+10);evxctau1=buffer1(index1+11)
852 18776 : ABI_FREE(buffer1)
853 :
854 : ! sync eigen1 and resid as well. No need to pack in buffer1, it duplicates a huge chunk of memory
855 18776 : call timab(48,1,tsec)
856 18776 : buffer_size=mbd2kpsp
857 18776 : call xmpi_sum(eigen1,buffer_size,spaceworld,ierr)
858 18776 : call timab(48,2,tsec)
859 :
860 18776 : call timab(48,1,tsec)
861 18776 : buffer_size=mbdkpsp
862 18776 : call xmpi_sum(resid,buffer_size,spaceworld,ierr)
863 18776 : call timab(48,2,tsec)
864 :
865 : ! sync rhor1 or rho1wfr as well. No need to pack in buffer1, it duplicates a huge chunk of memory
866 18776 : if(iscf_mod>0) then
867 15118 : buffer_size = cplex*dtset%nfft
868 15118 : if (psps%usepaw==0) then
869 12471 : call timab(48,1,tsec)
870 : ! TODO: add a pointer or something to get this into a single call to xmpisum
871 26695 : do isppol=1,nspden
872 26695 : call xmpi_sum(rhor1(:,isppol),buffer_size,spaceworld,ierr)
873 : end do
874 12471 : call timab(48,2,tsec)
875 : else
876 2647 : call timab(48,1,tsec)
877 5695 : do isppol=1,nspden
878 5695 : call xmpi_sum(rho1wfr(:,isppol),buffer_size,spaceworld,ierr)
879 : end do
880 2647 : call timab(48,2,tsec)
881 : end if
882 : end if
883 :
884 : ! Accumulate PAW occupancies
885 18776 : if (psps%usepaw==1.and.iscf_mod>0) then
886 2647 : call pawrhoij_mpisum_unpacked(pawrhoij1_unsym,spaceworld)
887 : end if
888 :
889 37552 : call timab(129,2,tsec)
890 : end if ! if kpt parallel
891 :
892 18776 : call timab(127,1,tsec)
893 :
894 : !If needed, compute rhog1, and symmetrize the density
895 18776 : if (iscf_mod > 0) then
896 :
897 : ! In order to have the symrhg working in parallel on FFT coefficients, the size
898 : ! of irzzon1 and phnons1 should be set to nfftot. Therefore, nsym\=1 does not work.
899 :
900 15118 : if(nspden==4) then
901 : ! FR symrhg will manage correctly this rearrangement
902 2884087 : rhor1(:,2)=rhor1(:,2)+(rhor1(:,1)+rhor1(:,4)) ! (n+mx)
903 2884087 : rhor1(:,3)=rhor1(:,3)+(rhor1(:,1)+rhor1(:,4)) ! (n+my)
904 : end if
905 : !
906 15118 : if (psps%usepaw==0) then
907 : call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,&
908 12471 : nspden,nsppol,nsym1,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,tnons1)
909 : else
910 : call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,&
911 2647 : nspden,nsppol,nsym1,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,tnons1)
912 : end if
913 : ! We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
914 : ! we also have the spin-up density, symmetrized, in rhor1(:,2).
915 : end if
916 :
917 18776 : ABI_FREE(grad_berry)
918 :
919 : !Find largest residual over bands, k points, and spins except for nbdbuf highest bands
920 18776 : ibdkpt=1
921 18776 : residm=zero
922 39139 : do isppol=1,nsppol
923 1064596 : do ikpt=1,nkpt_rbz
924 1025457 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
925 1025457 : nband_eff=max(1,nband_k-dtset%nbdbuf)
926 6752804 : residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
927 1045820 : ibdkpt=ibdkpt+nband_k
928 : end do
929 : end do
930 :
931 18776 : call timab(127,2,tsec)
932 :
933 18776 : if (iscf_mod>0) then
934 :
935 : ! PAW: Build new 1st-order rhoij quantities then symetrize them
936 : ! Compute and add the 1st-order compensation density to rho1wfr
937 : ! to get the total 1st-order density
938 15118 : if (psps%usepaw==1) then
939 : call pawmkrho(1,arg,cplex,gprimd,idir,indsy1,ipert,mpi_enreg,&
940 : my_natom,natom,nspden,nsym1,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
941 : dtset%pawprtvol,pawrhoij1,pawrhoij1_unsym,pawtab,dtset%qptn,rho1wfg,rho1wfr,&
942 : rhor1,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,&
943 2647 : pawang_sym=pawang1,pawnhat=nhat1,pawrhoij0=pawrhoij,rhog=rhog1)
944 2647 : ABI_FREE(rho1wfr)
945 2647 : ABI_FREE(rho1wfg)
946 2647 : if (paral_atom) then
947 1254 : call pawrhoij_free(pawrhoij1_unsym)
948 1254 : ABI_FREE(pawrhoij1_unsym)
949 : end if
950 : end if
951 :
952 : ! Compute density residual (if required) and its squared norm
953 15118 : if (optres==1) then
954 0 : nvresid1=rhor1-nvresid1
955 0 : call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid1)
956 : end if
957 : end if ! iscf>0
958 :
959 18776 : call timab(121,2,tsec)
960 : ABI_NVTX_END_RANGE()
961 :
962 : DBG_EXIT('COLL')
963 :
964 37552 : end subroutine dfpt_vtorho
965 : !!***
966 :
967 : end module m_dfpt_vtorho
968 : !!***
|
