Line data Source code
1 : !!****m* ABINIT/m_dfpt_vtowfk
2 : !! NAME
3 : !! m_dfpt_vtowfk
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1999-2026 ABINIT group (XG, AR, DRH, MB, MVer,XW, MT, GKA)
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_vtowfk
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_errors
29 : use m_xmpi
30 : use m_mpinfo
31 : use m_cgtools
32 : use m_wfk
33 : use m_rf2
34 : use m_dtset
35 : use m_dtfil
36 : use m_pawcprj
37 :
38 : use defs_datatypes, only : pseudopotential_type
39 : use defs_abitypes, only : MPI_type
40 : use m_rf2_init, only : rf2_init
41 : use m_time, only : timab
42 : use m_pawrhoij, only : pawrhoij_type
43 : use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type, KPRIME_H_KPRIME
44 : use m_spacepar, only : meanvalue_g
45 : use m_dfpt_mkrho, only : dfpt_accrho
46 : use m_dfpt_cgwf, only : dfpt_cgwf, full_active_wf1
47 : use m_getghc, only : getgsc, getghc_nucdip, getghc_mGGA
48 : use m_getgh1c, only : getgh1ndc, getgh1c_mGGA
49 :
50 : #if defined(HAVE_GPU_MARKERS)
51 : use m_nvtx_data
52 : #endif
53 :
54 : implicit none
55 :
56 : private
57 : !!***
58 :
59 : public :: dfpt_vtowfk
60 : !!***
61 :
62 : contains
63 : !!***
64 :
65 : !!****f* ABINIT/dfpt_vtowfk
66 : !! NAME
67 : !! dfpt_vtowfk
68 : !!
69 : !! FUNCTION
70 : !! This routine compute the partial density at a given k-point,
71 : !! for a given spin-polarization, from a fixed potential (vlocal1).
72 : !!
73 : !! INPUTS
74 : !! cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions
75 : !! cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
76 : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
77 : !! cplex=1 if rhoaug1 is real, 2 if rhoaug1 is complex
78 : !TODO MJV: PAW mband_mem
79 : !! cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
80 : !! projected with non-local projectors: cprj=<p_i|Cnk>
81 : !! cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q>
82 : !! dim_eig2rf = dimension for the second order eigenvalues
83 : !! dtfil <type(datafiles_type)>=variables related to files
84 : !! dtset <type(dataset_type)>=all input variables for this dataset
85 : !! eig0_k(nband_k)=GS eigenvalues at k (hartree)
86 : !! eig0_kq(nband_k)=GS eigenvalues at k+Q (hartree)
87 : !! fermie1=derivative of fermi energy wrt (strain) perturbation
88 : !! grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
89 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
90 : !! ibg=shift to be applied on the location of data in the array cprj
91 : !! ibgq=shift to be applied on the location of data in the array cprjq
92 : !! ibg1=shift to be applied on the location of data in the array cprj1
93 : !! icg=shift to be applied on the location of data in the array cg
94 : !! icgq=shift to be applied on the location of data in the array cgq
95 : !! icg1=shift to be applied on the location of data in the array cg1
96 : !! idir=direction of the current perturbation
97 : !! ikpt=k-point index number
98 : !! ipert=type of the perturbation
99 : !! isppol=1 index of current spin component
100 : !! mband=maximum number of bands
101 : !! mband_mem=maximum number of bands on this cpu
102 : !! mcgq=second dimension of the cgq array
103 : !! mcprjq=second dimension of the cprjq array
104 : !! mkmem =number of k points trated by this node (GS data).
105 : !! mk1mem =number of k points treated by this node (RF data)
106 : !! mpi_enreg=information about MPI parallelization
107 : !! mpw=maximum dimensioned size of npw or wfs at k
108 : !! mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
109 : !! natom=number of atoms in cell.
110 : !! nband_k=number of bands at this k point for that spin polarization
111 : !! ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
112 : !! nnsclo_now=number of non-self-consistent loops for the current vtrial
113 : !! (often 1 for SCF calculation, =nstep for non-SCF calculations)
114 : !! npw_k=number of plane waves at this k point
115 : !! npw1_k=number of plane waves at this k+q point
116 : !! nspinor=number of spinorial components of the wavefunctions
117 : !! nsppol=1 for unpolarized, 2 for spin-polarized
118 : !! n4,n5,n6 used for dimensioning real space arrays
119 : !! omega= frequency of the perturbation
120 : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
121 : !! prtvol=control print volume and debugging output
122 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
123 : !! rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
124 : !! rf_hamk_dir2 <type(rf_hamiltonian_type)>= (used only when ipert=natom+11, so q=0)
125 : !! same as rf_hamkq, but the direction of the perturbation is different
126 : !! rhoaug1(cplex*n4,n5,n6,nspden)= density in electrons/bohr**3,
127 : !! on the augmented fft grid. (cumulative, so input as well as output)
128 : !! rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
129 : !! if this ratio has been attributed to the band n (second argument), zero otherwise
130 : !! ddk<wfk_t>=struct info for DDK file.
131 : !! wtk_k=weight assigned to the k point.
132 : !!
133 : !! OUTPUT
134 : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF
135 : !! wavefunctions at k,q. They are orthogonalized to the occupied states.
136 : !! cg1_active(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)=pw coefficients of RF
137 : !! wavefunctions at k,q. They are orthogonalized to the active. Only needed for ieigrf/=0
138 : !! edocc_k(nband_k)=correction to 2nd-order total energy coming
139 : !! from changes of occupation
140 : !! eeig0_k(nband_k)=zero-order eigenvalues contribution to 2nd-order total
141 : !! energy from all bands at this k point.
142 : !! eig1_k(2*nband_k**2)=first-order eigenvalues (hartree)
143 : !! ek0_k(nband_k)=0-order kinetic energy contribution to 2nd-order total
144 : !! energy from all bands at this k point.
145 : !! ek1_k(nband_k)=1st-order kinetic energy contribution to 2nd-order total
146 : !! energy from all bands at this k point.
147 : !! eloc0_k(nband_k)=zero-order local contribution to 2nd-order total energy
148 : !! from all bands at this k point.
149 : !! end0_k(nband_k)=0-order nuclear dipole energy contribution to 2nd-order total
150 : !! energy from all bands at this k point.
151 : !! end1_k(nband_k)=1st-order nuclear dipole energy contribution to 2nd-order total
152 : !! energy from all bands at this k point.
153 : !! enl0_k(nband_k)=zero-order non-local contribution to 2nd-order total energy
154 : !! from all bands at this k point.
155 : !! enl1_k(nband_k)=first-order non-local contribution to 2nd-order total energy
156 : !! from all bands at this k point.
157 : !! evxctau0_k(nband_k)=0-order vxctau energy contribution to 2nd-order total
158 : !! energy from all bands at this k point.
159 : !! evxctau1_k(nband_k)=1-order vxctau energy contribution to 2nd-order total
160 : !! energy from all bands at this k point.
161 : !! gh1c_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
162 : !! gh0c1_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}k+q-eig^{(0)}nk|\Psi^{(1)}kq>
163 : !! The wavefunction is orthogonal to the active space (for metals). It is not coherent with cg1.
164 : !! resid_k(nband_k)=residuals for each band over all k points,
165 : !! rhoaug1(cplex*n4,n5,n6,nspden)= density in electrons/bohr**3,
166 : !! on the augmented fft grid. (cumulative, so input as well as output).
167 : !! ==== if (gs_hamkq%usepaw==1) ====
168 : !TODO MJV: PAW mband_mem
169 : !! cprj1(natom,nspinor*mband*mk1mem*nsppol*usecprj)=
170 : !! 1st-order wave functions at k,q projected with non-local projectors:
171 : !! cprj1=<p_i|C1nk,q> where p_i is a non-local projector
172 : !! pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
173 : !! (cumulative, so input as well as output)
174 : !!
175 : !! SOURCE
176 :
177 1689674 : subroutine dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,&
178 : & dim_eig2rf,dtfil,dtset,&
179 1689674 : & edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,&
180 844837 : & ek0_k,ek1_k,eloc0_k,end0_k,end1_k,enl0_k,enl1_k,evxctau0_k,evxctau1_k,&
181 844837 : & fermie1,ffnl1,ffnl1_test,gh0c1_set,gh1c_set,grad_berry,gs_hamkq,&
182 : & ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,&
183 : & isppol,mband,mband_mem,mcgq,mcprjq,mkmem,mk1mem,&
184 : & mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,&
185 : & nnsclo_now,npw_k,npw1_k,nspinor,nsppol,&
186 844837 : & n4,n5,n6,occ_k,pawrhoij1,prtvol,psps,resid_k,rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,&
187 : & ddk_f,wtk_k,nlines_done,cg1_out,eta,omega)
188 :
189 : !Arguments ------------------------------------
190 : !scalars
191 : integer,intent(in) :: cplex,dim_eig2rf,ibg
192 : integer,intent(in) :: ibg1,ibgq,icg,icg1,icgq,idir,ikpt,ipert,isppol
193 : integer,intent(in) :: mband,mcgq,mcprjq,mk1mem,mkmem
194 : integer,intent(in) :: mband_mem
195 : integer,intent(in) :: mpw,mpw1,n4,n5,n6,natom,ncpgr
196 : integer,intent(in) :: nnsclo_now,nspinor,nsppol,prtvol
197 : integer,optional,intent(in) :: cg1_out
198 : integer,intent(in) :: nband_k,npw1_k,npw_k
199 : integer,intent(inout) :: nlines_done
200 : real(dp),intent(in) :: fermie1,wtk_k
201 : real(dp),optional,intent(in) :: eta,omega
202 : type(MPI_type),intent(in) :: mpi_enreg
203 : type(datafiles_type),intent(in) :: dtfil
204 : type(dataset_type),intent(in) :: dtset
205 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
206 : type(rf_hamiltonian_type),intent(inout) :: rf_hamkq,rf_hamk_dir2
207 : type(pseudopotential_type),intent(in) :: psps
208 : !arrays
209 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol),cgq(2,mcgq)
210 : real(dp),intent(in) :: eig0_k(nband_k),eig0_kq(nband_k)
211 : real(dp),intent(in) :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
212 : real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband_k)
213 : real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k)
214 : real(dp),intent(inout) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)
215 : real(dp),intent(inout) :: rhoaug1(cplex*n4,n5,n6,gs_hamkq%nvloc)
216 : real(dp),intent(inout) :: cg1_active(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)
217 : real(dp),intent(inout) :: gh1c_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)
218 : real(dp),intent(inout) :: gh0c1_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)
219 : real(dp),intent(inout) :: edocc_k(nband_k),eeig0_k(nband_k),eig1_k(2*nband_k**2)
220 : real(dp),intent(out) :: ek0_k(nband_k),eloc0_k(nband_k)
221 : real(dp),intent(inout) :: ek1_k(nband_k)
222 : real(dp),intent(out) :: end0_k(nband_k),end1_k(nband_k),enl0_k(nband_k),enl1_k(nband_k)
223 : real(dp),intent(out) :: evxctau0_k(nband_k),evxctau1_k(nband_k)
224 : real(dp),intent(out) :: resid_k(nband_k)
225 : !TODO: PAW distrib bands mband_mem
226 : type(pawcprj_type),intent(in) :: cprj(natom,nspinor*mband_mem*mkmem*nsppol*gs_hamkq%usecprj)
227 : type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq)
228 : type(pawcprj_type),intent(inout) :: cprj1(natom,nspinor*mband_mem*mk1mem*nsppol*gs_hamkq%usecprj)
229 : type(pawrhoij_type),intent(inout) :: pawrhoij1(natom*gs_hamkq%usepaw)
230 : type(wfk_t),intent(inout) :: ddk_f(4)
231 :
232 : !Local variables-------------------------------
233 : !scalars
234 : integer,parameter :: level=14,tim_fourwf=5
235 : integer,save :: nskip=0
236 : integer :: iband,idir0,ierr,igs,igscq,ii,dim_dcwf,inonsc
237 : integer :: iband_me,nband_me !, unit_me
238 : integer :: iorder_cprj,iorder_cprj1,ipw,iscf_mod,ispinor,me,mgscq,nkpt_max
239 : integer :: option,opt_gvnlx1,quit,test_ddk,ndat
240 : integer :: tocceig,usedcwavef,ptr,shift_band
241 : real(dp) :: aa,ai,ar,eig0nk,eta_,omega_,resid,residk,scprod,energy_factor
242 : character(len=500) :: msg
243 844837 : type(rf2_t) :: rf2
244 : !arrays
245 844837 : logical,allocatable :: cycle_bands(:)
246 2534511 : integer :: rank_band(nband_k), bands_treated_now(nband_k)
247 : real(dp) :: tsec(2)
248 844837 : real(dp),allocatable :: cwave0(:,:),cwave1(:,:),cwavef(:,:)
249 844837 : real(dp),allocatable :: dcwavef(:,:),gh1c_n(:,:),gh0c1(:,:),ghc_vectornd(:,:)
250 844837 : real(dp),allocatable :: ghc_vxctau(:,:)
251 844837 : real(dp),allocatable :: gsc(:,:),gscq(:,:),gvnlx1(:,:),gvnlxc(:,:)
252 844837 : real(dp),pointer :: kinpw1(:)
253 844837 : type(pawcprj_type),allocatable :: cwaveprj(:,:),cwaveprj0(:,:),cwaveprj1(:,:)
254 : ! *********************************************************************
255 :
256 : DBG_ENTER('COLL')
257 :
258 : ABI_NVTX_START_RANGE(NVTX_DFPT_VTOWFK)
259 :
260 : !Keep track of total time spent in dfpt_vtowfk
261 844837 : call timab(128,1,tsec)
262 :
263 844837 : nkpt_max=50; if (xmpi_paral==1) nkpt_max=-1
264 :
265 844837 : if(prtvol>2 .or. ikpt<=nkpt_max)then
266 51182 : write(msg,'(2a,i5,2x,a,3f9.5,2x,a)')ch10,' Non-SCF iterations; k pt #',ikpt,'k=',gs_hamkq%kpt_k(:),'band residuals:'
267 51182 : call wrtout(std_out,msg)
268 : end if
269 :
270 : !Initializations and allocations
271 844837 : me=mpi_enreg%me_kpt
272 844837 : quit=0
273 :
274 : !The value of iscf must be modified if ddk perturbation
275 844837 : iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3
276 :
277 844837 : kinpw1 => gs_hamkq%kinpw_kp
278 2534511 : ABI_MALLOC(gh0c1,(2,npw1_k*nspinor))
279 1689674 : ABI_MALLOC(gvnlxc,(2,npw1_k*nspinor))
280 1689674 : ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor))
281 2534511 : ABI_MALLOC(cwave0,(2,npw_k*nspinor))
282 1689674 : ABI_MALLOC(cwavef,(2,npw1_k*nspinor))
283 1689674 : ABI_MALLOC(cwave1,(2,npw1_k*nspinor))
284 1689674 : ABI_MALLOC(gh1c_n,(2,npw1_k*nspinor))
285 844837 : if (gs_hamkq%usepaw==1) then
286 225530 : ABI_MALLOC(gsc,(2,npw1_k*nspinor))
287 : else
288 732072 : ABI_MALLOC(gsc,(0,0))
289 : end if
290 :
291 844837 : omega_=zero; if(present(omega))omega_=omega
292 844837 : eta_=zero; if(present(eta))eta_=eta
293 :
294 : !Read the npw and kg records of wf files
295 844837 : test_ddk=0
296 : if ((ipert==natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and.&
297 : & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
298 : & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.&
299 3379348 : & ipert==natom+10.or.ipert==natom+11) then
300 104202 : test_ddk=1
301 844837 : if(ipert==natom+10.or.ipert==natom+11) test_ddk=0
302 : end if
303 :
304 : !Additional stuff for PAW
305 844837 : ABI_MALLOC(cwaveprj0,(0,0))
306 844837 : if (gs_hamkq%usepaw==1) then
307 : ! 1-Compute all <g|S|Cnk+q>
308 112765 : igscq=0
309 : !TODO MJV: PAW mband_mem
310 112765 : mgscq=mpw1*nspinor*mband_mem
311 338295 : ABI_MALLOC_OR_DIE(gscq,(2,mgscq), ierr)
312 : #ifdef HAVE_OPENMP_OFFLOAD
313 : !$OMP TARGET ENTER DATA MAP(alloc:gscq) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
314 : #endif
315 :
316 : ABI_NVTX_START_RANGE(NVTX_GETGSC)
317 : call getgsc(cgq,cprjq,gs_hamkq,gscq,ibgq,icgq,igscq,ikpt,isppol,mcgq,mcprjq,&
318 112765 : & mgscq,mpi_enreg,dtset%bandpp,natom,nband_k,npw1_k,dtset%nspinor,select_k=KPRIME_H_KPRIME)
319 : ABI_NVTX_END_RANGE()
320 : ! 2-Initialize additional scalars/arrays
321 112765 : iorder_cprj=0;iorder_cprj1=0
322 112765 : dim_dcwf=npw1_k*nspinor;if (ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) dim_dcwf=0
323 338295 : ABI_MALLOC(dcwavef,(2,dim_dcwf))
324 112765 : if (gs_hamkq%usecprj==1) then
325 112765 : ABI_FREE(cwaveprj0)
326 804010 : ABI_MALLOC(cwaveprj0,(natom,nspinor))
327 112765 : call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj)
328 : end if
329 804010 : ABI_MALLOC(cwaveprj,(natom,nspinor))
330 691245 : ABI_MALLOC(cwaveprj1,(natom,nspinor))
331 112765 : call pawcprj_alloc(cwaveprj ,0,gs_hamkq%dimcprj)
332 112765 : call pawcprj_alloc(cwaveprj1,0,gs_hamkq%dimcprj)
333 : else
334 732072 : igscq=0;mgscq=0;dim_dcwf=0
335 732072 : ABI_MALLOC(gscq,(0,0))
336 732072 : ABI_MALLOC(dcwavef,(0,0))
337 732072 : ABI_MALLOC(cwaveprj,(0,0))
338 732072 : ABI_MALLOC(cwaveprj1,(0,0))
339 : end if
340 :
341 844837 : energy_factor=two
342 844837 : if(ipert==natom+10.or.ipert==natom+11) energy_factor=six
343 :
344 : !For rf2 perturbation :
345 844837 : if(ipert==natom+10.or.ipert==natom+11) then
346 : call rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,ffnl1,ffnl1_test,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,&
347 14616 : mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f)
348 : end if
349 :
350 844837 : call timab(139,1,tsec)
351 :
352 : !======================================================================
353 : !================== LOOP OVER BANDS ==================================
354 : !======================================================================
355 :
356 844837 : ndat=1;if (mpi_enreg%paral_kgb==1) ndat=mpi_enreg%bandpp
357 : call proc_distrb_band(rank_band,mpi_enreg%proc_distrb,ikpt,isppol,mband,&
358 844837 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
359 :
360 844837 : iband_me = 0
361 5280849 : do iband=1,nband_k
362 :
363 : ! Skip bands not treated by current proc
364 4436012 : if( (mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me)) cycle
365 4429964 : iband_me = iband_me + 1
366 :
367 : !unit_me = 300+iband
368 : !Get ground-state wavefunctions
369 4429964 : ptr = 1+(iband_me-1)*npw_k*nspinor+icg
370 4429964 : call cg_zcopy(npw_k*nspinor,cg(1,ptr),cwave0)
371 :
372 : ! Get PAW ground state projected WF (cprj)
373 4429964 : if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.ipert/=natom+10.and.ipert/=natom+11) then
374 589604 : idir0 = idir
375 589604 : if(ipert==natom+3.or.ipert==natom+4) idir0 =1
376 : ! PAW distributes cprj by band and k
377 : call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,natom,iband_me,ibg,ikpt,iorder_cprj,&
378 : & isppol,mband_mem,mkmem,natom,1,nband_me,nspinor,nsppol,dtfil%unpaw,&
379 : !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,&
380 589604 : & icpgr=idir0,ncpgr=ncpgr)
381 : end if
382 :
383 : ! Get first-order wavefunctions
384 4429964 : ptr = 1+(iband_me-1)*npw1_k*nspinor+icg1
385 4429964 : call cg_zcopy(npw1_k*nspinor,cg1(1,ptr),cwavef)
386 :
387 : ! Read PAW projected 1st-order WF (cprj)
388 : ! Unuseful for the time being (will be recomputed in dfpt_cgwf)
389 : ! if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
390 : !TODO MJV: PAW
391 : ! call pawcprj_get(gs_hamkq%atindx1,cwaveprj,cprj1,natom,iband,ibg1,ikpt,iorder_cprj1,&
392 : ! & isppol,mband,mk1mem,natom,1,nband_k,nspinor,nsppol,dtfil%unpaw1,
393 : ! & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
394 : ! end if
395 :
396 : ! Filter the wavefunctions for large modified kinetic energy
397 : ! The GS wavefunctions should already be non-zero
398 9003366 : do ispinor=1,nspinor
399 4573402 : igs=(ispinor-1)*npw1_k
400 849303984 : do ipw=1+igs,npw1_k+igs
401 844874020 : if(kinpw1(ipw-igs)>huge(zero)*1.d-20)then
402 0 : cwavef(1,ipw)=zero
403 0 : cwavef(2,ipw)=zero
404 : end if
405 : end do
406 : end do
407 :
408 :
409 : ! If electric field, the derivative of the wf should be read, and multiplied by i.
410 4429964 : if(test_ddk==1) then
411 414611 : ii = ddk_f(1)%findk(gs_hamkq%kpt_k)
412 414611 : ABI_CHECK(ii == ikpt, "ii != ikpt, something is wrong with k-point, check kptopt/ngkpt, etc")
413 : !TODO MJV: check if this iband should be _me
414 414611 : call ddk_f(1)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=gvnlx1)
415 :
416 : ! Multiplication by -i
417 : ! MVeithen 021212 : use + i instead,
418 : ! See X. Gonze, Phys. Rev. B 55, 10337 (1997) [[cite:Gonze1997]] Eq. (79)
419 : ! the operator used to compute the first-order derivative
420 : ! of the wavefunctions with respect to an electric field
421 : ! is $+i \frac{d}{dk}$
422 : ! This change will affect the computation of the 2dtes from non
423 : ! stationary expressions, see dfpt_nstdy.f and dfpt_nstwf.f
424 57915906 : do ipw=1,npw1_k*nspinor
425 : ! aa=gvnlx1(1,ipw)
426 : ! gvnlx1(1,ipw)=gvnlx1(2,ipw)
427 : ! gvnlx1(2,ipw)=-aa
428 57501295 : aa=gvnlx1(1,ipw)
429 57501295 : gvnlx1(1,ipw)=-gvnlx1(2,ipw)
430 57915906 : gvnlx1(2,ipw)=aa
431 : end do
432 : end if
433 :
434 : ! Unlike in GS calculations, the inonsc loop is inside the band loop
435 : ! nnsclo_now=number of non-self-consistent loops for the current vtrial
436 : ! (often 1 for SCF calculation, =nstep for non-SCF calculations)
437 8407022 : do inonsc=1,nnsclo_now
438 :
439 : ! Note that the following translation occurs in the called routine :
440 : ! iband->band, nband_k->nband, npw_k->npw, npw1_k->npw1
441 4531327 : eig0nk=eig0_k(iband)
442 4531327 : usedcwavef=gs_hamkq%usepaw;if (dim_dcwf==0) usedcwavef=0
443 4531327 : if (inonsc==1) usedcwavef=2*usedcwavef
444 4531327 : opt_gvnlx1=0;if (ipert==natom+2) opt_gvnlx1=1
445 4531327 : if (ipert==natom+2.and.gs_hamkq%usepaw==1.and.inonsc==1) opt_gvnlx1=2
446 :
447 4531327 : if ( (ipert/=natom+10 .and. ipert/=natom+11) .or. abs(occ_k(iband))>tol8 ) then
448 4524927 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
449 :
450 39879865 : bands_treated_now = 0
451 4524927 : bands_treated_now(iband) = 1
452 4524927 : call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr)
453 :
454 4524927 : if (dtset%rf2_dkdk==2 .and. (idir==1 .or. idir==2 .or. idir==3)) then
455 1584 : eig1_k = zero
456 48 : resid = zero
457 : else
458 : call dfpt_cgwf(iband,iband_me,rank_band,bands_treated_now,dtset%berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,&
459 : & rf2,dcwavef,&
460 : & eig0_k,eig0_kq,eig1_k,gh0c1,gh1c_n,grad_berry,gsc,gscq,gs_hamkq,gvnlxc,gvnlx1,icgq,&
461 : & idir,ipert,igscq,mcgq,mgscq,mpi_enreg,mpw1,natom,nband_k,nband_me,dtset%nbdbuf,dtset%nline,&
462 : & npw_k,npw1_k,nspinor,opt_gvnlx1,prtvol,quit,resid,rf_hamkq,dtset%dfpt_sciss,dtset%tolrde,&
463 4524879 : & dtset%tolwfr,usedcwavef,dtset%wfoptalg,nlines_done,eta=eta_,omega=omega)
464 : end if
465 4524927 : resid_k(iband)=resid
466 :
467 : else
468 6400 : resid_k(iband)=zero
469 : end if
470 :
471 4531327 : if (ipert/=natom+10 .and. ipert/= natom+11) then
472 : ! At this stage, the 1st order function cwavef is orthogonal to cgq (unlike
473 : ! when it is input to dfpt_cgwf). Here, restore the "active space" content
474 : ! of the first-order wavefunction, to give cwave1.
475 : ! PAW: note that dcwavef (1st-order change of WF due to overlap change)
476 : ! remains in the subspace orthogonal to cgq
477 4465167 : call proc_distrb_cycle_bands(cycle_bands, mpi_enreg%proc_distrb,ikpt,isppol,me)
478 4465167 : if (dtset%prtfull1wf > 0) then
479 : call full_active_wf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,eig1_k,fermie1,&
480 : eig0nk,eig0_kq,dtset%elph2_imagden,iband,ibgq,icgq,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,&
481 573 : 0,gs_hamkq%usepaw)
482 14325 : edocc_k=zero
483 573 : tocceig=1
484 : else
485 : call corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,edocc_k,eig1_k,fermie1,gh0c1,&
486 : iband,ibgq,icgq,gs_hamkq%istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,&
487 4464594 : occ_k,rocceig,0,gs_hamkq%usepaw,tocceig)
488 : end if
489 4465167 : ABI_FREE (cycle_bands)
490 : else
491 66160 : tocceig=0
492 66160 : call cg_zcopy(npw1_k*nspinor,cwavef,cwave1)
493 66160 : if (gs_hamkq%usepaw==1) then
494 8960 : call pawcprj_copy(cwaveprj,cwaveprj1)
495 : end if
496 : end if
497 :
498 4531327 : if (abs(occ_k(iband))<= tol8) then
499 1165082 : ek0_k(iband)=zero
500 1165082 : ek1_k(iband)=zero
501 1165082 : eeig0_k(iband)=zero
502 1165082 : end0_k(iband)=zero
503 1165082 : end1_k(iband)=zero
504 1165082 : enl0_k(iband)=zero
505 1165082 : enl1_k(iband)=zero
506 1165082 : eloc0_k(iband)=zero
507 1165082 : evxctau0_k(iband)=zero
508 1165082 : evxctau1_k(iband)=zero
509 1165082 : nskip=nskip+1
510 : else
511 : ! Compute the 0-order kinetic operator contribution (with cwavef)
512 : call meanvalue_g(ar,kinpw1,0,gs_hamkq%istwf_k,mpi_enreg,npw1_k,nspinor,cwavef,cwavef,0,&
513 3366245 : & gpu_thread_limit=dtset%gpu_thread_limit)
514 : ! There is an additional factor of 2 with respect to the bare matrix element
515 3366245 : ek0_k(iband)=energy_factor*ar
516 : ! Compute the 1-order kinetic operator contribution (with cwave1 and cwave0), if needed.
517 : ! Note that this is called only for ddk or strain, so that npw1_k=npw_k
518 3366245 : if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4)then
519 : call matrixelmt_g(ai,ar,rf_hamkq%dkinpw_k,gs_hamkq%istwf_k,0,npw_k,nspinor,cwave1,cwave0,&
520 961896 : & mpi_enreg%me_g0, mpi_enreg%comm_fft)
521 : ! There is an additional factor of 4 with respect to the bare matrix element
522 961896 : ek1_k(iband)=two*energy_factor*ar
523 : end if
524 :
525 : ! Compute the 0-order nuclear dipole contribution (with cwavef)
526 : ! only relevant for DDK
527 3366245 : if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(gs_hamkq%vectornd)) ) then
528 222804 : ABI_MALLOC(ghc_vectornd,(2,npw_k*nspinor))
529 : ! ndat hard-coded as 1
530 : call getghc_nucdip(cwavef,ghc_vectornd,gs_hamkq%gbound_k,gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kpt_k,&
531 : & gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw_k,gs_hamkq%nvloc,&
532 74268 : & gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,nspinor,gs_hamkq%vectornd,gs_hamkq%vlocal,gs_hamkq%zora,gs_hamkq%gpu_option)
533 : ! There is an additional factor of 2 with respect to the bare matrix element
534 : end0_k(iband)=energy_factor*(DOT_PRODUCT(cwavef(1,1:npw_k*nspinor),ghc_vectornd(1,1:npw_k*nspinor))+&
535 33233892 : & DOT_PRODUCT(cwavef(2,1:npw_k*nspinor),ghc_vectornd(2,1:npw_k*nspinor)))
536 74268 : ABI_FREE(ghc_vectornd)
537 : end if
538 :
539 : ! Compute the 0-order vxctau contribution (with cwavef)
540 : ! only relevant for DDK
541 3366245 : if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(gs_hamkq%vxctaulocal)) ) then
542 24192 : ABI_MALLOC(ghc_vxctau,(2,npw_k*nspinor))
543 : ! ndat hard-coded as 1
544 : call getghc_mGGA(cwavef,ghc_vxctau,gs_hamkq%gbound_k,gs_hamkq%gprimd,gs_hamkq%istwf_k,&
545 : & gs_hamkq%kg_k,gs_hamkq%kpt_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw_k,&
546 : & gs_hamkq%nvloc,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,nspinor,gs_hamkq%vxctaulocal,&
547 8064 : & gs_hamkq%gpu_option)
548 : ! There is an additional factor of 2 with respect to the bare matrix element
549 : evxctau0_k(iband)=energy_factor*(DOT_PRODUCT(cwavef(1,1:npw_k*nspinor),ghc_vxctau(1,1:npw_k*nspinor))+&
550 3939264 : & DOT_PRODUCT(cwavef(2,1:npw_k*nspinor),ghc_vxctau(2,1:npw_k*nspinor)))
551 8064 : ABI_FREE(ghc_vxctau)
552 : end if
553 :
554 : ! Compute the 1-order nuclear dipole contribution (with cwave1 and cwave0), if needed.
555 : ! only relevant for DDK
556 3366245 : if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(rf_hamkq%vectornd)) ) then
557 222804 : ABI_MALLOC(ghc_vectornd,(2,npw1_k*nspinor))
558 : ! ndat hard-coded as 1
559 : call getgh1ndc(cwave1,ghc_vectornd,gs_hamkq%gbound_k,gs_hamkq%istwf_k,gs_hamkq%kg_k,&
560 : & gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw1_k,gs_hamkq%nvloc,&
561 : & gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,nspinor,rf_hamkq%vectornd,&
562 74268 : & gs_hamkq%vlocal,rf_hamkq%zora,gs_hamkq%gpu_option)
563 : ! There is an additional factor of 4 with respect to the bare matrix element
564 : end1_k(iband)=two*energy_factor*(DOT_PRODUCT(cwave0(1,1:npw_k*nspinor),ghc_vectornd(1,1:npw_k*nspinor))+&
565 33233892 : & DOT_PRODUCT(cwave0(2,1:npw_k*nspinor),ghc_vectornd(2,1:npw_k*nspinor)))
566 74268 : ABI_FREE(ghc_vectornd)
567 : end if
568 :
569 : ! Compute the 1-order vxctau contribution (with cwave1 and cwave0), if needed.
570 : ! only relevant for DDK
571 3366245 : if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(rf_hamkq%vxctaulocal)) ) then
572 24192 : ABI_MALLOC(ghc_vxctau,(2,npw1_k*nspinor))
573 : ! ndat hard-coded as 1
574 : call getgh1c_mGGA(cwave1,gs_hamkq%gbound_k,ghc_vxctau,gs_hamkq%gmet,gs_hamkq%gprimd,idir,gs_hamkq%istwf_k,&
575 : & gs_hamkq%kg_k,gs_hamkq%kpt_k,gs_hamkq%mgfft,mpi_enreg,nspinor,gs_hamkq%n4,gs_hamkq%n5,&
576 8064 : & gs_hamkq%n6,1,gs_hamkq%ngfft,npw_k,gs_hamkq%nvloc,rf_hamkq%vxctaulocal,gs_hamkq%gpu_option)
577 : ! There is an additional factor of 4 with respect to the bare matrix element
578 : evxctau1_k(iband)=two*energy_factor*(DOT_PRODUCT(cwave0(1,1:npw_k*nspinor),ghc_vxctau(1,1:npw_k*nspinor))+&
579 3939264 : & DOT_PRODUCT(cwave0(2,1:npw_k*nspinor),ghc_vxctau(2,1:npw_k*nspinor)))
580 8064 : ABI_FREE(ghc_vxctau)
581 : end if
582 : !
583 : ! Compute eigenvalue part of total energy (with cwavef)
584 3366245 : if (gs_hamkq%usepaw==1) then
585 : call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gsc,mpi_enreg%me_g0,&
586 501319 : & mpi_enreg%comm_spinorfft)
587 : else
588 : call sqnorm_g(scprod,gs_hamkq%istwf_k,npw1_k*nspinor,cwavef,mpi_enreg%me_g0,&
589 2864926 : & mpi_enreg%comm_fft)
590 : end if
591 3366245 : eeig0_k(iband)=-energy_factor*(eig0_k(iband)- (dtset%dfpt_sciss) +omega_)*scprod
592 :
593 : ! Compute nonlocal psp contributions to nonlocal energy:
594 : ! <G|Vnl+VFockACE|C1nk(perp)> is contained in gvnlxc (with cwavef)
595 : call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gvnlxc,mpi_enreg%me_g0,&
596 3366245 : & mpi_enreg%comm_spinorfft)
597 3366245 : enl0_k(iband)=energy_factor*scprod
598 :
599 3366245 : if(ipert/=natom+10.and.ipert/=natom+11) then
600 : ! <G|Vnl1|Cnk> is contained in gvnlx1 (with cwave1)
601 : ! gvnlx1 contains at this stage first order kinetic energy, first order nuclear dipole,
602 : ! first order vxctau1
603 : call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,gvnlx1,mpi_enreg%me_g0,&
604 3306485 : & mpi_enreg%comm_spinorfft)
605 3306485 : enl1_k(iband)=two*energy_factor*scprod
606 : end if
607 :
608 : ! Removal of the 1st-order kinetic energy, dipole, and vxctau1 from the 1st-order non-local part.
609 : ! note that in getgh1c, the "nonlocal" piece, gvnlx1, contains also first order
610 : ! kinetic, nuclear dipole, and vxctau1
611 3366245 : if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4) then
612 961896 : enl1_k(iband)=enl1_k(iband)-ek1_k(iband)-end1_k(iband)-evxctau1_k(iband)
613 : end if
614 :
615 : ! Accumulate 1st-order density (only at the last inonsc)
616 : ! Accumulate zero-order potential part of the 2nd-order total energy
617 : ! BUGFIX from Max Stengel: need to initialize eloc at each inonsc iteration, in case nnonsc > 1
618 3366245 : eloc0_k(iband) = zero
619 3366245 : option=2;if (iscf_mod>0.and.inonsc==nnsclo_now) option=3
620 : call dfpt_accrho(cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,eloc0_k(iband),&
621 : & gs_hamkq,iband,idir,ipert,isppol,dtset%kptopt,mpi_enreg,1,natom,nband_k,ncpgr,&
622 3366245 : & npw_k,npw1_k,nspinor,occ_k,option,pawrhoij1,rhoaug1,tim_fourwf,tocceig,wtk_k)
623 3366245 : if(ipert==natom+10.or.ipert==natom+11) eloc0_k(iband)=energy_factor*eloc0_k(iband)/two
624 :
625 3366245 : if(ipert==natom+10.or.ipert==natom+11) then
626 59760 : shift_band=(iband-1)*npw1_k*nspinor
627 : call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,&
628 59760 : & rf2%RHS_Stern(:,1+shift_band:npw1_k*nspinor+shift_band),mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
629 59760 : ek1_k(iband)=two*energy_factor*scprod
630 : end if
631 :
632 : end if ! End of non-zero occupation
633 :
634 : ! Exit loop over inonsc if converged and non-self-consistent
635 8407022 : if (iscf_mod<0 .and. resid<dtset%tolwfr) exit
636 :
637 : end do ! End loop over inonsc
638 :
639 : ! Get first-order eigenvalues and wavefunctions
640 4429964 : ptr = 1+(iband_me-1)*npw1_k*nspinor+icg1
641 4429964 : if (.not. present(cg1_out)) then
642 4428812 : call cg_zcopy(npw1_k*nspinor,cwave1,cg1(1,ptr))
643 : end if
644 4429964 : if(dim_eig2rf > 0) then
645 319143 : if (.not. present(cg1_out)) then
646 407758863 : cg1_active(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)=cwavef(:,:)
647 : end if
648 408145791 : gh1c_set(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)=gh1c_n(:,:)
649 408145791 : gh0c1_set(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)=gh0c1(:,:)
650 : end if
651 :
652 : ! PAW: write first-order projected wavefunctions
653 5274801 : if (psps%usepaw==1.and.gs_hamkq%usecprj==1) then
654 : call pawcprj_put(gs_hamkq%atindx,cwaveprj,cprj1,natom,iband_me,ibg1,ikpt,iorder_cprj1,isppol,&
655 598564 : & mband_mem,mk1mem,natom,1,nband_me,gs_hamkq%dimcprj,nspinor,nsppol,dtfil%unpaw1)
656 : !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,to_be_gathered=.true.)
657 : end if
658 :
659 : end do
660 :
661 : !======================================================================
662 : !================== END LOOP OVER BANDS ==============================
663 : !======================================================================
664 :
665 : ! select eig1 matrix for only my slice of bands at present k-point
666 : ! this enables global xmpi_sum in dfpt_vtorho without double counting
667 844837 : ii = 0
668 5280849 : do iband=1, nband_k
669 4436012 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) then
670 296352 : eig1_k(ii+1:ii+2*nband_k) = zero
671 : end if
672 5280849 : ii = ii + 2*nband_k
673 : end do
674 :
675 : ! call xmpi_sum(eig1_k,mpi_enreg%comm_band,ierr)
676 :
677 : ! NB: no need to sum eXX_k over band communicator here, as it is a sub-comm of kpt,
678 : ! and full mpi_sum will be done higher up.
679 :
680 : !For rf2 perturbation
681 844837 : if(ipert==natom+10.or.ipert==natom+11) call rf2_destroy(rf2)
682 :
683 : !Find largest resid over bands at this k point
684 6125686 : residk=maxval(resid_k(:))
685 844837 : if (prtvol>2 .or. ikpt<=nkpt_max) then
686 108110 : do ii=0,(nband_k-1)/8
687 56928 : write(msg,'(1p,8e10.2)')(resid_k(iband),iband=1+ii*8,min(nband_k,8+ii*8))
688 901765 : call wrtout(std_out,msg)
689 : end do
690 : end if
691 :
692 844837 : call timab(139,2,tsec)
693 844837 : call timab(130,1,tsec)
694 :
695 844837 : ABI_FREE(cwave0)
696 844837 : ABI_FREE(cwavef)
697 844837 : ABI_FREE(cwave1)
698 844837 : ABI_FREE(gh0c1)
699 844837 : ABI_FREE(gvnlxc)
700 844837 : ABI_FREE(gvnlx1)
701 844837 : ABI_FREE(gh1c_n)
702 :
703 844837 : if (gs_hamkq%usepaw==1) then
704 112765 : call pawcprj_free(cwaveprj)
705 112765 : call pawcprj_free(cwaveprj1)
706 112765 : if (gs_hamkq%usecprj==1) then
707 112765 : call pawcprj_free(cwaveprj0)
708 : end if
709 : end if
710 844837 : ABI_FREE(dcwavef)
711 : #ifdef HAVE_OPENMP_OFFLOAD
712 : !$OMP TARGET EXIT DATA MAP(delete:gscq) IF(gs_hamkq%usepaw==1 .and. gs_hamkq%gpu_option==ABI_GPU_OPENMP)
713 : #endif
714 844837 : ABI_FREE(gscq)
715 844837 : ABI_FREE(gsc)
716 1081908 : ABI_FREE(cwaveprj0)
717 1081908 : ABI_FREE(cwaveprj)
718 1081908 : ABI_FREE(cwaveprj1)
719 :
720 :
721 : !###################################################################
722 :
723 : ! Write the number of one-way 3D ffts skipped until now (in case of fixed occupation numbers)
724 844837 : call xmpi_sum(nskip,mpi_enreg%comm_band,ierr)
725 844837 : if (iscf_mod>0 .and. (prtvol>2 .or. ikpt<=nkpt_max)) then
726 49882 : write(msg,'(a,i0)')' dfpt_vtowfk : number of one-way 3D ffts skipped in vtowfk3 until now =',nskip
727 49882 : call wrtout(std_out,msg)
728 : end if
729 :
730 844837 : if (prtvol<=2 .and. ikpt==nkpt_max+1) then
731 0 : write(msg,'(3a)') ch10,' dfpt_vtowfk : prtvol=0, 1 or 2, do not print more k-points.',ch10
732 0 : call wrtout(std_out,msg)
733 : end if
734 :
735 844837 : if (residk>dtset%tolwfr .and. iscf_mod<=0 .and. iscf_mod/=-3) then
736 5016 : write(msg,'(a,2i0,a,es13.5)')'Wavefunctions not converged for nnsclo,ikpt=',nnsclo_now,ikpt,' max resid=',residk
737 5016 : ABI_WARNING(msg)
738 : end if
739 :
740 844837 : call timab(130,2,tsec)
741 844837 : call timab(128,2,tsec)
742 :
743 : ABI_NVTX_END_RANGE()
744 :
745 : DBG_EXIT('COLL')
746 :
747 2534511 : end subroutine dfpt_vtowfk
748 : !!***
749 :
750 : !!****f* ABINIT/corrmetalwf1
751 : !!
752 : !! NAME
753 : !! corrmetalwf1
754 : !!
755 : !! FUNCTION
756 : !! Response function calculation only:
757 : !! Correct 1st-order wave-function, taking into account "metallic" occupations.
758 : !! 1st-order WF orthogonal to C_n,k+q, restore the "active space" content of the first-order WF.
759 : !! receives a single band at k as input, and works on all bands at k+q
760 : !!
761 : !! INPUTS
762 : !! cgq(2,mcgq)=planewave coefficients of wavefunctions at k+q
763 : !! cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
764 : !! cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
765 : !! cwaveprj(natom,nspinor)= 1st-order wave-function before correction
766 : !! projected on NL projectors (PAW)
767 : !! cycle_bands(nband)=array of logicals for bands we have on this cpu
768 : !! eig1(2*nband**2)=first-order eigenvalues (hartree)
769 : !! fermie1=derivative of fermi energy wrt (strain) perturbation
770 : !! ghc(2,npw1*nspinor)=<G|H0-eig0_k.I|C1 band,k> (NCPP) or <G|H0-eig0_k.S0|C1 band,k> (PAW)
771 : !! (C1 before correction)
772 : !! iband=index of current band
773 : !! ibgq=shift to be applied on the location of data in the array cprjq
774 : !! icgq=shift to be applied on the location of data in the array cgq
775 : !! istwf_k=option parameter that describes the storage of wfs
776 : !! mcgq=second dimension of the cgq array
777 : !! mcprjq=second dimension of the cprjq array
778 : !! mpi_enreg=information about MPI parallelization
779 : !! natom=number of atoms in cell
780 : !! nband=number of bands
781 : !! npw1=number of plane waves at this k+q point
782 : !! nspinor=number of spinorial components of the wavefunctions
783 : !! occ(nband)=occupation number for each band for each k.
784 : !! rocceig(nband,nband)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
785 : !! if this ratio has been attributed to the band n (second argument), zero otherwise
786 : !! timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
787 : !! usepaw=flag for PAW
788 : !!
789 : !! OUTPUT
790 : !! cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
791 : !! cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
792 : !! projected on NL projectors (PAW)
793 : !! edocc(nband)=correction to 2nd-order total energy coming from changes of occupations
794 : !! wf_corrected=flag put to 1 if input cwave1 is effectively different from output cwavef
795 : !!
796 : !! NOTES
797 : !! Was part of dfpt_vtowfk before.
798 : !!
799 : !! SOURCE
800 :
801 4464594 : subroutine corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,edocc,eig1,fermie1,ghc,iband, &
802 4464594 : & ibgq,icgq,istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,nspinor,occ,rocceig,timcount,&
803 : & usepaw,wf_corrected)
804 :
805 : !Arguments ------------------------------------
806 : !scalars
807 : integer,intent(in) :: iband,ibgq,icgq,istwf_k,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw
808 : integer,intent(out) :: wf_corrected
809 : real(dp),intent(in) :: fermie1
810 : type(MPI_type),intent(in) :: mpi_enreg
811 : !arrays
812 : logical,intent(in) :: cycle_bands(nband)
813 : real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor)
814 : real(dp),intent(in) :: eig1(2*nband**2),ghc(2,npw1*nspinor),occ(nband),rocceig(nband,nband)
815 : real(dp),intent(out) :: cwave1(2,npw1*nspinor),edocc(nband)
816 : type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw)
817 : type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw) !vz_i
818 :
819 : !Local variables-------------------------------
820 : !scalars
821 : integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii
822 : integer :: ibandkq_me, ierr, iband_
823 : integer :: wf_corrected_
824 : real(dp) :: facti,factr,invocc
825 : real(dp) :: edocc_tmp
826 : !arrays
827 8929188 : integer :: bands_treated_now(nband)
828 4464594 : integer, allocatable :: nlmn(:)
829 : real(dp) :: tsec(2)
830 4464594 : real(dp),allocatable :: cwcorr(:,:)
831 15976316 : type(pawcprj_type) :: cwaveprj1_corr(natom,nspinor*usepaw)
832 : ! *********************************************************************
833 :
834 : DBG_ENTER("COLL")
835 :
836 4464594 : call timab(214+timcount,1,tsec)
837 :
838 39512340 : bands_treated_now = 0
839 4464594 : bands_treated_now(iband) = 1
840 4464594 : call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr)
841 :
842 : !At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf).
843 : !Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 .
844 :
845 13393782 : ABI_MALLOC(cwcorr,(2,npw1*nspinor))
846 :
847 4464594 : wf_corrected=0
848 :
849 4464594 : if(usepaw==1) then
850 1768812 : ABI_MALLOC(nlmn,(natom))
851 1822052 : nlmn(:)=cwaveprj1(:,1)%nlmn
852 589604 : call pawcprj_alloc(cwaveprj1_corr, cwaveprj1(1,1)%ncpgr, nlmn)
853 589604 : ABI_FREE(nlmn)
854 : end if
855 :
856 : ! loop iband_ over all bands being treated for the moment
857 : ! all procs in pool of bands should be working on the same iband_ at a given time
858 : ! I will save in _my_ array cwave1, if iband==iband_
859 39512340 : do iband_ = 1, nband
860 35047746 : if (bands_treated_now(iband_) == 0) cycle
861 :
862 2796839541 : cwcorr = zero
863 4470642 : if (usepaw==1) then
864 592772 : call pawcprj_set_zero(cwaveprj1_corr)
865 : end if
866 4470642 : edocc_tmp = zero
867 4470642 : wf_corrected_=0
868 :
869 : !Correct WF only for occupied states
870 4470642 : if (abs(occ(iband_)) > tol8) then
871 3312437 : invocc=one/occ(iband_)
872 :
873 :
874 : ! Loop over WF at k+q subspace
875 3312437 : ibandkq_me = 0
876 24928918 : do ibandkq=1,nband
877 21616481 : if(cycle_bands(ibandkq)) cycle
878 21471329 : ibandkq_me = ibandkq_me + 1
879 :
880 : ! Select bands with variable occupation
881 24783766 : if (abs(rocceig(ibandkq,iband_))>tol8) then
882 :
883 4226013 : wf_corrected_=1
884 :
885 4226013 : index_eig1=2*ibandkq-1+(iband_-1)*2*nband
886 4226013 : index_cgq=npw1*nspinor*(ibandkq_me-1)+icgq
887 :
888 4226013 : if(ibandkq==iband_) then
889 160043 : factr=rocceig(ibandkq,iband_)*invocc*(eig1(index_eig1)-fermie1)
890 : else
891 4065970 : factr=rocceig(ibandkq,iband_)*invocc*eig1(index_eig1)
892 : end if
893 4226013 : facti = rocceig(ibandkq,iband_)*invocc*eig1(index_eig1+1)
894 :
895 : ! Apply correction to 1st-order WF
896 : !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwcorr,facti,factr,index_cgq,npw1,nspinor)
897 1609527787 : do ii=1,npw1*nspinor
898 1605301774 : cwcorr(1,ii)=cwcorr(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq))
899 1609527787 : cwcorr(2,ii)=cwcorr(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq))
900 : end do
901 :
902 : ! In the PAW case, also apply correction to projected WF
903 4226013 : if (usepaw==1) then
904 517773 : index_cprjq=nspinor*(ibandkq_me-1)+ibgq
905 1553319 : call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1_corr)
906 : end if
907 :
908 : ! The factor of two is needed because we compute the 2DTE, and not E(2)
909 4226013 : edocc_tmp = edocc_tmp-two*(factr*eig1(index_eig1)+facti*eig1(index_eig1+1))
910 :
911 : end if ! Variable occupations
912 : end do ! Loop over k+q subspace
913 : end if ! if occupied states
914 :
915 : ! 3) reduce over bands to get all contributions to correction
916 : ! need MPI reduce over band communicator only
917 4470642 : call xmpi_sum(cwcorr, mpi_enreg%comm_band, ierr)
918 4470642 : if (usepaw==1) then
919 592772 : call pawcprj_mpi_sum(cwaveprj1_corr, mpi_enreg%comm_band, ierr)
920 : end if
921 :
922 : ! this sums over the k+q contributions to the present iband_
923 4470642 : call xmpi_sum(edocc_tmp, mpi_enreg%comm_band, ierr)
924 4470642 : call xmpi_sum(wf_corrected_, mpi_enreg%comm_band, ierr)
925 :
926 :
927 : ! 4) add correction to the cwave1
928 : ! if I have iband_, correct my cwave1
929 13405878 : if (iband_==iband) then
930 4464594 : if (wf_corrected_ > 0) wf_corrected = 1
931 4464594 : edocc(iband) = edocc_tmp
932 2743112565 : cwave1 = cwcorr
933 4464594 : call cg_zaxpy(npw1*nspinor,(/one,zero/),cwavef,cwave1)
934 : !Idem for cprj
935 4464594 : if (usepaw==1) then
936 589604 : call pawcprj_copy(cwaveprj,cwaveprj1)
937 589604 : call pawcprj_zaxpby((/one,zero/),(/one,zero/),cwaveprj1_corr,cwaveprj1)
938 : end if
939 : end if
940 : end do ! loop over all bands presently running in parallel
941 :
942 :
943 : !In the PAW case, compute <Psi^(1)_ortho|H-Eig0_k.S|Psi^(1)_parallel> contribution to 2DTE
944 4464594 : if (usepaw==1.and.wf_corrected==1) then
945 : !$OMP WORKSHARE
946 79365317 : cwcorr(:,:)=cwave1(:,:)-cwavef(:,:)
947 : !$OMP END WORKSHARE
948 114155 : call dotprod_g(factr,facti,istwf_k,npw1*nspinor,1,cwcorr,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
949 114155 : edocc(iband)=edocc(iband)+four*factr
950 : end if
951 :
952 :
953 4464594 : ABI_FREE(cwcorr)
954 4464594 : if (usepaw==1) then
955 589604 : call pawcprj_free(cwaveprj1_corr)
956 : end if
957 :
958 4464594 : call timab(214+timcount,2,tsec)
959 :
960 : DBG_EXIT("COLL")
961 :
962 4464594 : end subroutine corrmetalwf1
963 : !!***
964 :
965 : end module m_dfpt_vtowfk
966 : !!***
|
