Branch data 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 : 1690778 : subroutine dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,&
178 : : & dim_eig2rf,dtfil,dtset,&
179 : 1690778 : & edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,&
180 : 845389 : & ek0_k,ek1_k,eloc0_k,end0_k,end1_k,enl0_k,enl1_k,evxctau0_k,evxctau1_k,&
181 [ - + + + ]: 845389 : & 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 : 845389 : & 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 : 845389 : type(rf2_t) :: rf2
244 : : !arrays
245 : 845389 : logical,allocatable :: cycle_bands(:)
246 : 2536167 : integer :: rank_band(nband_k), bands_treated_now(nband_k)
247 : : real(dp) :: tsec(2)
248 : 845389 : real(dp),allocatable :: cwave0(:,:),cwave1(:,:),cwavef(:,:)
249 : 845389 : real(dp),allocatable :: dcwavef(:,:),gh1c_n(:,:),gh0c1(:,:),ghc_vectornd(:,:)
250 : 845389 : real(dp),allocatable :: ghc_vxctau(:,:)
251 : 845389 : real(dp),allocatable :: gsc(:,:),gscq(:,:),gvnlx1(:,:),gvnlxc(:,:)
252 : 845389 : real(dp),pointer :: kinpw1(:)
253 : 845389 : 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 : 845389 : call timab(128,1,tsec)
262 : :
263 : 845389 : nkpt_max=50; if (xmpi_paral==1) nkpt_max=-1
264 : :
265 [ + + - + ]: 845389 : 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 : 845389 : me=mpi_enreg%me_kpt
272 : 845389 : quit=0
273 : :
274 : : !The value of iscf must be modified if ddk perturbation
275 [ + + + + : 845389 : iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3
+ + ]
276 : :
277 : 845389 : kinpw1 => gs_hamkq%kinpw_kp
278 [ - + + - : 2536167 : ABI_MALLOC(gh0c1,(2,npw1_k*nspinor))
+ - + - -
+ - + ]
279 [ - + + - : 1690778 : ABI_MALLOC(gvnlxc,(2,npw1_k*nspinor))
- + - + ]
280 [ - + + - : 1690778 : ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor))
- + - + ]
281 [ - + + - : 2536167 : ABI_MALLOC(cwave0,(2,npw_k*nspinor))
+ - + - -
+ - + ]
282 [ - + + - : 1690778 : ABI_MALLOC(cwavef,(2,npw1_k*nspinor))
- + - + ]
283 [ - + + - : 1690778 : ABI_MALLOC(cwave1,(2,npw1_k*nspinor))
- + - + ]
284 [ - + + - : 1690778 : ABI_MALLOC(gh1c_n,(2,npw1_k*nspinor))
- + - + ]
285 [ + + ]: 845389 : if (gs_hamkq%usepaw==1) then
286 [ - + + - : 226634 : ABI_MALLOC(gsc,(2,npw1_k*nspinor))
- + - + ]
287 : : else
288 [ - + ]: 732072 : ABI_MALLOC(gsc,(0,0))
289 : : end if
290 : :
291 [ + - ]: 845389 : omega_=zero; if(present(omega))omega_=omega
292 [ + - ]: 845389 : eta_=zero; if(present(eta))eta_=eta
293 : :
294 : : !Read the npw and kg records of wf files
295 : 845389 : 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 [ + + + + : 3381556 : & ipert==natom+10.or.ipert==natom+11) then
+ - + + +
+ + + ]
300 : 104202 : test_ddk=1
301 [ + + + + ]: 845389 : if(ipert==natom+10.or.ipert==natom+11) test_ddk=0
302 : : end if
303 : :
304 : : !Additional stuff for PAW
305 [ - + ]: 845389 : ABI_MALLOC(cwaveprj0,(0,0))
306 [ + + ]: 845389 : if (gs_hamkq%usepaw==1) then
307 : : ! 1-Compute all <g|S|Cnk+q>
308 : 113317 : igscq=0
309 : : !TODO MJV: PAW mband_mem
310 : 113317 : mgscq=mpw1*nspinor*mband_mem
311 [ - + + - : 339951 : 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 : 113317 : & 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 : 113317 : iorder_cprj=0;iorder_cprj1=0
322 [ + + + + : 113317 : dim_dcwf=npw1_k*nspinor;if (ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) dim_dcwf=0
+ + ]
323 [ + + + - : 339951 : ABI_MALLOC(dcwavef,(2,dim_dcwf))
+ - + + -
+ - + ]
324 [ + - ]: 113317 : if (gs_hamkq%usecprj==1) then
325 [ - + ]: 113317 : ABI_FREE(cwaveprj0)
326 [ - + + - : 808762 : ABI_MALLOC(cwaveprj0,(natom,nspinor))
- + + - +
- + - - +
- + + + +
+ ]
327 : 113317 : call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj)
328 : : end if
329 [ - + + - : 808762 : ABI_MALLOC(cwaveprj,(natom,nspinor))
- + + - +
- + - - +
- + + + +
+ ]
330 [ - + + - : 695445 : ABI_MALLOC(cwaveprj1,(natom,nspinor))
- + + - -
+ - + + +
+ + ]
331 : 113317 : call pawcprj_alloc(cwaveprj ,0,gs_hamkq%dimcprj)
332 : 113317 : 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 : 845389 : energy_factor=two
342 [ + + + + ]: 845389 : if(ipert==natom+10.or.ipert==natom+11) energy_factor=six
343 : :
344 : : !For rf2 perturbation :
345 [ + + + + ]: 845389 : 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 : 845389 : call timab(139,1,tsec)
351 : :
352 : : !======================================================================
353 : : !================== LOOP OVER BANDS ==================================
354 : : !======================================================================
355 : :
356 : 845389 : 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 : 845389 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
359 : :
360 : 845389 : iband_me = 0
361 [ + + ]: 5296089 : do iband=1,nband_k
362 : :
363 : : ! Skip bands not treated by current proc
364 [ + + ]: 4450700 : if( (mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me)) cycle
365 : 4444436 : iband_me = iband_me + 1
366 : :
367 : : !unit_me = 300+iband
368 : : !Get ground-state wavefunctions
369 : 4444436 : ptr = 1+(iband_me-1)*npw_k*nspinor+icg
370 : 4444436 : call cg_zcopy(npw_k*nspinor,cg(1,ptr),cwave0)
371 : :
372 : : ! Get PAW ground state projected WF (cprj)
373 [ + + + + : 4444436 : if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.ipert/=natom+10.and.ipert/=natom+11) then
+ + ]
374 : 604076 : idir0 = idir
375 [ + + + + ]: 604076 : 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 : 604076 : & icpgr=idir0,ncpgr=ncpgr)
381 : : end if
382 : :
383 : : ! Get first-order wavefunctions
384 : 4444436 : ptr = 1+(iband_me-1)*npw1_k*nspinor+icg1
385 : 4444436 : 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 [ + + ]: 9032310 : do ispinor=1,nspinor
399 : 4587874 : igs=(ispinor-1)*npw1_k
400 [ + + ]: 857902296 : do ipw=1+igs,npw1_k+igs
401 [ - + ]: 853457860 : 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 [ + + ]: 4444436 : 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 [ + + ]: 8435951 : 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 : 4545799 : eig0nk=eig0_k(iband)
442 [ + + ]: 4545799 : usedcwavef=gs_hamkq%usepaw;if (dim_dcwf==0) usedcwavef=0
443 [ + + ]: 4545799 : if (inonsc==1) usedcwavef=2*usedcwavef
444 [ + + ]: 4545799 : opt_gvnlx1=0;if (ipert==natom+2) opt_gvnlx1=1
445 [ + + + + : 4545799 : if (ipert==natom+2.and.gs_hamkq%usepaw==1.and.inonsc==1) opt_gvnlx1=2
+ - ]
446 : :
447 [ + + + + : 4545799 : if ( (ipert/=natom+10 .and. ipert/=natom+11) .or. abs(occ_k(iband))>tol8 ) then
+ + ]
448 : 4539399 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
449 : :
450 [ + + ]: 40326625 : bands_treated_now = 0
451 : 4539399 : bands_treated_now(iband) = 1
452 : 4539399 : call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr)
453 : :
454 [ + + + + ]: 4539399 : 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 : 4539351 : & dtset%tolwfr,usedcwavef,dtset%wfoptalg,nlines_done,eta=eta_,omega=omega)
464 : : end if
465 : 4539399 : resid_k(iband)=resid
466 : :
467 : : else
468 : 6400 : resid_k(iband)=zero
469 : : end if
470 : :
471 [ + + + + ]: 4545799 : 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 [ - + ]: 4479639 : call proc_distrb_cycle_bands(cycle_bands, mpi_enreg%proc_distrb,ikpt,isppol,me)
478 [ + + ]: 4479639 : 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 : 4479066 : occ_k,rocceig,0,gs_hamkq%usepaw,tocceig)
488 : : end if
489 [ - + ]: 4479639 : 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 [ + + ]: 4545799 : if (abs(occ_k(iband))<= tol8) then
499 : 1167962 : ek0_k(iband)=zero
500 : 1167962 : ek1_k(iband)=zero
501 : 1167962 : eeig0_k(iband)=zero
502 : 1167962 : end0_k(iband)=zero
503 : 1167962 : end1_k(iband)=zero
504 : 1167962 : enl0_k(iband)=zero
505 : 1167962 : enl1_k(iband)=zero
506 : 1167962 : eloc0_k(iband)=zero
507 : 1167962 : evxctau0_k(iband)=zero
508 : 1167962 : evxctau1_k(iband)=zero
509 : 1167962 : 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 [ - + ]: 3377837 : & gpu_thread_limit=dtset%gpu_thread_limit)
514 : : ! There is an additional factor of 2 with respect to the bare matrix element
515 : 3377837 : 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 [ + + + + : 3377837 : 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 : 961968 : & mpi_enreg%me_g0, mpi_enreg%comm_fft)
521 : : ! There is an additional factor of 4 with respect to the bare matrix element
522 : 961968 : 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 [ + + + + ]: 3377837 : if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(gs_hamkq%vectornd)) ) then
528 [ - + + - : 223020 : 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 [ + - - - : 74340 : & 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 [ + + + + ]: 36419100 : & DOT_PRODUCT(cwavef(2,1:npw_k*nspinor),ghc_vectornd(2,1:npw_k*nspinor)))
536 [ - + ]: 74340 : ABI_FREE(ghc_vectornd)
537 : : end if
538 : :
539 : : ! Compute the 0-order vxctau contribution (with cwavef)
540 : : ! only relevant for DDK
541 [ + + + + ]: 3377837 : 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 [ + + + + ]: 3377837 : if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(rf_hamkq%vectornd)) ) then
557 [ - + + - : 223020 : 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 [ + - - - : 74340 : & 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 [ + + + + ]: 36419100 : & DOT_PRODUCT(cwave0(2,1:npw_k*nspinor),ghc_vectornd(2,1:npw_k*nspinor)))
566 [ - + ]: 74340 : 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 [ + + + + ]: 3377837 : 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 [ + + ]: 3377837 : 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 : 512911 : & 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 : 3377837 : 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 : 3377837 : & mpi_enreg%comm_spinorfft)
597 : 3377837 : enl0_k(iband)=energy_factor*scprod
598 : :
599 [ + + + + ]: 3377837 : 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 : 3318077 : & mpi_enreg%comm_spinorfft)
605 : 3318077 : 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 [ + + + + : 3377837 : if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4) then
+ + ]
612 : 961968 : 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 : 3377837 : eloc0_k(iband) = zero
619 [ + + + - ]: 3377837 : 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 : 3377837 : & npw_k,npw1_k,nspinor,occ_k,option,pawrhoij1,rhoaug1,tim_fourwf,tocceig,wtk_k)
623 [ + + + + ]: 3377837 : if(ipert==natom+10.or.ipert==natom+11) eloc0_k(iband)=energy_factor*eloc0_k(iband)/two
624 : :
625 [ + + + + ]: 3377837 : 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 [ + + + + ]: 8435951 : 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 : 4444436 : ptr = 1+(iband_me-1)*npw1_k*nspinor+icg1
641 [ + + ]: 4444436 : if (.not. present(cg1_out)) then
642 : 4443284 : call cg_zcopy(npw1_k*nspinor,cwave1,cg1(1,ptr))
643 : : end if
644 [ + + ]: 4444436 : 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 [ + + - + ]: 5289825 : 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 : 613036 : & 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 : 845389 : ii = 0
668 [ + + ]: 5296089 : do iband=1, nband_k
669 [ + + ]: 4450700 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) then
670 [ + + ]: 298296 : eig1_k(ii+1:ii+2*nband_k) = zero
671 : : end if
672 : 5296089 : 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 [ + + + + ]: 845389 : if(ipert==natom+10.or.ipert==natom+11) call rf2_destroy(rf2)
682 : :
683 : : !Find largest resid over bands at this k point
684 [ + - - + : 6141478 : residk=maxval(resid_k(:))
- - + + ]
685 [ + + - + ]: 845389 : 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 : 902317 : call wrtout(std_out,msg)
689 : : end do
690 : : end if
691 : :
692 : 845389 : call timab(139,2,tsec)
693 : 845389 : call timab(130,1,tsec)
694 : :
695 [ - + ]: 845389 : ABI_FREE(cwave0)
696 [ - + ]: 845389 : ABI_FREE(cwavef)
697 [ - + ]: 845389 : ABI_FREE(cwave1)
698 [ - + ]: 845389 : ABI_FREE(gh0c1)
699 [ - + ]: 845389 : ABI_FREE(gvnlxc)
700 [ - + ]: 845389 : ABI_FREE(gvnlx1)
701 [ - + ]: 845389 : ABI_FREE(gh1c_n)
702 : :
703 [ + + ]: 845389 : if (gs_hamkq%usepaw==1) then
704 : 113317 : call pawcprj_free(cwaveprj)
705 : 113317 : call pawcprj_free(cwaveprj1)
706 [ + - ]: 113317 : if (gs_hamkq%usecprj==1) then
707 : 113317 : call pawcprj_free(cwaveprj0)
708 : : end if
709 : : end if
710 [ - + ]: 845389 : 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 [ - + ]: 845389 : ABI_FREE(gscq)
715 [ - + ]: 845389 : ABI_FREE(gsc)
716 [ + - + + : 1084452 : ABI_FREE(cwaveprj0)
- + - + -
+ ]
717 [ + - + + : 1084452 : ABI_FREE(cwaveprj)
- + - + -
+ ]
718 [ + - + + : 1084452 : 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 : 845389 : call xmpi_sum(nskip,mpi_enreg%comm_band,ierr)
725 [ + + + + : 845389 : 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 [ + + - + ]: 845389 : 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 [ + + + + ]: 845389 : 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 : 845389 : call timab(130,2,tsec)
741 : 845389 : call timab(128,2,tsec)
742 : :
743 : : ABI_NVTX_END_RANGE()
744 : :
745 : : DBG_EXIT('COLL')
746 : :
747 [ - + - + : 2536167 : 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 : 4479066 : subroutine corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,edocc,eig1,fermie1,ghc,iband, &
802 : 4479066 : & 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 : 8958132 : integer :: bands_treated_now(nband)
828 : 4479066 : integer, allocatable :: nlmn(:)
829 : : real(dp) :: tsec(2)
830 : 4479066 : real(dp),allocatable :: cwcorr(:,:)
831 [ + + + + : 16308092 : type(pawcprj_type) :: cwaveprj1_corr(natom,nspinor*usepaw)
+ + + + +
- + - ]
832 : : ! *********************************************************************
833 : :
834 : : DBG_ENTER("COLL")
835 : :
836 : 4479066 : call timab(214+timcount,1,tsec)
837 : :
838 [ + + ]: 39959100 : bands_treated_now = 0
839 : 4479066 : bands_treated_now(iband) = 1
840 : 4479066 : 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 [ - + + - : 13437198 : ABI_MALLOC(cwcorr,(2,npw1*nspinor))
+ - + - -
+ - + ]
846 : :
847 : 4479066 : wf_corrected=0
848 : :
849 [ + + ]: 4479066 : if(usepaw==1) then
850 [ - + + - : 1812228 : ABI_MALLOC(nlmn,(natom))
+ - + - -
+ - + ]
851 [ + + ]: 1894196 : nlmn(:)=cwaveprj1(:,1)%nlmn
852 : 604076 : call pawcprj_alloc(cwaveprj1_corr, cwaveprj1(1,1)%ncpgr, nlmn)
853 [ - + ]: 604076 : 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 [ + + ]: 39959100 : do iband_ = 1, nband
860 [ + + ]: 35480034 : if (bands_treated_now(iband_) == 0) cycle
861 : :
862 [ + + + + ]: 2836895445 : cwcorr = zero
863 [ + + ]: 4485330 : if (usepaw==1) then
864 : 607460 : call pawcprj_set_zero(cwaveprj1_corr)
865 : : end if
866 : 4485330 : edocc_tmp = zero
867 : 4485330 : wf_corrected_=0
868 : :
869 : : !Correct WF only for occupied states
870 [ + + ]: 4485330 : if (abs(occ(iband_)) > tol8) then
871 : 3324245 : invocc=one/occ(iband_)
872 : :
873 : :
874 : : ! Loop over WF at k+q subspace
875 : 3324245 : ibandkq_me = 0
876 [ + + ]: 25287478 : do ibandkq=1,nband
877 [ + + ]: 21963233 : if(cycle_bands(ibandkq)) cycle
878 : 21817217 : ibandkq_me = ibandkq_me + 1
879 : :
880 : : ! Select bands with variable occupation
881 [ + + ]: 25141462 : if (abs(rocceig(ibandkq,iband_))>tol8) then
882 : :
883 : 4295133 : wf_corrected_=1
884 : :
885 : 4295133 : index_eig1=2*ibandkq-1+(iband_-1)*2*nband
886 : 4295133 : index_cgq=npw1*nspinor*(ibandkq_me-1)+icgq
887 : :
888 [ + + ]: 4295133 : if(ibandkq==iband_) then
889 : 160043 : factr=rocceig(ibandkq,iband_)*invocc*(eig1(index_eig1)-fermie1)
890 : : else
891 : 4135090 : factr=rocceig(ibandkq,iband_)*invocc*eig1(index_eig1)
892 : : end if
893 : 4295133 : 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 [ + + ]: 1643085547 : do ii=1,npw1*nspinor
898 : 1638790414 : cwcorr(1,ii)=cwcorr(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq))
899 : 1643085547 : 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 [ + + ]: 4295133 : if (usepaw==1) then
904 : 586893 : index_cprjq=nspinor*(ibandkq_me-1)+ibgq
905 [ + + ]: 1760679 : 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 : 4295133 : 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 : 4485330 : call xmpi_sum(cwcorr, mpi_enreg%comm_band, ierr)
918 [ + + ]: 4485330 : if (usepaw==1) then
919 : 607460 : 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 : 4485330 : call xmpi_sum(edocc_tmp, mpi_enreg%comm_band, ierr)
924 : 4485330 : 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 [ + + ]: 13449726 : if (iband_==iband) then
930 [ + + ]: 4479066 : if (wf_corrected_ > 0) wf_corrected = 1
931 : 4479066 : edocc(iband) = edocc_tmp
932 [ + + + + ]: 2768835141 : cwave1 = cwcorr
933 : 4479066 : call cg_zaxpy(npw1*nspinor,(/one,zero/),cwavef,cwave1)
934 : : !Idem for cprj
935 [ + + ]: 4479066 : if (usepaw==1) then
936 : 604076 : call pawcprj_copy(cwaveprj,cwaveprj1)
937 : 604076 : 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 [ + + + + ]: 4479066 : if (usepaw==1.and.wf_corrected==1) then
945 : : !$OMP WORKSHARE
946 [ + + + + ]: 96121157 : cwcorr(:,:)=cwave1(:,:)-cwavef(:,:)
947 : : !$OMP END WORKSHARE
948 : 125675 : call dotprod_g(factr,facti,istwf_k,npw1*nspinor,1,cwcorr,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
949 : 125675 : edocc(iband)=edocc(iband)+four*factr
950 : : end if
951 : :
952 : :
953 [ - + ]: 4479066 : ABI_FREE(cwcorr)
954 [ + + ]: 4479066 : if (usepaw==1) then
955 : 604076 : call pawcprj_free(cwaveprj1_corr)
956 : : end if
957 : :
958 : 4479066 : call timab(214+timcount,2,tsec)
959 : :
960 : : DBG_EXIT("COLL")
961 : :
962 [ - + ]: 4479066 : end subroutine corrmetalwf1
963 : : !!***
964 : :
965 : : end module m_dfpt_vtowfk
966 : : !!***
|
