Branch data 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 : 18888 : subroutine dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dbl_nnsclo,&
225 : 18888 : & dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,qphon,&
226 : 18888 : & edocc,eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,end0,end1,enl0,enl1,&
227 : 18888 : & evxctau0,evxctau1,fermie1,gh0c1_set,gh1c_set,gmet,gprimd,idir,indsy1,&
228 : 18888 : & ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,mband_mem,&
229 : : & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,&
230 : 18888 : & natom,nband_rbz,ncpgr,nfftf,nhat1,nkpt_rbz,npwarr,npwar1,nres2,nspden,&
231 : 18888 : & nsppol,nsym1,ntypat,nvresid1,occkq,occ_rbz,optres,&
232 : 18888 : & paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,pawrhoij1,pawtab,&
233 : 18888 : & phnons1,ph1d,prtvol,psps,pwindall,qmat,resid,residm,rhog1,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,ucvol,&
234 : 18888 : & 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 : 18888 : type(gs_hamiltonian_type) :: gs_hamkq
316 : 18888 : type(rf_hamiltonian_type) :: rf_hamkq,rf_hamk_dir2
317 : : !arrays
318 : 18888 : integer,allocatable :: kg1_k(:,:),kg_k(:,:)
319 : 18888 : integer, pointer :: my_atmtab(:)
320 : : real(dp) :: kpoint(3),kpq(3)
321 : : real(dp) :: tsec(2)
322 : 18888 : real(dp),allocatable :: buffer1(:)
323 : 18888 : real(dp),allocatable :: ddkinpw(:),dkinpw(:),dkinpw2(:)
324 : 18888 : real(dp),allocatable :: doccde_k(:),doccde_kq(:)
325 : 18888 : real(dp),allocatable :: edocc_k(:),eeig0_k(:),eig0_k(:),eig0_kq(:),eig1_k(:)
326 : 18888 : real(dp),allocatable :: ek0_k(:),ek1_k(:),eloc0_k(:),end0_k(:),end1_k(:),enl0_k(:),enl1_k(:)
327 : 18888 : real(dp),allocatable :: evxctau0_k(:),evxctau1_k(:)
328 : 18888 : real(dp),allocatable :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:),ffnlk(:,:,:,:)
329 : 18888 : real(dp),allocatable :: grad_berry(:,:,:),kinpw1(:),kpg1_k(:,:)
330 : 18888 : real(dp),allocatable :: kpg_k(:,:),occ_k(:),occ_kq(:)
331 : 18888 : real(dp),allocatable :: ph3d(:,:,:),ph3d1(:,:,:),resid_k(:)
332 : 18888 : real(dp),allocatable :: rho1wfg(:,:),rho1wfr(:,:),rhoaug1(:,:,:,:),rocceig(:,:)
333 : 18888 : real(dp),allocatable :: vectornd_pac(:,:,:,:,:),vectornd_pac_idir(:,:,:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:)
334 : 18888 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
335 : 18888 : real(dp),allocatable :: ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:)
336 : 18888 : 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 : 18888 : call timab(121,1,tsec)
345 : 18888 : call timab(124,1,tsec)
346 : :
347 : : !Retrieve parallelism data
348 : 18888 : spaceworld=mpi_enreg%comm_cell
349 : 18888 : me=mpi_enreg%me_kpt
350 : 18888 : paral_atom=(my_natom/=natom)
351 : 18888 : my_atmtab=>mpi_enreg%my_atmtab
352 : :
353 : : if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
354 [ + + ]: 18888 : & 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 [ - + ]: 18791 : 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 [ + + + - : 18888 : 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 [ + + + + : 18888 : iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3
+ + ]
367 : :
368 : 18888 : edocc=zero ; eeig0=zero ; ek0=zero ; ek1=zero
369 : 18888 : eloc0=zero ; end0=zero ; end1=zero ; enl0=zero ; enl1=zero
370 : 18888 : evxctau0=zero; evxctau1=zero
371 : 18888 : bdtot_index=0
372 : 18888 : bd2tot_index=0
373 : 18888 : ibg=0;icg=0
374 : 18888 : ibgq=0;icgq=0
375 : 18888 : ibg1=0;icg1=0
376 : 18888 : mbdkpsp=mband*nkpt_rbz*nsppol
377 : 18888 : mbd2kpsp=2*mband**2*nkpt_rbz*nsppol
378 : :
379 : 18888 : n1=dtset%ngfft(1); n2=dtset%ngfft(2); n3=dtset%ngfft(3)
380 : 18888 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6)
381 : 18888 : qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=tol14)
382 : :
383 [ + + ]: 18888 : omega_=zero; if(present(omega))omega_=omega
384 [ + + ]: 18888 : 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 [ + + ]: 18888 : if (iscf_mod>0) then
389 [ - + - - : 15158 : if (optres==1) nvresid1=rhor1
- - ]
390 [ + + ]: 15158 : if (psps%usepaw==0) then
391 [ + + + + ]: 160367912 : rhor1(:,:)=zero
392 : : else
393 [ - + + - : 10748 : ABI_MALLOC(rho1wfr,(cplex*dtset%nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
394 [ - + + - : 8061 : ABI_MALLOC(rho1wfg,(2,dtset%nfft))
+ - + - -
+ - + ]
395 [ + + + + ]: 27774285 : 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 [ + + ]: 18888 : if(iscf_mod<=0 .and. iscf_mod/=-3)then
401 : 124 : nnsclo_now=dtset%nstep
402 : : else
403 [ - + ]: 18764 : if(dtset%nnsclo>0)then
404 : 0 : nnsclo_now=dtset%nnsclo
405 : : else
406 : 18764 : nnsclo_now=1
407 : : end if
408 [ - + ]: 18764 : if(dbl_nnsclo==1) nnsclo_now=nnsclo_now*2
409 : : end if
410 : :
411 : : !Prepare GS k+q wf
412 : 18888 : mcgq=mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol;mcgq_disk=0
413 : :
414 : : !Prepare RF PAW files
415 [ + + ]: 18888 : if (psps%usepaw==1) then
416 : 3608 : 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 : 18888 : 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 : 18888 : 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 : 18888 : & 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 : 18888 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
443 [ + + + + : 18888 : 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 : 18888 : pawrhoij1_unsym => pawrhoij1
450 [ + + + + ]: 18888 : if (psps%usepaw==1.and.iscf_mod>0) then
451 [ + + ]: 2687 : 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 : 1433 : pawrhoij1_unsym => pawrhoij1
460 : 1433 : call pawrhoij_init_unpacked(pawrhoij1_unsym)
461 : : end if
462 : : end if
463 : :
464 [ - + + - : 113328 : ABI_MALLOC(rhoaug1,(cplex*n4,n5,n6,gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
465 [ - + + - : 113328 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
466 [ - + + - : 94440 : ABI_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - - + -
+ - + ]
467 : :
468 : 18888 : with_vxctau = ( usevxctau > 0 )
469 [ + + ]: 18888 : if(with_vxctau) then
470 [ - + + - : 294 : ABI_MALLOC(vxctaulocal,(n4,n5,n6,gs_hamkq%nvloc,4))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + - + ]
471 : : end if
472 : :
473 : 18888 : has_vectornd = (with_vectornd .EQ. 1)
474 [ + + ]: 18888 : if(has_vectornd) then
475 [ - + + - : 3066 : ABI_MALLOC(vectornd_pac,(n4,n5,n6,gs_hamkq%nvloc,3))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + - + ]
476 [ - + + - : 2190 : ABI_MALLOC(vectornd_pac_idir,(n4,n5,n6,gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - - + -
+ - + ]
477 : : end if
478 : :
479 : 18888 : nlines_done = 0
480 [ + + ]: 5425496 : resid = zero
481 : :
482 : : !LOOP OVER SPINS
483 [ + + ]: 39403 : do isppol=1,nsppol
484 : :
485 : : ! Rewind kpgsph data file if needed:
486 : 20515 : 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 : 20515 : gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1)
493 : :
494 : : ! Continue to initialize the Hamiltonian
495 : 20515 : call gs_hamkq%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.)
496 : 20515 : call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.)
497 [ + + + + : 20515 : 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 [ + + + + : 20515 : 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 [ + + + + : 277379659 : 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 [ + + ]: 20515 : if(has_vectornd) then
515 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
516 : 462 : & dtset%nspden, gs_hamkq%nvloc, 3, pawfgr, mpi_enreg, vectornd, vectornd_pac)
517 : 462 : call gs_hamkq%load_spin(isppol, vectornd=vectornd_pac)
518 [ + + + + : 40282665 : vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,idir)
+ + + + ]
519 : 462 : 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 [ + + ]: 20515 : 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 : 20515 : call timab(125,1,tsec)
531 : :
532 : : !======================================================================
533 : : !============== BIG FAT K POINT LOOP ================================
534 : : !======================================================================
535 : :
536 [ + + ]: 1046524 : do ikpt=1,nkpt_rbz
537 : 1026009 : counter=100*ikpt+isppol
538 : :
539 : 1026009 : nband_k = nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
540 : : ! enables variable nband less than the block size in serial case
541 : 1026009 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
542 : 1026009 : istwf_k = istwfk_rbz(ikpt)
543 : 1026009 : npw_k = npwarr(ikpt,1)
544 : 1026009 : npw1_k = npwar1(ikpt,1)
545 : 1026009 : wtk_k = wtk_rbz(ikpt)
546 : :
547 [ + + ]: 1026009 : 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 [ + + ]: 3381556 : kpoint(:)=kpt_rbz(:,ikpt)
557 : 845389 : kpq(:)=kpoint(:)
558 [ + + + + : 845389 : if (ipert<natom+3.or.ipert==natom+5.or.ipert==natom+6.or.&
+ + + + +
- ]
559 [ + + ]: 3107008 : &(ipert>natom+11.and.ipert<=2*natom+11)) kpq(:)=kpq(:)+qphon(1:3)
560 [ - + + - : 2536167 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + -
+ ]
561 [ - + + - : 2536167 : ABI_MALLOC(kg1_k,(3,npw1_k))
+ - + - -
+ - + -
+ ]
562 [ - + + - : 3381556 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + ]
563 [ - + + - : 3381556 : ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + ]
564 [ - + + - : 4226945 : ABI_MALLOC(ylmgr1_k,(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
+ - - + +
- + + + -
+ - + + -
+ - + ]
565 [ - + + - : 2536167 : ABI_MALLOC(doccde_k,(nband_k))
+ - + - -
+ - + ]
566 [ - + + - : 1690778 : ABI_MALLOC(doccde_kq,(nband_k))
- + - + ]
567 [ - + + - : 1690778 : ABI_MALLOC(eig0_k,(nband_k))
- + - + ]
568 [ - + + - : 1690778 : ABI_MALLOC(eig0_kq,(nband_k))
- + - + ]
569 [ + - - + : 2536167 : ABI_MALLOC(eig1_k,(2*nband_k**2))
+ - + - +
- - + -
+ ]
570 [ - + + - : 1690778 : ABI_MALLOC(edocc_k,(nband_k))
- + - + ]
571 [ - + + - : 1690778 : ABI_MALLOC(eeig0_k,(nband_k))
- + - + ]
572 [ - + + - : 1690778 : ABI_MALLOC(ek0_k,(nband_k))
- + - + ]
573 [ - + + - : 1690778 : ABI_MALLOC(ek1_k,(nband_k))
- + - + ]
574 [ - + + - : 1690778 : ABI_MALLOC(eloc0_k,(nband_k))
- + - + ]
575 [ - + + - : 1690778 : ABI_MALLOC(end0_k,(nband_k))
- + - + ]
576 [ - + + - : 1690778 : ABI_MALLOC(end1_k,(nband_k))
- + - + ]
577 [ - + + - : 1690778 : ABI_MALLOC(enl0_k,(nband_k))
- + - + ]
578 [ - + + - : 1690778 : ABI_MALLOC(enl1_k,(nband_k))
- + - + ]
579 [ - + + - : 1690778 : ABI_MALLOC(evxctau0_k,(nband_k))
- + - + ]
580 [ - + + - : 1690778 : ABI_MALLOC(evxctau1_k,(nband_k))
- + - + ]
581 [ - + + - : 1690778 : ABI_MALLOC(occ_k,(nband_k))
- + - + ]
582 [ - + + - : 1690778 : ABI_MALLOC(occ_kq,(nband_k))
- + - + ]
583 [ - + + - : 1690778 : ABI_MALLOC(resid_k,(nband_k))
- + - + ]
584 [ - + + - : 3381556 : ABI_MALLOC(rocceig,(nband_k,nband_k))
- + + - +
- + - - +
- + ]
585 : :
586 [ + + ]: 70401985 : eig1_k(:)=zero
587 [ + + ]: 5296089 : eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
588 [ + + ]: 5296089 : eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
589 [ + + ]: 5296089 : edocc_k(:)=zero
590 [ + + + + : 14197489 : eeig0_k(:)=zero ; ek0_k(:)=zero ; ek1_k(:)=zero
+ + ]
591 [ + + + + : 14197489 : eloc0_k(:)=zero ; end0_k(:)=zero ; end1_k(:)=zero
+ + ]
592 [ + + + + ]: 9746789 : enl0_k(:)=zero ; enl1_k(:)=zero
593 [ + + + + ]: 9746789 : evxctau0_k(:)=zero; evxctau1_k(:)=zero
594 [ + + ]: 5296089 : occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
595 [ + + ]: 5296089 : occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index)
596 [ + + ]: 5296089 : doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index)
597 [ + + ]: 5296089 : doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index)
598 [ + + ]: 5296089 : 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 : 845389 : 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 [ - + ]: 845389 : ABI_FREE(doccde_k)
607 [ - + ]: 845389 : ABI_FREE(doccde_kq)
608 [ - + ]: 845389 : ABI_FREE(occ_kq)
609 : :
610 [ + + + + ]: 477810477 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
611 [ + + ]: 845389 : if (psps%useylm==1) then
612 [ + + ]: 2617916 : do ilm=1,psps%mpsang*psps%mpsang
613 [ + + ]: 324775045 : 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 [ + + + + ]: 477843797 : kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
619 [ + + ]: 845389 : if (psps%useylm==1) then
620 [ + + ]: 2617916 : do ilm=1,psps%mpsang*psps%mpsang
621 [ + + ]: 324775045 : ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
622 : : end do
623 [ + + ]: 268562 : if (useylmgr1==1) then
624 [ + + ]: 782746 : do ilm=1,psps%mpsang*psps%mpsang
625 [ + + ]: 3662242 : do ii=1,3+6*((ipert-natom)/10)
626 [ + + ]: 413706564 : 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 [ - + - + : 845389 : ffnl1_test=ffnl1_test) ! Optional
- + - + -
+ - + -
+ ]
639 : :
640 [ - + ]: 845389 : 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 [ - + ]: 845389 : 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 [ - + ]: 845389 : ABI_FREE(ylm_k)
669 [ - + ]: 845389 : ABI_FREE(ylm1_k)
670 [ - + ]: 845389 : 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 : 845389 : 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 [ - + ]: 845389 : & rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,ddk_f,wtk_k,nlines_done,cg1_out,eta=eta_,omega=omega_)
684 : :
685 [ - + ]: 845389 : 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 [ - + ]: 845389 : ABI_FREE(kinpw1)
691 [ - + ]: 845389 : ABI_FREE(kg_k)
692 [ - + ]: 845389 : ABI_FREE(kg1_k)
693 [ - + ]: 845389 : ABI_FREE(kpg_k)
694 [ - + ]: 845389 : ABI_FREE(kpg1_k)
695 [ - + ]: 845389 : ABI_FREE(dkinpw)
696 [ + + ]: 845389 : 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 [ - + ]: 845389 : ABI_FREE(ffnlk)
703 [ - + ]: 845389 : ABI_FREE(ffnl1)
704 [ + + ]: 845389 : ABI_SFREE(ffnl1_test)
705 [ - + ]: 845389 : ABI_FREE(eig0_k)
706 [ - + ]: 845389 : ABI_FREE(eig0_kq)
707 [ - + ]: 845389 : ABI_FREE(rocceig)
708 [ - + ]: 845389 : ABI_FREE(ph3d)
709 [ + + ]: 845389 : ABI_SFREE(ph3d1)
710 : :
711 : : ! Save eigenvalues (hartree), residuals (hartree**2)
712 [ + + ]: 70401985 : eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:)
713 [ + + ]: 5296089 : 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 [ + + + - ]: 845389 : if (iscf_mod>0 .or. iscf_mod==-3 .or. iscf_mod==-2)then
718 [ + + ]: 5296089 : do iband=1,nband_k
719 : 4450700 : edocc=edocc+wtk_k*occ_k(iband)*edocc_k(iband)
720 : 4450700 : eeig0=eeig0+wtk_k*occ_k(iband)*eeig0_k(iband)
721 : 4450700 : ek0=ek0+wtk_k*occ_k(iband)*ek0_k(iband)
722 : 4450700 : ek1=ek1+wtk_k*occ_k(iband)*ek1_k(iband)
723 : 4450700 : eloc0=eloc0+wtk_k*occ_k(iband)*eloc0_k(iband)
724 : 4450700 : end0=end0+wtk_k*occ_k(iband)*end0_k(iband)
725 : 4450700 : end1=end1+wtk_k*occ_k(iband)*end1_k(iband)
726 : 4450700 : enl0=enl0+wtk_k*occ_k(iband)*enl0_k(iband)
727 : 4450700 : enl1=enl1+wtk_k*occ_k(iband)*enl1_k(iband)
728 : 4450700 : evxctau0=evxctau0+wtk_k*occ_k(iband)*evxctau0_k(iband)
729 : 5296089 : evxctau1=evxctau1+wtk_k*occ_k(iband)*evxctau1_k(iband)
730 : : end do
731 : : end if
732 : :
733 [ - + ]: 845389 : ABI_FREE(eig1_k)
734 [ - + ]: 845389 : ABI_FREE(occ_k)
735 [ - + ]: 845389 : ABI_FREE(resid_k)
736 [ - + ]: 845389 : ABI_FREE(edocc_k)
737 [ - + ]: 845389 : ABI_FREE(eeig0_k)
738 [ - + ]: 845389 : ABI_FREE(ek0_k)
739 [ - + ]: 845389 : ABI_FREE(ek1_k)
740 [ - + ]: 845389 : ABI_FREE(eloc0_k)
741 [ - + ]: 845389 : ABI_FREE(end0_k)
742 [ - + ]: 845389 : ABI_FREE(end1_k)
743 [ - + ]: 845389 : ABI_FREE(enl0_k)
744 [ - + ]: 845389 : ABI_FREE(enl1_k)
745 [ - + ]: 845389 : ABI_FREE(evxctau0_k)
746 [ - + ]: 845389 : 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 : 845389 : bdtot_index=bdtot_index+nband_k
750 : 845389 : bd2tot_index=bd2tot_index+2*nband_k**2
751 : :
752 : : ! Shift array memory
753 [ + - ]: 845389 : if (mkmem/=0) then
754 : 845389 : ibg=ibg+dtset%nspinor*nband_me
755 : 845389 : icg=icg+npw_k*dtset%nspinor*nband_me
756 : 845389 : ikg=ikg+npw_k
757 : : end if
758 [ + - ]: 845389 : if (mkqmem/=0) then
759 : 845389 : ibgq=ibgq+dtset%nspinor*nband_me
760 : 845389 : icgq=icgq+npw1_k*dtset%nspinor*nband_me
761 : : end if
762 [ + - ]: 865904 : if (mk1mem/=0) then
763 : 845389 : ibg1=ibg1+dtset%nspinor*nband_me
764 : 845389 : icg1=icg1+npw1_k*dtset%nspinor*nband_me
765 : 845389 : ikg1=ikg1+npw1_k
766 : : end if
767 : :
768 : : end do !ikpt loop
769 : :
770 : : !======================================================================
771 : : !================== END BIG K POINT LOOP ============================
772 : : !======================================================================
773 : :
774 : 20515 : 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 [ + + ]: 39403 : if(iscf_mod>0) then
779 [ + + ]: 16593 : 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 : 3128 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,1),1)
788 [ - + ]: 3128 : 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 : 18888 : call gs_hamkq%free()
800 : 18888 : call rf_hamkq%free()
801 [ + + + + : 18888 : if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) call rf_hamk_dir2%free()
+ + ]
802 : :
803 [ - + ]: 18888 : ABI_FREE(rhoaug1)
804 [ - + ]: 18888 : ABI_FREE(vlocal)
805 [ - + ]: 18888 : ABI_FREE(vlocal1)
806 [ + + ]: 18888 : ABI_SFREE(vxctaulocal)
807 [ + + ]: 18888 : ABI_SFREE(vectornd_pac)
808 [ + + ]: 18888 : ABI_SFREE(vectornd_pac_idir)
809 : :
810 : 18888 : call timab(124,2,tsec)
811 : :
812 : : !=== MPI communications ==================
813 : : if(xmpi_paral==1)then
814 : 18888 : 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 : 18888 : buffer_size=11
826 [ - + - + ]: 18888 : ABI_MALLOC(buffer1,(buffer_size))
827 : :
828 : : ! Pack edocc,eeig0,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,eigen1,resid
829 : 18888 : index1=0
830 : 18888 : buffer1(index1+1)=edocc;buffer1(index1+2)=eeig0
831 : 18888 : buffer1(index1+3)=ek0 ;buffer1(index1+4)=ek1
832 : 18888 : buffer1(index1+5)=eloc0;buffer1(index1+6)=enl0
833 : 18888 : buffer1(index1+7)=enl1
834 : 18888 : buffer1(index1+8)=end0;buffer1(index1+9)=end1
835 : 18888 : buffer1(index1+10)=evxctau0;buffer1(index1+11)=evxctau1
836 : 18888 : index1=index1+11
837 : : if(index1<buffer_size)buffer1(index1+1:buffer_size)=zero
838 : :
839 : : ! Build sum of everything
840 : 18888 : call timab(48,1,tsec)
841 : 18888 : call xmpi_sum(buffer1,buffer_size,spaceworld,ierr)
842 : 18888 : call timab(48,2,tsec)
843 : :
844 : : ! Unpack the final result
845 : 18888 : index1=0
846 : 18888 : edocc=buffer1(index1+1);eeig0=buffer1(index1+2)
847 : 18888 : ek0=buffer1(index1+3) ;ek1=buffer1(index1+4)
848 : 18888 : eloc0=buffer1(index1+5);enl0=buffer1(index1+6)
849 : 18888 : enl1=buffer1(index1+7)
850 : 18888 : end0=buffer1(index1+8);end1=buffer1(index1+9)
851 : 18888 : evxctau0=buffer1(index1+10);evxctau1=buffer1(index1+11)
852 : 18888 : 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 : 18888 : call timab(48,1,tsec)
856 : 18888 : buffer_size=mbd2kpsp
857 : 18888 : call xmpi_sum(eigen1,buffer_size,spaceworld,ierr)
858 : 18888 : call timab(48,2,tsec)
859 : :
860 : 18888 : call timab(48,1,tsec)
861 : 18888 : buffer_size=mbdkpsp
862 : 18888 : call xmpi_sum(resid,buffer_size,spaceworld,ierr)
863 : 18888 : 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 [ + + ]: 18888 : if(iscf_mod>0) then
867 : 15158 : buffer_size = cplex*dtset%nfft
868 [ + + ]: 15158 : 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 : 2687 : call timab(48,1,tsec)
877 [ + + ]: 5815 : do isppol=1,nspden
878 : 5815 : call xmpi_sum(rho1wfr(:,isppol),buffer_size,spaceworld,ierr)
879 : : end do
880 : 2687 : call timab(48,2,tsec)
881 : : end if
882 : : end if
883 : :
884 : : ! Accumulate PAW occupancies
885 [ + + + + ]: 18888 : if (psps%usepaw==1.and.iscf_mod>0) then
886 : 2687 : call pawrhoij_mpisum_unpacked(pawrhoij1_unsym,spaceworld)
887 : : end if
888 : :
889 : 37776 : call timab(129,2,tsec)
890 : : end if ! if kpt parallel
891 : :
892 : 18888 : call timab(127,1,tsec)
893 : :
894 : : !If needed, compute rhog1, and symmetrize the density
895 [ + + ]: 18888 : 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 [ + + ]: 15158 : 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 [ + + ]: 15158 : 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 : 2687 : 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 [ - + ]: 18888 : ABI_FREE(grad_berry)
918 : :
919 : : !Find largest residual over bands, k points, and spins except for nbdbuf highest bands
920 : 18888 : ibdkpt=1
921 : 18888 : residm=zero
922 [ + + ]: 39403 : do isppol=1,nsppol
923 [ + + ]: 1065412 : do ikpt=1,nkpt_rbz
924 : 1026009 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
925 : 1026009 : nband_eff=max(1,nband_k-dtset%nbdbuf)
926 [ + - - + : 6765716 : residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
- - + + ]
927 : 1046524 : ibdkpt=ibdkpt+nband_k
928 : : end do
929 : : end do
930 : :
931 : 18888 : call timab(127,2,tsec)
932 : :
933 [ + + ]: 18888 : 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 [ + + ]: 15158 : 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 : 2687 : pawang_sym=pawang1,pawnhat=nhat1,pawrhoij0=pawrhoij,rhog=rhog1)
944 [ - + ]: 2687 : ABI_FREE(rho1wfr)
945 [ - + ]: 2687 : ABI_FREE(rho1wfg)
946 [ + + ]: 2687 : 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 [ - + ]: 15158 : 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 : 18888 : call timab(121,2,tsec)
960 : : ABI_NVTX_END_RANGE()
961 : :
962 : : DBG_EXIT('COLL')
963 : :
964 [ - + - + : 37776 : end subroutine dfpt_vtorho
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
965 : : !!***
966 : :
967 : : end module m_dfpt_vtorho
968 : : !!***
|
