Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfpt_cgwf
2 : : !! NAME
3 : : !! m_dfpt_cgwf
4 : : !!
5 : : !! FUNCTION
6 : : !! Update one single wavefunction (cwavef), non self-consistently.
7 : : !! Uses a conjugate-gradient algorithm.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 1999-2026 ABINIT group (XG,DRH,XW,FJ,MT,LB,MG)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : ! nvtx related macro definition
24 : : #include "nvtx_macros.h"
25 : :
26 : : module m_dfpt_cgwf
27 : :
28 : : use, intrinsic :: iso_c_binding
29 : : use defs_basis
30 : : use m_abicore
31 : : use m_errors
32 : : use m_xomp
33 : : use m_xmpi
34 : : use m_cgtools
35 : : use m_rf2
36 : :
37 : : use m_fstrings, only : sjoin, ftoa, itoa, ktoa
38 : : use defs_abitypes, only : MPI_type
39 : : use m_dtset, only : dataset_type
40 : : use m_time, only : timab
41 : : use m_mpinfo, only : init_mpi_enreg, destroy_mpi_enreg, copy_mpi_enreg
42 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_set_zero, pawcprj_axpby, &
43 : : pawcprj_mpi_sum, pawcprj_zaxpby
44 : : use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type, KPRIME_H_KPRIME
45 : : use m_getghc, only : getghc
46 : : use m_getgh1c, only : getgh1c, getdc1
47 : : use m_fft, only : fft_ug
48 : :
49 : : #if defined(HAVE_GPU_MARKERS)
50 : : use m_nvtx_data
51 : : #endif
52 : :
53 : : implicit none
54 : :
55 : : private
56 : : !!***
57 : :
58 : : public :: dfpt_cgwf ! Update one single wavefunction (cwavef), non self-consistently.
59 : : public :: full_active_wf1 ! Restore the full "active space" contribution to the 1st-order wavefunctions.
60 : : !!***
61 : :
62 : : !!****t* m_dfpt_cgwf/stern_t
63 : : !! NAME
64 : : !! stern_t
65 : : !!
66 : : !! FUNCTION
67 : : !! Simplified interface to the Sternheimer solver.
68 : : !! Wrapper around the dfpt_cgwf routine.
69 : : !!
70 : : !! SOURCE
71 : :
72 : : type,public :: stern_t
73 : :
74 : : integer :: npw_k = -1
75 : : integer :: npw_kq = -1
76 : : integer :: nspinor = -1
77 : : ! Number of plane-waves at k and k+q and number of spinors.
78 : :
79 : : integer :: mcgq = -1
80 : : ! second dimension of the cgq array.
81 : : integer :: mgscq = -1
82 : : ! second dimension of gscq, with only mband_mem bands.
83 : : integer :: mcprjq = -1
84 : : ! second dimension of the cprjq array.
85 : : integer :: nband = -1
86 : : ! Number of bands.
87 : : integer :: nband_me = -1
88 : : ! Number of bands treated by this MPI proc.
89 : : integer :: nline_in = -1
90 : : ! Max number of line minimization.
91 : : integer :: nlines_done = -1
92 : : ! Number of line minimization in stern_solve.
93 : : integer :: usedcwavef = 0
94 : : ! flag controlling the use of dcwavef array (PAW only):
95 : : integer :: usepaw = -1
96 : : ! 1 if PAW is used.
97 : :
98 : : integer :: work_ngfft(18)
99 : : logical :: use_cache
100 : : logical :: has_band_para
101 : :
102 : : type(dataset_type),pointer :: dtset => null()
103 : : type(mpi_type) :: mpi_enreg
104 : : !type(u1_cache_t) :: u1c
105 : :
106 : : integer,allocatable :: bands_treated_now(:)
107 : : integer,allocatable :: rank_band(:)
108 : :
109 : : real(dp),allocatable :: fermie1_idir_ipert(:,:)
110 : : real(dp),allocatable :: eig1_k(:,:,:)
111 : : ! (2, nband, nband)
112 : : ! matrix of first-order eigenvalues (hartree)
113 : : ! The array eig1_k contains:
114 : : ! <u_(jband,k+q)^(0)|H_(k+q,k)^(1)|u_(iband,k)^(0)> (NC psps)
115 : : ! <u_(jband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(iband,k)^(0)> (PAW)
116 : :
117 : : real(dp),allocatable :: dcwavef(:, :), gh1c_n(:, :), ghc(:,:), gsc(:,:), gvnlxc(:,:), gvnlx1(:,:)
118 : : real(dp),allocatable :: cgq(:,:,:), gscq(:,:,:), work(:,:,:,:)
119 : :
120 : : type(pawcprj_type),allocatable :: cprjq(:,:)
121 : : ! (natom, mcprjq)
122 : :
123 : : type(pawcprj_type),allocatable :: cwaveprj1(:,:)
124 : : ! (natom, nspinor*usepaw)
125 : :
126 : : contains
127 : :
128 : : procedure :: init => stern_init
129 : : ! Initialize the object.
130 : :
131 : : procedure :: solve => stern_solve
132 : : ! Solves the NSCF Sternheimer equation. Simplified wrapper around dfpt_cgwf.
133 : :
134 : : procedure :: free => stern_free
135 : : ! Free dynamic memory.
136 : :
137 : : end type stern_t
138 : : !!***
139 : :
140 : : contains
141 : : !!***
142 : :
143 : : !!****f* ABINIT/dfpt_cgwf
144 : : !! NAME
145 : : !! dfpt_cgwf
146 : : !!
147 : : !! FUNCTION
148 : : !! Update one single wavefunction (cwavef), non self-consistently.
149 : : !! Uses a conjugate-gradient algorithm.
150 : : !! Try to keep close to the formulas in PRB 55, 10337 (1997) [[cite:Gonze1997]], for the
151 : : !! non-self-consistent case, except that we are computing here
152 : : !! the second-derivative of the total energy, and not E(2). There
153 : : !! is a factor of 2 between the two quantities.
154 : : !! The wavefunction that is generated is always orthogonal to cgq.
155 : : !! It is orthogonal to the active Hilbert space, and will be complemented
156 : : !! by contributions from the active space in the calling routine, if needed.
157 : : !!
158 : : !! As concerns the MPI algorithm: cgq and gscq are distributed inside comm_band and
159 : : !! each proc has nband_me non-overlapping blocks. Each proc in comm_band calls dfpt_cgwf
160 : : !! with a different u^1_{band} state (the band index is therefore LOCAL) but then
161 : : !! we need to communicate every time we call projbd to orthogonalize wrt the MPI-distributed cgq.
162 : : !! Other arrays such as rank_band and bands_treated_now are GLOBAL i.e. all procs in comm_band
163 : : !! are supposed to call dfpt_cgwf with the same values.
164 : : !!
165 : : !! INPUTS
166 : : !! u1_band_=which particular band we are converging (LOCAL)
167 : : !! A negative value is used when the routine is called in band-mode with MPI-distributed cgq
168 : : !! to indicate that this proc is not optimizing abs(u1_band). Used when calling dfpt_cgw in EPH.
169 : : !! band_me=cpu-local index in cgq array of band which we are converging.
170 : : !! berryopt=option for Berry phase
171 : : !! cgq(2,mcgq)=wavefunction coefficients for MY bands at k+Q
172 : : !! cwave0(2,npw*nspinor)=GS wavefunction at k, in reciprocal space
173 : : !! cwaveprj0(natom,nspinor*usecprj)=GS wave function at k projected with nl projectors
174 : : !! rank_band(nband)=rank of processor in band_comm which have the other bands for cgq below (GLOBAL)
175 : : !! bands_treated_now(nband) (GLOBAL)
176 : : !! eig0_k=0-order eigenvalues for the present wavefunction at k
177 : : !! eig0_kq(nband)=GS eigenvalues at k+q (hartree)
178 : : !! grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
179 : : !! gscq(2,mgscq)=<g|S|Cnk+q> coefficients for MY bands (PAW) at k+Q
180 : : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+Q
181 : : !! icgq=shift to be applied on the location of data in the array cgq
182 : : !! igscq=shift to be applied on the location of data in the array gscq
183 : : !! idir=direction of the perturbation
184 : : !! ipert=type of the perturbation
185 : : !! mcgq=second dimension of the cgq array
186 : : !! mgscq=second dimension of gscq, with only mband_mem bands
187 : : !! mpi_enreg=information about MPI parallelization
188 : : !! mpw1=maximum number of planewave for first-order wavefunctions
189 : : !! natom=number of atoms in cell.
190 : : !! nband=number of bands.
191 : : !! nband_me=number of bands on this cpu
192 : : !! nbdbuf=number of buffer bands for the minimisation
193 : : !! nline=number of line minimizations per band.
194 : : !! npw=number of planewaves in basis sphere at given k.
195 : : !! npw1=number of planewaves in basis sphere at k+Q
196 : : !! nspinor=number of spinorial components of the wavefunctions
197 : : !! omega= frequency of the perturbation
198 : : !! opt_gvnlx1=option controlling the use of gvnlx1 array:
199 : : !! 0: used as an output
200 : : !! 1: used as an input: - used only for ipert=natom+2
201 : : !! NCPP: contains the ddk 1-st order WF
202 : : !! PAW: contains frozen part of 1st-order hamiltonian
203 : : !! 2: used as input/output: - used only for PAW and ipert=natom+2
204 : : !! At input: contains the ddk 1-st order WF (times i)
205 : : !! At output: contains frozen part of 1st-order hamiltonian
206 : : !! prtvol=control print volume and debugging output
207 : : !! quit= if 1, proceeds to smooth ending of the job.
208 : : !! dfpt_sciss=scissor shift (Ha)
209 : : !! tolrde=tolerance on the ratio of differences of energies (for the line minimisation)
210 : : !! tolwfr=tolerance on largest wf residual
211 : : !! usedcwavef=flag controlling the use of dcwavef array (PAW only):
212 : : !! 0: not used (not allocated)
213 : : !! 1: used as input
214 : : !! 2: used as output
215 : : !! wfoptalg=govern the choice of algorithm for wf optimisation (0 or 10, at present)
216 : : !! [usetolrde]: 0 if the reduction in trial energy deltae should not be checked (see below). Default 1.
217 : : !! 0 is used in NSCF mode to prevent an early return before nline iterations.
218 : : !!
219 : : !! OUTPUT
220 : : !! eig1_k(2*nband**2)=matrix of first-order eigenvalues (hartree)
221 : : !! eig1(:,ii,jj)=<C0 ii|H1|C0 jj> for norm-conserving psps
222 : : !! eig1(:,ii,jj)=<C0 ii|H1-(eig0_k+eig0_k+q)/2.S(1)|C0 jj> for PAW
223 : : !! ghc(2,npw1*nspinor)=<G|H0-eig0_k.I|C1 band,k> (NCPP) or <G|H0-eig0_k.S0|C1 band,k> (PAW)
224 : : !! gvnlxc(2,npw1*nspinor)=<G|Vnl+VFockACE|C1 band,k>
225 : : !! gvnlx1(2,npw1*nspinor)= part of <G|K1+Vnl1|C0 band,k> not depending on VHxc1 (NCPP)
226 : : !! or part of <G|K1+Vnl1-eig0k.S1|C0 band,k> not depending on VHxc1 (PAW)
227 : : !! resid=wf residual for current band
228 : : !! gh1c_n= <G|H1|C0 band,k> (NCPP) or <G|H1-eig0k.S1|C0 band,k> (PAW).
229 : : !! This vector is not projected on the subspace orthogonal to the WF.
230 : : !! === if gs_hamkq%usepaw==1 ===
231 : : !! gsc(2,npw1*nspinor*usepaw)=<G|S0|C1 band,k>
232 : : !!
233 : : !! SIDE EFFECTS
234 : : !! Input/Output:
235 : : !! cwavef(2,npw1*nspinor)=first-order wavefunction at k,q, in reciprocal space (updated)
236 : : !!
237 : : !! === if gs_hamkq%usepaw==1 ===
238 : : !! cwaveprj(natom,nspinor)= wave functions at k projected with nl projectors
239 : : !!
240 : : !! === if also usedcwavef>0 ===
241 : : !! dcwavef(2,npw1*nspinor)=change of wavefunction due to change of overlap:
242 : : !! dcwavef is delta_Psi(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
243 : : !! see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)
244 : : !! input if usedcwavef=1, output if usedcwavef=2
245 : : !!
246 : : !! SOURCE
247 : :
248 : 4643781 : subroutine dfpt_cgwf(u1_band_,band_me,rank_band,bands_treated_now,berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,&
249 : 4643781 : & rf2,dcwavef,&
250 : 4643781 : & eig0_k,eig0_kq,eig1_k,ghc,gh1c_n,grad_berry,gsc,gscq,&
251 : 4643781 : & gs_hamkq,gvnlxc,gvnlx1,icgq,idir,ipert,igscq,&
252 : : & mcgq,mgscq,mpi_enreg,mpw1,natom,nband,nband_me,nbdbuf,nline_in,npw,npw1,nspinor,&
253 : : & opt_gvnlx1,prtvol,quit,resid,rf_hamkq,dfpt_sciss,tolrde,tolwfr,&
254 : : & usedcwavef,wfoptalg,nlines_done, &
255 : : eta,omega,usetolrde) ! optional
256 : :
257 : : !Arguments ------------------------------------
258 : : !scalars
259 : : integer,intent(in) :: u1_band_,berryopt
260 : : integer,intent(in) :: band_me, nband_me
261 : : integer,intent(in) :: icgq,idir,igscq,ipert,mcgq,mgscq,mpw1,natom,nband
262 : : integer,intent(in) :: nbdbuf,nline_in,npw,npw1,nspinor,opt_gvnlx1
263 : : integer,intent(in) :: prtvol,quit,usedcwavef,wfoptalg
264 : : integer,intent(inout) :: nlines_done
265 : : integer,optional,intent(in) :: usetolrde
266 : : real(dp),intent(in) :: dfpt_sciss,tolrde,tolwfr
267 : : real(dp),optional,intent(in) :: eta,omega
268 : : real(dp),intent(out) :: resid
269 : : type(MPI_type),intent(in) :: mpi_enreg
270 : : type(rf2_t), intent(in) :: rf2
271 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
272 : : type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
273 : : !arrays
274 : : integer,intent(in) :: rank_band(nband)
275 : : integer,intent(in) :: bands_treated_now (nband)
276 : : real(dp),intent(in) :: cgq(2,mcgq),eig0_kq(nband)
277 : : real(dp),intent(in) :: eig0_k(nband)
278 : : real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband),gscq(2,mgscq)
279 : : real(dp),intent(inout) :: cwave0(2,npw*nspinor),cwavef(2,npw1*nspinor)
280 : : real(dp),intent(inout) :: dcwavef(2,npw1*nspinor*((usedcwavef+1)/2))
281 : : real(dp),intent(inout) :: eig1_k(2*nband**2)
282 : : real(dp),intent(out) :: gh1c_n(2,npw1*nspinor)
283 : : real(dp),intent(out) :: ghc(2,npw1*nspinor)
284 : : real(dp),intent(out) :: gsc(2,npw1*nspinor*gs_hamkq%usepaw)
285 : : real(dp),intent(inout) :: gvnlx1(2,npw1*nspinor),gvnlxc(2,npw1*nspinor)
286 : : type(pawcprj_type),intent(inout) :: cwaveprj(natom,nspinor)
287 : : type(pawcprj_type),intent(inout) :: cwaveprj0(natom,nspinor*gs_hamkq%usecprj)
288 : :
289 : : !Local variables-------------------------------
290 : : !scalars
291 : : integer,parameter :: level=15,tim_getgh1c=1,tim_getghc=2,tim_projbd=2
292 : : integer,save :: nskip=0
293 : : integer :: cpopt,iband,igs,iline,indx_cgq,ipw,me_g0,comm_fft
294 : : integer :: iband_me, jband_me, ierr, me_band, np_band, band_off, u1_band
295 : : integer :: ipws,ispinor,istwf_k,jband,nline,optlocal,optnl,dc_shift_band,sij_opt
296 : : integer :: test_is_ok,useoverlap,usepaw,usevnl,usetolrde__
297 : : real(dp) :: d2edt2,d2te,d2teold,dedt,deltae,deold,dotgg
298 : : real(dp) :: dotgp,doti,dotr,eta_,eshift,eshiftkq,gamma,omega_,optekin,prod1,prod2
299 : : real(dp) :: theta,tol_restart,u1h0me0u1
300 : : logical :: gen_eigenpb
301 : 9287562 : integer :: skipme, bands_skipped_now(nband)
302 : : character(len=500) :: msg
303 : : !arrays
304 : : real(dp) :: dummy(0,0),tsec(2)
305 : 9287562 : real(dp) :: eig1_k_loc(2,nband,nband)
306 : 9287562 : real(dp),allocatable :: conjgr(:,:),cwaveq(:,:),cwwork(:,:),direc(:,:)
307 : 4643781 : real(dp),allocatable :: gberry(:,:),gh1c(:,:),gh_direc(:,:),gresid(:,:),gvnlx1_saved(:,:)
308 : 4643781 : real(dp),allocatable :: gs1c(:,:),gvnlx_direc(:,:),pcon(:),sconjgr(:,:)
309 : 4643781 : real(dp),allocatable :: scprod(:,:),work(:,:),work1(:,:),work2(:,:)
310 : 4643781 : real(dp),pointer :: kinpw1(:)
311 : 4643781 : type(pawcprj_type),allocatable :: conjgrprj(:,:)
312 [ + + + + : 32506467 : type(pawcprj_type) :: cprj_dummy(1,1)
+ + + - +
- ]
313 : : ! *********************************************************************
314 : :
315 : : DBG_ENTER("COLL")
316 : :
317 : : ABI_NVTX_START_RANGE(NVTX_DFPT_CGWF)
318 : :
319 : 4643781 : call timab(122,1,tsec)
320 : :
321 : : !==================================================================
322 : : !========= LOCAL VARIABLES DEFINITIONS AND ALLOCATIONS ============
323 : : !==================================================================
324 : :
325 : 4643781 : u1_band = abs(u1_band_)
326 : :
327 : : ! LB-23/04/17:
328 : : ! For ipert=natom+10 or ipert=natom+11, the Sternheimer equation is non-self-consistent, so we have
329 : : ! to solve a true linear problem (A.X = B) for each kpoint and band. In this case, the conjugate
330 : : ! gradient algorithm can find the exact solution (within the numerical precision) with only ONE call
331 : : ! of dfpt_cgwf (per kpoint and band...). The solution is found with at most N iterations, N being the dimension of X.
332 : : ! In order to avoid useless scfcv loops (and many calls of rf2_init, which can be time consuming),
333 : : ! we want to leave this routine only if 'tolwfr' is reached. Consequently, 'tolrde' is not used and 'nline' is set to 100.
334 : : ! One could use nline=npw1*nspinor (>> 100 !) instead, but when the method cannot converge (i.e when tolwfr is lower than the
335 : : ! numerical noise) the program could be stuck here for a very long time.
336 : : ! NOTE: This is also true for ipert==natom+1, but a lot of references in the test suite have to be changed...
337 : 4643781 : nline = nline_in
338 : 4643781 : usetolrde__ = 1
339 [ + + + + ]: 4643781 : if(ipert==natom+10.or.ipert==natom+11) then
340 : 59712 : nline = 100 ! The default value is only 4... This should be sufficient to converge with nstep=1 or 2
341 [ - + ]: 59712 : if (nline_in > 100) nline = nline_in ! Keep the possibility to increase nline
342 : : usetolrde__ = 0 ! see below
343 : : end if
344 [ + + ]: 4643781 : if (present(usetolrde)) usetolrde__ = usetolrde
345 : :
346 [ + + ]: 4643781 : if (prtvol>=10) then
347 : : !Tell us what is going on:
348 : 278931 : write(msg,'(a,i0,2x,a,i0,a)')' --- dfpt_cgwf is called for band: ', u1_band,' for: ',nline,' lines'
349 : 278931 : call wrtout(std_out,msg)
350 : : end if
351 : :
352 : 4643781 : me_g0 = mpi_enreg%me_g0
353 : 4643781 : comm_fft = mpi_enreg%comm_fft
354 : 4643781 : me_band = mpi_enreg%me_band
355 : 4643781 : np_band = mpi_enreg%nproc_band
356 : :
357 : 4643781 : skipme = 0
358 : :
359 : : ! if PAW, one has to solve a generalized eigenproblem
360 : 4643781 : usepaw=gs_hamkq%usepaw
361 : 4643781 : gen_eigenpb=(usepaw==1)
362 [ + + ]: 4643781 : useoverlap=0;if (gen_eigenpb) useoverlap=1
363 : :
364 : : ! Use scissor shift on 0-order eigenvalue
365 : 4643781 : eshift=eig0_k(u1_band)-dfpt_sciss
366 : :
367 : : ! Remove omega for a finite-frequency calculation
368 [ + + ]: 4643781 : omega_=zero ; if (present(omega)) omega_=omega
369 [ + + ]: 4643781 : eta_=zero ; if (present(eta)) eta_=eta
370 : 4643781 : eshift=eshift+omega_
371 : :
372 : : ! Additional initializations
373 : 4643781 : istwf_k=gs_hamkq%istwf_k
374 : 4643781 : optekin=0;if (wfoptalg>=10) optekin=1
375 [ + + ]: 4643781 : tol_restart=tol12;if (gen_eigenpb) tol_restart=tol8
376 [ + + + + ]: 4643781 : if (ipert == natom+10 .or. ipert == natom+11) tol_restart = tol7
377 : : !tol_restart = tol7
378 : :
379 : 4643781 : kinpw1 => gs_hamkq%kinpw_kp
380 : :
381 : : ! Memory allocations
382 [ - + + - : 13931343 : ABI_MALLOC(gh1c,(2,npw1*nspinor))
+ - + - -
+ - + ]
383 [ - + + - : 13931343 : ABI_MALLOC(pcon,(npw1))
+ - + - -
+ - + ]
384 [ - + + - : 13931343 : ABI_MALLOC(scprod,(2,nband_me))
+ - + - -
+ - + ]
385 : :
386 [ + + ]: 4643781 : if (berryopt== 4.or.berryopt== 6.or.berryopt== 7.or. berryopt==14.or.berryopt==16.or.berryopt==17) then
387 [ - + + - : 6208 : ABI_MALLOC(gberry,(2,npw1*nspinor))
- + - + ]
388 [ + + + + ]: 389552 : gberry(:,1:npw1*nspinor)=grad_berry(:,1:npw1*nspinor,u1_band)
389 : : else
390 [ - + ]: 4640677 : ABI_MALLOC(gberry,(0,0))
391 : : end if
392 : :
393 : : !TODO MJV: this should probably be adjusted as well for natom+10/11 perts: set to band_me instead of u1_band
394 : 4643781 : dc_shift_band=(u1_band-1)*npw1*nspinor
395 : :
396 : : ! this is used many times - no use de and re allocating
397 [ - + + - : 9287562 : ABI_MALLOC(work,(2,npw1*nspinor))
- + - + -
+ ]
398 : : #ifdef HAVE_OPENMP_OFFLOAD
399 : : !$OMP TARGET ENTER DATA MAP(alloc:work,scprod) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
400 : : #endif
401 : :
402 : : !DEBUG!! Several checking statements
403 [ + - - + ]: 4643781 : if (prtvol==-level.or.prtvol==-19.or.prtvol==-20) then
404 : 0 : write(msg,'(a)') " ** cgwf3 : debugging mode, tests will be done"
405 : : ! Search CGWF3_WARNING in the log file to find errors (if any)
406 : 0 : call wrtout(std_out,msg)
407 [ # # # # : 0 : ABI_MALLOC(work1,(2,npw1*nspinor))
# # # # #
# ]
408 : : ! ===== Check <Psi_k+q^(0)|S(0)|Psi_k+q^(0)>=delta_{ij}
409 [ # # # # : 0 : if (.not.gen_eigenpb) work1(:,:)=cgq(:,1+npw1*nspinor*(band_me-1)+icgq:npw1*nspinor*band_me+icgq)
# # ]
410 [ # # # # ]: 0 : if ( gen_eigenpb) work1(:,:)=gscq(:,1+npw1*nspinor*(band_me-1)+igscq:npw1*nspinor*band_me+igscq)
411 : : ! NB: this loop is not band-block diagonal
412 : : ! the present logic does a lot of communication: 1 for each band.
413 : : ! Could be grouped outside the jband loop into 1 big one, but if the wf are big this is a waste. Tradeoffs...
414 : 0 : jband_me = 0
415 [ # # ]: 0 : do jband=1,nband
416 [ # # ]: 0 : if (bands_treated_now(jband) == 0) cycle
417 [ # # ]: 0 : if (rank_band(jband) == me_band) then
418 : 0 : jband_me = jband_me + 1
419 [ # # # # ]: 0 : work(:,:)=cgq(:,1+npw1*nspinor*(jband_me-1)+icgq:npw1*nspinor*jband_me+icgq)
420 : : end if
421 : : ! send to everyone else, who is also working on jband right now
422 : 0 : call xmpi_bcast(work,rank_band(jband),mpi_enreg%comm_band,ierr)
423 : :
424 : 0 : call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,work1,work,me_g0,mpi_enreg%comm_spinorfft)
425 : 0 : test_is_ok=1
426 [ # # ]: 0 : if(jband==u1_band) then
427 [ # # ]: 0 : if(abs(dotr-one)>tol12) test_is_ok=0
428 : : else
429 [ # # ]: 0 : if(abs(dotr)>tol12) test_is_ok=0
430 : : end if
431 [ # # ]: 0 : if(abs(doti)>tol12) test_is_ok=0
432 [ # # ]: 0 : if(test_is_ok/=1) then
433 : 0 : write(msg,'(a,i3,a,2es22.15)') "CGWF3_WARNING : <Psi_k+q,i^(0)|S(0)|Psi_k+q,j^(0)> for band j=",jband," is ",dotr,doti
434 : 0 : call wrtout(std_out,msg)
435 : : end if
436 : : end do
437 : :
438 : : ! ===== Check Pc.Psi_k+q^(0)=0
439 : : ! each jband is checked by everybody, against the bands attributed to present cpu
440 : : ! NB - this does not depend on the "band" input
441 : 0 : jband_me = 0
442 [ # # ]: 0 : do jband=1,nband
443 [ # # ]: 0 : if (bands_treated_now(jband) == 0) cycle
444 [ # # ]: 0 : if (rank_band(jband) == me_band) then
445 : 0 : jband_me = jband_me + 1
446 [ # # # # ]: 0 : work(:,:)=cgq(:,1+npw1*nspinor*(jband_me-1)+icgq:npw1*nspinor*jband_me+icgq)
447 : : end if
448 : : ! send to everyone else, who is also working on jband right now
449 : 0 : call xmpi_bcast(work,rank_band(jband),mpi_enreg%comm_band,ierr)
450 [ # # # # : 0 : work1 = work
# # # # #
# # # # #
# # ]
451 : :
452 : : call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
453 : 0 : gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
454 : :
455 : : ! if bands are parallelized, I have only projected against bands on my cpu
456 : : ! Pc|work> = |work> - Sum_l <psi_{k+q, l}|work> |psi_{k+q, l}>
457 : : ! = Sum_nproc_band (|work> - Sum_{my l} <psi_{k+q, l}|work> |psi_{k+q, l}>) - (nproc_band-1) |work>
458 : : !
459 [ # # ]: 0 : if (mpi_enreg%nproc_band > 1) then
460 : 0 : call xmpi_sum(work,mpi_enreg%comm_band,ierr)
461 : : !TODO: make this a blas call? zaxpy
462 [ # # # # : 0 : work = work - (mpi_enreg%nproc_band-1)*work1
# # ]
463 : : end if
464 : :
465 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
466 [ # # ]: 0 : if(sqrt(dotr)>tol12) then
467 : 0 : write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc.Psi_k+q_j^(0) for band j=",jband," is ",sqrt(dotr)
468 : 0 : call wrtout(std_out,msg)
469 : : end if
470 : : end do
471 : :
472 : : ! ===== Check Pc.Psi_k^(0)=0
473 : : ! NB: this _does_ depend on the input band "band" stored in cwave0
474 [ # # ]: 0 : do iband = 1, nband
475 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
476 [ # # # # : 0 : if (rank_band(iband) == me_band) work(:,:)=cwave0(:,:)
# # ]
477 : :
478 : : ! send to everyone else, who is also working on jband right now
479 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
480 [ # # # # : 0 : work1 = work
# # # # #
# # # # #
# # ]
481 : : call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
482 : 0 : gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
483 : :
484 [ # # ]: 0 : if (mpi_enreg%nproc_band > 1) then
485 : 0 : call xmpi_sum(work,mpi_enreg%comm_band,ierr)
486 : : !TODO: make this a blas call? zaxpy
487 [ # # # # : 0 : work = work - (mpi_enreg%nproc_band-1)*work1
# # ]
488 : : end if
489 : :
490 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
491 [ # # ]: 0 : if(sqrt(dotr)>tol12) then
492 : 0 : write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc.Psi_k^(0) for band ",u1_band," is ",sqrt(dotr)
493 : 0 : call wrtout(std_out,msg)
494 : : end if
495 : : end do
496 : :
497 : : ! ===== Check Pc.dcwavef=0 (for 2nd order only)
498 [ # # # # ]: 0 : if(ipert==natom+10.or.ipert==natom+11) then
499 [ # # ]: 0 : do iband = 1, nband
500 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
501 [ # # ]: 0 : if (rank_band(iband) == me_band) then
502 [ # # # # ]: 0 : work(:,:)=rf2%dcwavef(:,1+dc_shift_band:npw1*nspinor+dc_shift_band)
503 : : end if
504 : : ! send to everyone else, who is also working on jband right now
505 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
506 [ # # # # : 0 : work1 = work
# # # # #
# # # # #
# # ]
507 : : call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
508 : 0 : gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
509 : :
510 [ # # ]: 0 : if (mpi_enreg%nproc_band > 1) then
511 : 0 : call xmpi_sum(work,mpi_enreg%comm_band,ierr)
512 : : !TODO: make this a blas call? zaxpy
513 [ # # # # : 0 : work = work - (mpi_enreg%nproc_band-1)*work1
# # ]
514 : : end if
515 : :
516 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
517 [ # # ]: 0 : if(sqrt(dotr)>tol10) then
518 : 0 : write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc.dcwavef for band ",u1_band," is ",sqrt(dotr)
519 : 0 : call wrtout(std_out,msg)
520 : : end if
521 : : end do
522 : : end if
523 : :
524 : : ! ===== Check Pc^*.S(0).Psi_k+q^(0)=0
525 : : ! NB: here again does not depend on input "u1_band"
526 [ # # ]: 0 : if (gen_eigenpb) then
527 : 0 : jband_me = 0
528 [ # # ]: 0 : do jband=1,nband
529 [ # # ]: 0 : if (bands_treated_now(jband) == 0) cycle
530 [ # # ]: 0 : if (rank_band(jband) == me_band) then
531 : 0 : jband_me = jband_me + 1
532 [ # # # # ]: 0 : work(:,:)=gscq(:,1+npw1*nspinor*(jband_me-1)+igscq:npw1*nspinor*jband_me+igscq)
533 : : end if
534 : : ! send to everyone else, who is also working on jband right now
535 : 0 : call xmpi_bcast(work,rank_band(jband),mpi_enreg%comm_band,ierr)
536 [ # # # # : 0 : work1 = work
# # # # #
# # # # #
# # ]
537 : :
538 : : call projbd(gscq,work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
539 : 0 : cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
540 : :
541 [ # # ]: 0 : if (mpi_enreg%nproc_band > 1) then
542 : 0 : call xmpi_sum(work,mpi_enreg%comm_band,ierr)
543 : : !TODO: make this a blas call? zaxpy
544 [ # # # # : 0 : work = work - (mpi_enreg%nproc_band-1)*work1
# # ]
545 : : end if
546 : :
547 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
548 [ # # ]: 0 : if(sqrt(dotr)>tol12) then
549 : 0 : write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc^*.S(0).Psi_k+q_j^(0) for band j=",jband," is ",sqrt(dotr)
550 : 0 : call wrtout(std_out,msg)
551 : : end if
552 : : end do
553 : : end if
554 [ # # ]: 0 : ABI_FREE(work1)
555 : : end if
556 : : !ENDDEBUG!! Several checking statements
557 : :
558 : :
559 : : !======================================================================
560 : : !========== INITIALISATION OF MINIMIZATION ITERATIONS =================
561 : : !======================================================================
562 : :
563 [ + + + + ]: 4643781 : if (ipert/=natom+10.and.ipert/=natom+11) then
564 : : ! The following is needed for first order perturbations only
565 : : ! Otherwise, the work is already done in rf2_init (called in dfpt_vtowfk.F90)
566 : :
567 : : ! Compute H(1) applied to GS wavefunction Psi(0)
568 [ + + ]: 4584069 : if (gen_eigenpb) then
569 : 604076 : sij_opt=1
570 [ - + + - : 1208152 : ABI_MALLOC(gs1c,(2,npw1*nspinor))
- + - + ]
571 : : else
572 [ - + ]: 3979993 : ABI_MALLOC(gs1c,(0,0))
573 : 3979993 : sij_opt=0
574 : : end if
575 : 4584069 : usevnl=1; optlocal=1; optnl=2
576 [ - + ]: 4584069 : if (prtvol==-level.or.prtvol==-19) then
577 [ # # # # : 0 : ABI_MALLOC(gvnlx1_saved,(2,npw1*nspinor))
# # # # ]
578 [ # # # # ]: 0 : gvnlx1_saved(:,:) = gvnlx1(:,:)
579 : : end if
580 : : call getgh1c(berryopt,cwave0,cwaveprj0,gh1c,gberry,gs1c,gs_hamkq,gvnlx1,idir,ipert,[eshift],&
581 [ + + ]: 9168138 : mpi_enreg,1,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
582 : :
583 [ + + ]: 4584069 : if (gen_eigenpb) then
584 [ + + ]: 604076 : if (ipert/=natom+2) then ! S^(1) is zero for ipert=natom+2
585 : : !$OMP PARALLEL
586 : : !$OMP DO
587 [ + + ]: 83342118 : do ipw=1,npw1*nspinor
588 [ + + ]: 248927962 : gh1c (1:2,ipw)=gh1c (1:2,ipw)-eshift*gs1c(1:2,ipw)
589 : : end do
590 : : !$OMP END DO NOWAIT
591 [ + - ]: 549196 : if (opt_gvnlx1/=1) then
592 : : !$OMP DO
593 [ + + ]: 83342118 : do ipw=1,npw1*nspinor
594 [ + + ]: 248927962 : gvnlx1(1:2,ipw)=gvnlx1(1:2,ipw)-eshift*gs1c(1:2,ipw)
595 : : end do
596 : : !$OMP END DO NOWAIT
597 : : end if
598 : : !$OMP END PARALLEL
599 : : end if
600 : :
601 : : ! If generalized eigenPb and dcwavef requested, compute it:
602 : : ! dcwavef is delta_Psi(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
603 : : ! see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)
604 [ + + ]: 604076 : if (usedcwavef==2) then
605 : : call getdc1(u1_band,rank_band,bands_treated_now,cgq,cprj_dummy,dcwavef,cprj_dummy,&
606 : : & 0,icgq,istwf_k,mcgq,0,&
607 : 549196 : & mpi_enreg,1,natom,nband,nband_me,npw1,nspinor,0,gs1c,gpu_option=ABI_GPU_DISABLED)
608 : : end if
609 : : end if ! gen_eigenpb
610 : :
611 : : else
612 : : ! 2nd order case (wrt k perturbation)
613 : : ! Copy RHS_Stern(:,:) of the given band in gh1c
614 [ + + + + ]: 22080132 : gh1c(:,:)=rf2%RHS_Stern(:,1+dc_shift_band:npw1*nspinor+dc_shift_band)
615 : : end if
616 : :
617 [ - + - - ]: 4643781 : if (prtvol==-level.and.usedcwavef==2) then
618 : : !Check that Pc^*.(H^(0)-E.S^(0)).delta_Psi^(1) is zero ! This is a consequence of P_c delta_Psi^(1) = 0
619 [ # # # # : 0 : ABI_MALLOC(cwwork,(2,npw1*nspinor))
# # # # #
# ]
620 [ # # ]: 0 : do iband = 1, nband
621 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
622 [ # # ]: 0 : if (rank_band(iband) == me_band) then
623 [ # # # # : 0 : cwwork=dcwavef
# # # # #
# # # ]
624 : : ! - Apply H^(0)-E.S^(0) to delta_Psi^(1)
625 [ # # ]: 0 : sij_opt=0;if (gen_eigenpb) sij_opt=-1
626 : 0 : cpopt=-1
627 [ # # # # : 0 : ABI_MALLOC(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
# # # # #
# # # #
# ]
628 [ # # # # : 0 : ABI_MALLOC(work2,(2,npw1*nspinor))
# # # # #
# ]
629 : : call getghc(cpopt,cwwork,conjgrprj,work,work1,gs_hamkq,work2,eshift,mpi_enreg,&
630 : 0 : 1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
631 [ # # ]: 0 : ABI_FREE(work1)
632 [ # # ]: 0 : ABI_FREE(work2)
633 : : end if
634 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
635 [ # # # # : 0 : cwwork=work
# # # # #
# # # # #
# # ]
636 : :
637 : : ! -Apply Pc^*
638 : : call projbd(gscq,cwwork,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
639 : 0 : cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
640 : :
641 : 0 : call xmpi_sum(cwwork,mpi_enreg%comm_band,ierr)
642 : :
643 [ # # ]: 0 : if (mpi_enreg%nproc_band > 1) then
644 : : !TODO: make this a blas call? zaxpy
645 [ # # # # : 0 : cwwork = cwwork - (mpi_enreg%nproc_band-1)*work
# # ]
646 : : end if
647 : :
648 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwwork,me_g0,comm_fft)
649 [ # # ]: 0 : if(sqrt(dotr)>tol12) then
650 : 0 : write(msg,'(a,i3,a,es22.15)') 'CGWF3_WARNING : |Pc^*.(H^(0)-E.S^(0)).delta_Psi^(1)| (band ',u1_band,')=',sqrt(dotr)
651 : 0 : call wrtout(std_out,msg)
652 : : end if
653 : : end do
654 [ # # ]: 0 : ABI_FREE(cwwork)
655 : : end if ! prtvol==-level.and.usedcwavef==2
656 : :
657 : 4643781 : call cg_zcopy(npw1*nspinor,gh1c,gh1c_n)
658 : :
659 : : ! Projecting out all bands
660 : : ! While we could avoid calculating all the eig1_k to obtain the perturbed density,
661 : : ! we do need all of the matrix elements when outputting the full 1st-order wfn.
662 : : ! Note the subtlety:
663 : : ! -For the generalized eigenPb, S|cgq> is used in place of |cgq>,
664 : : ! in order to apply P_c+ projector (see PRB 73, 235101 (2006) [[cite:Audouze2006]], Eq. (71), (72))
665 [ + + + + : 1756034683 : eig1_k_loc = zero
+ + ]
666 [ + + ]: 41606779 : do iband = 1, nband
667 [ + + ]: 36962998 : if (bands_treated_now(iband) == 0) cycle
668 [ + + + - : 3166586520 : if (rank_band(iband) == me_band) work = gh1c
+ - - + -
- + + +
+ ]
669 : 4650045 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
670 : :
671 : : #ifdef HAVE_OPENMP_OFFLOAD
672 : : !$OMP TARGET UPDATE TO (work) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
673 : : #endif
674 [ + + ]: 4650045 : if(gen_eigenpb)then
675 : : call projbd(gscq,work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
676 : 611300 : cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft,gpu_option=gs_hamkq%gpu_option)
677 : : else
678 : : call projbd(cgq,work,-1,icgq,0,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
679 : 4038745 : dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft,gpu_option=gs_hamkq%gpu_option)
680 : : end if
681 : : #ifdef HAVE_OPENMP_OFFLOAD
682 : : !$OMP TARGET UPDATE FROM (work,scprod) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
683 : : #endif
684 : :
685 : : ! sum projections against all bands k+q
686 : 4650045 : call xmpi_sum_master(work, rank_band(iband), mpi_enreg%comm_band, ierr)
687 : :
688 : : ! scprod now contains scalar products of band iband (runs over all bands in current queue) with local bands j
689 : 4650045 : jband_me = 0
690 [ + + ]: 41759059 : do jband=1,nband
691 [ + + ]: 37109014 : if (rank_band(jband) /= me_band) cycle
692 : 36962998 : jband_me = jband_me + 1
693 [ + + ]: 115539039 : eig1_k_loc(:,jband,iband)=scprod(:,jband_me)
694 : : end do
695 : :
696 : : ! save this for me only
697 : : !TODO: make this a blas call? zaxpy
698 [ + + ]: 18593916 : if (rank_band(iband) == me_band) then
699 [ + - + - : 3166580256 : gh1c = work - (mpi_enreg%nproc_band-1)*gh1c
- + - - +
+ + + ]
700 : : end if
701 : :
702 : : end do !iband
703 : :
704 : :
705 [ + + + + ]: 4643781 : if(ipert/=natom+10.and.ipert/=natom+11) then
706 : : ! For ipert=natom+10 or natom+11, this is done in rf2_init
707 : :
708 : : ! The array eig1_k contains:
709 : : ! <u_(jband,k+q)^(0)|H_(k+q,k)^(1)|u_(iband,k)^(0)> (NC psps)
710 : : ! <u_(jband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(iband,k)^(0)> (PAW)
711 : : ! so in case of PAW need to add the overlap term below
712 : : !
713 : : ! NB: 2019 11 15: MJV: I swapped the names of jband and iband to be more consistent with other loops above
714 [ + + ]: 4584069 : if (gen_eigenpb) then
715 [ + + ]: 5073012 : do iband=1,nband
716 [ + + ]: 4468936 : if (bands_treated_now(iband) == 0) cycle
717 [ + + + - : 272376670 : if (rank_band(iband) == me_band) work = gs1c
+ - - + -
- + + +
+ ]
718 : : ! for iband on this proc, bcast to all others to get full line of iband,jband pairs
719 : 607460 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
720 : :
721 : : ! add PAW overlap correction term to present iband (all procs) and local jband elements
722 : 607460 : indx_cgq=icgq
723 [ + + ]: 5757368 : do jband=1,nband
724 [ + + ]: 4545832 : if (rank_band(jband) /= me_band) cycle
725 : :
726 : 4468936 : eshiftkq=half*(eig0_kq(jband)-eig0_k(iband))
727 : : call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,cgq(:,indx_cgq+1:indx_cgq+npw1*nspinor),work,&
728 : 4468936 : me_g0,mpi_enreg%comm_spinorfft)
729 : 4468936 : eig1_k_loc(1,jband,iband)=eig1_k_loc(1,jband,iband)-eshiftkq*dotr
730 : 4468936 : eig1_k_loc(2,jband,iband)=eig1_k_loc(2,jband,iband)-eshiftkq*doti
731 : 9014768 : indx_cgq=indx_cgq+npw1*nspinor
732 : : end do
733 : : end do ! iband
734 : : end if ! PAW and generalized eigenproblem
735 : :
736 : : ! No more need of gs1c
737 [ - + ]: 4584069 : ABI_FREE(gs1c)
738 : :
739 : : ! reduce over band procs to fill in the matrix for all jband (distributed over procs)
740 : : ! must only do this once for eig1_k_loc: now have all jband for current ibands on all procs
741 : 4584069 : call xmpi_sum(eig1_k_loc, mpi_enreg%comm_band, ierr)
742 : :
743 : : ! TODO: I think this is just a reshape
744 [ + + ]: 41253819 : do iband=1,nband
745 [ + + ]: 36669750 : if (bands_treated_now(iband) == 0) cycle
746 : 4590333 : band_off=(iband-1)*2*nband
747 [ + + ]: 45990168 : do jband=1,nband
748 : 36815766 : eig1_k(2*jband-1+band_off) = eig1_k_loc(1,jband,iband)
749 : 73485516 : eig1_k(2*jband +band_off) = eig1_k_loc(2,jband,iband)
750 : : end do
751 : : end do
752 : : end if ! ipert/=natom+10.and.ipert/=natom+11
753 : :
754 : : ! Filter the wavefunctions for large modified kinetic energy (see routine mkkin.f)
755 : : ! TODO: should this also be applied to cwaveq for the preconditioning with kinpw1 below?
756 [ + + ]: 9431000 : do ispinor=1,nspinor
757 : 4787219 : ipws=(ispinor-1)*npw1
758 : : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef,kinpw1,ipws,npw1)
759 [ + + ]: 1061861898 : do ipw=1+ipws,npw1+ipws
760 [ - + ]: 1057218117 : if(kinpw1(ipw-ipws)>huge(zero)*1.d-11)then
761 [ # # ]: 0 : cwavef(1:2,ipw)=zero
762 : : end if
763 : : end do
764 : : end do
765 : :
766 : : ! Apply the orthogonality condition: <C1 k,q|C0 k+q>=0 (NCPP) or <C1 k,q|S0|C0 k+q>=0 (PAW)
767 : : ! Project out all bands from cwavef, i.e. apply P_c projector on cwavef
768 : : ! (this is needed when there are some partially or unoccupied states)
769 [ + + ]: 41606779 : do iband = 1, nband
770 [ + + ]: 36962998 : if (bands_treated_now(iband) == 0) cycle
771 : :
772 [ + + + - : 3166586520 : if (rank_band(iband) == me_band) work = cwavef
+ - - + -
- + + +
+ ]
773 : 4650045 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
774 : :
775 : : #ifdef HAVE_OPENMP_OFFLOAD
776 : : !$OMP TARGET UPDATE TO (work) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
777 : : #endif
778 : : call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
779 : 4650045 : gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft,gpu_option=gs_hamkq%gpu_option)
780 : : #ifdef HAVE_OPENMP_OFFLOAD
781 : : !$OMP TARGET UPDATE FROM (work,scprod) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
782 : : #endif
783 : :
784 : 4650045 : call xmpi_sum_master(work, rank_band(iband), mpi_enreg%comm_band, ierr)
785 : :
786 : : ! save this for me_band only
787 : : !TODO: make this a blas call? zaxpy
788 [ + + ]: 13943871 : if (rank_band(iband) == me_band) then
789 [ + + + + ]: 3161936475 : cwavef = work - (mpi_enreg%nproc_band-1)*cwavef
790 : : end if
791 : : end do
792 : :
793 : :
794 [ + + + + ]: 4643781 : if(ipert/=natom+10.and.ipert/=natom+11) then
795 : : ! If PAW, the orthogonality condition is <C1 k,q|S0|C0 k+q>+1/2<C0 k|S1|C0 k+q>=0
796 [ + + + + ]: 4584069 : if (usepaw==1.and.usedcwavef>0) then
797 : : !$OMP PARALLEL DO
798 [ + + ]: 83342118 : do ipw=1,npw1*nspinor
799 [ + + ]: 248927962 : cwavef(1:2,ipw)=cwavef(1:2,ipw)+dcwavef(1:2,ipw)
800 : : end do
801 : : end if
802 : : else
803 : : ! In 2nd order case, dcwavef/=0 even in NC, and it is already computed in rf2_init (called in dfpt_vtowfk.F90)
804 [ + + ]: 7399852 : do ipw=1,npw1*nspinor
805 [ + + ]: 22080132 : cwavef(:,ipw)=cwavef(:,ipw)+rf2%dcwavef(:,ipw+dc_shift_band)
806 : : end do
807 : : end if
808 : :
809 [ + + ]: 4643781 : if (u1_band>max(1,nband-nbdbuf))then
810 : : ! Treat the case of buffer bands
811 [ + + + + ]: 280460796 : cwavef=zero
812 [ + + + + ]: 280460796 : ghc =zero
813 [ + + + + ]: 280460796 : gvnlxc =zero
814 [ + + + + : 23598714 : if (gen_eigenpb) gsc=zero
+ + ]
815 [ + + + + : 21721914 : if (usedcwavef==2) dcwavef=zero
+ + ]
816 [ + + ]: 398682 : if (usepaw==1) then
817 : 26932 : call pawcprj_set_zero(cwaveprj)
818 : : end if
819 [ + + + + : 21721914 : if (usedcwavef==2) dcwavef=zero
+ + ]
820 : : ! Number of one-way 3D ffts skipped
821 : 398682 : nskip=nskip+nline
822 : :
823 : : ! At the end of the treatment of a set of bands, write the number of one-way 3D ffts skipped
824 [ + + + + ]: 398682 : if (xmpi_paral==1 .and. u1_band==nband .and. prtvol>=10) then
825 : 41857 : write(msg,'(a,i0)')' dfpt_cgwf: number of one-way 3D ffts skipped in cgwf3 until now =',nskip
826 : 41857 : call wrtout(std_out,msg)
827 : : end if
828 : :
829 : : skipme = 1
830 : : end if
831 : :
832 : : ! If not a buffer band, perform the optimisation
833 : :
834 [ - + + - : 9287562 : ABI_MALLOC(conjgr,(2,npw1*nspinor))
- + - + -
+ ]
835 [ - + + - : 9287562 : ABI_MALLOC(direc,(2,npw1*nspinor))
- + - + ]
836 [ - + + - : 9287562 : ABI_MALLOC(gresid,(2,npw1*nspinor))
- + - + ]
837 [ - + + - : 9287562 : ABI_MALLOC(cwaveq,(2,npw1*nspinor))
- + - + ]
838 [ + + ]: 4643781 : if (usepaw==1) then
839 [ - + + - : 4417912 : ABI_MALLOC(conjgrprj,(natom,nspinor))
- + + - +
- + - - +
- + - + +
+ + + ]
840 : 607916 : call pawcprj_alloc(conjgrprj,0,gs_hamkq%dimcprj)
841 : : else
842 [ - + - + ]: 4035865 : ABI_MALLOC(conjgrprj,(0,0))
843 : : end if
844 : :
845 [ + + + + ]: 3161936475 : cwaveq(:,:)=cgq(:,1+npw1*nspinor*(band_me-1)+icgq:npw1*nspinor*band_me+icgq)
846 : 4643781 : dotgp=one
847 : :
848 : : ! Here apply H(0) at k+q to input orthogonalized 1st-order wfs
849 [ + + ]: 4643781 : sij_opt=0;if (gen_eigenpb) sij_opt=1
850 : 4643781 : cpopt=-1+usepaw
851 : : call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamkq,gvnlxc,eshift,mpi_enreg,1,&
852 : 4643781 : prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
853 : :
854 : : ! ghc also includes the eigenvalue shift
855 [ + + ]: 4643781 : if (gen_eigenpb) then
856 [ + + ]: 1823748 : call cg_zaxpy(npw1*nspinor, [-eshift, -eta_], gsc,ghc)
857 : : else
858 [ + + ]: 12107595 : call cg_zaxpy(npw1*nspinor, [-eshift, -eta_], cwavef,ghc)
859 : : end if
860 : :
861 : : ! Initialize resid, in case of nline==0
862 : 4643781 : resid=zero
863 : :
864 [ + + ]: 41606779 : bands_skipped_now = 0
865 : 4643781 : bands_skipped_now(u1_band) = skipme
866 : 4643781 : call xmpi_sum(bands_skipped_now,mpi_enreg%comm_band,ierr)
867 : :
868 : : ! ======================================================================
869 : : ! ====== BEGIN LOOP FOR A GIVEN BAND: MINIMIZATION ITERATIONS ==========
870 : : ! ======================================================================
871 : : ABI_NVTX_START_RANGE(NVTX_DFPT_CGWF_CORE)
872 : :
873 [ + + ]: 17308209 : do iline=1,nline
874 : : ! ======================================================================
875 : : ! ================= COMPUTE THE RESIDUAL ===============================
876 : : ! ======================================================================
877 : : ! Note that gresid (=steepest-descent vector, Eq.(26) of PRB 55, 10337 (1996) [[cite:Gonze1997]])
878 : : ! is precomputed to guarantee cancellation of errors
879 : : ! and allow residuals to reach values as small as 1.0d-24 or better.
880 [ + + ]: 15713133 : if (berryopt== 4.or.berryopt== 6.or.berryopt== 7.or. berryopt==14.or.berryopt==16.or.berryopt==17) then
881 [ + + ]: 9715 : if (ipert==natom+2) then
882 [ - + - - : 3625 : if (opt_gvnlx1/=1) gvnlx1=zero
- - ]
883 : : !$OMP PARALLEL DO
884 [ + + ]: 154236 : do ipw=1,npw1*nspinor
885 [ + + ]: 455458 : gresid(1:2,ipw)=-ghc(1:2,ipw)-gh1c(1:2,ipw)
886 : : end do
887 : : else
888 : : !$OMP PARALLEL DO
889 [ + + ]: 259094 : do ipw=1,npw1*nspinor
890 : 253004 : gresid(1,ipw)=-ghc(1,ipw)-gh1c(1,ipw)+gberry(2,ipw)
891 : 259094 : gresid(2,ipw)=-ghc(2,ipw)-gh1c(2,ipw)-gberry(1,ipw)
892 : : end do
893 : : end if
894 : : else
895 : : !$OMP PARALLEL DO
896 [ + + ]: 4485092416 : do ipw=1,npw1*nspinor
897 [ + + ]: 13423870412 : gresid(1:2,ipw)=-ghc(1:2,ipw)-gh1c(1:2,ipw)
898 : : end do
899 : : end if
900 : :
901 : : ! ======================================================================
902 : : ! =========== PROJECT THE STEEPEST DESCENT DIRECTION ===================
903 : : ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
904 : : ! ======================================================================
905 : : ! Project all bands from gresid into direc:
906 : : ! The following projection over the subspace orthogonal to occupied bands
907 : : ! is not optional in the RF case, unlike the GS case.
908 : : ! However, the order of operations could be changed, so that
909 : : ! as to make it only applied at the beginning, to H(1) psi(0),
910 : : ! so, THIS IS TO BE REEXAMINED
911 : : ! Note the subtlety:
912 : : ! -For the generalized eigenPb, S|cgq> is used in place of |cgq>,
913 : : ! in order to apply P_c+ projector (see PRB 73, 235101 (2006) [[cite:Audouze2006]], Eq. (71), (72)
914 [ + + ]: 137842034 : do iband = 1, nband
915 [ + + ]: 122128901 : if (bands_treated_now(iband) == 0) cycle
916 [ + + + - : 13440829125 : if (rank_band(iband) == me_band) work = gresid
+ - - + -
- + + +
+ ]
917 : 15738153 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
918 : :
919 : : #ifdef HAVE_OPENMP_OFFLOAD
920 : : !$OMP TARGET UPDATE TO (work) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
921 : : #endif
922 [ + + ]: 15738153 : if(gen_eigenpb)then
923 : : call projbd(gscq,work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
924 : 1877900 : cgq, scprod,0,tim_projbd,useoverlap,me_g0,comm_fft,gpu_option=gs_hamkq%gpu_option)
925 : : else
926 : : call projbd( cgq,work,-1, icgq, 0,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
927 : 13860253 : dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft,gpu_option=gs_hamkq%gpu_option)
928 : : end if
929 : : #ifdef HAVE_OPENMP_OFFLOAD
930 : : !$OMP TARGET UPDATE FROM (work,scprod) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
931 : : #endif
932 : :
933 : 15738153 : call xmpi_sum_master(work, rank_band(iband), mpi_enreg%comm_band, ierr)
934 : :
935 : : ! save this for me_band only
936 : : !TODO: make this a blas call? zaxpy
937 [ + + ]: 47189439 : if (rank_band(iband) == me_band) then
938 [ + - + - : 13440804105 : gresid = work - (mpi_enreg%nproc_band-1)*gresid
- + - - +
+ + + ]
939 : : end if
940 : : end do
941 : :
942 : 15713133 : call cg_zcopy(npw1*nspinor,gresid,direc)
943 : :
944 : : ! ======================================================================
945 : : ! ============== CHECK FOR CONVERGENCE CRITERIA ========================
946 : : ! ======================================================================
947 : :
948 : : ! Compute second-order derivative of the energy using a variational expression
949 : 15713133 : call dotprod_g(prod1,doti,istwf_k,npw1*nspinor,1,cwavef,gresid,me_g0,mpi_enreg%comm_spinorfft)
950 : 15713133 : call dotprod_g(prod2,doti,istwf_k,npw1*nspinor,1,cwavef,gh1c,me_g0,mpi_enreg%comm_spinorfft)
951 : 15713133 : d2te=two*(-prod1+prod2)
952 : : ! write(std_out,'(a,f14.6,a,f14.6)') 'prod1 = ',prod1,' prod2 = ',prod2 ! Keep this debugging feature!
953 : :
954 : : ! Compute <u_m(1)|H(0)-e_m(0)|u_m(1)>
955 : : ! (<u_m(1)|H(0)-e_m(0).S|u_m(1)> if gen. eigenPb),
956 : : ! that should be positive,
957 : : ! except when the eigenvalue eig_mk(0) is higher than
958 : : ! the lowest non-treated eig_mk+q(0). For insulators, this
959 : : ! has no influence, but for metallic occupations,
960 : : ! the conjugate gradient algorithm breaks down. The solution adopted here
961 : : ! is very crude, and rely upon the fact that occupancies of such
962 : : ! levels should be smaller and smaller with increasing nband, so that
963 : : ! a convergence study will give the right result.
964 : : ! The same trick is also used later.
965 : 15713133 : u1h0me0u1=-prod1-prod2
966 : :
967 : : ! Some tolerance is allowed, to account for very small numerical inaccuracies and cancellations.
968 [ + + ]: 15713133 : if(u1h0me0u1<-tol_restart) then ! .and. skipme == 0)then
969 [ - + ]: 178 : if (prtvol==-level.or.prtvol==-19) then
970 : 0 : write(msg,'(a,es22.13e3)') ' cgwf3: u1h0me0u1 = ',u1h0me0u1
971 : 0 : call wrtout(std_out,msg)
972 : : end if
973 [ + + + + ]: 1045306 : cwavef =zero
974 [ + + + + ]: 1045306 : ghc =zero
975 [ + + + + ]: 1045306 : gvnlxc =zero
976 [ - + - - : 178 : if (gen_eigenpb) gsc(:,:)=zero
- - ]
977 : 0 : if (usepaw==1) call pawcprj_set_zero(cwaveprj)
978 : : ! A negative residual will be the signal of this problem ...
979 : 178 : resid=-one
980 [ + + ]: 178 : if (prtvol > 0) call wrtout(std_out,' dfpt_cgwf: problem of minimisation (likely metallic), set resid to -1')
981 : : ! Number of one-way 3D ffts skipped
982 : 178 : nskip=nskip+(nline-iline+1)
983 : :
984 : : skipme = 1
985 : : end if
986 : :
987 [ + + ]: 15712955 : if (skipme == 0) then
988 : : ! Compute residual (squared) norm
989 : 15318206 : call sqnorm_g(resid,istwf_k,npw1*nspinor,gresid,me_g0,comm_fft)
990 [ - + ]: 15318206 : if (prtvol==-level.or.prtvol==-19)then
991 : 0 : write(msg,'(a,a,i3,f14.6,a,a,4es12.4)') ch10,&
992 : 0 : ' dfpt_cgwf: iline,eshift =',iline,eshift,ch10,&
993 : 0 : ' resid,prod1,prod2,d2te=',resid,prod1,prod2,d2te
994 : 0 : call wrtout(std_out,msg)
995 : : end if
996 : : end if
997 : :
998 : : ! If residual sufficiently small stop line minimizations
999 [ + + + + ]: 15713133 : if (resid<tolwfr .and. skipme == 0) then
1000 [ + + ]: 705428 : if(prtvol>=10)then
1001 : : write(msg,'(a,i4,a,i2,a,es12.4)')&
1002 : 7793 : ' dfpt_cgwf: band',u1_band,' converged after ',iline,' line minimizations: resid = ',resid
1003 : 7793 : call wrtout(std_out,msg)
1004 : : end if
1005 : 705428 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
1006 : 705428 : skipme = 1
1007 : : end if
1008 : :
1009 : : ! If user require exiting the job, stop line minimisations
1010 [ - + ]: 15713133 : if (quit==1) then
1011 : 0 : write(msg,'(a,i0)')' dfpt_cgwf: user require exiting => skip update of band ',u1_band
1012 : 0 : call wrtout(std_out,msg)
1013 : 0 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
1014 : 0 : exit ! Exit from the loop on iline
1015 : : end if
1016 : :
1017 : : ! Check that d2te is decreasing on succeeding lines:
1018 [ + + ]: 15713133 : if (iline/=1) then
1019 [ + + + - : 11104920 : if (d2te>d2teold+tol6 .and. u1_band_ > 0 .and. prtvol > 0) then
- + ]
1020 : 0 : write(msg,'(a,i0,a,e16.8,a,e16.9)')'New trial energy at iline ',iline,' = ',d2te,' is higher than former: ',d2teold
1021 : 0 : ABI_WARNING(msg)
1022 : : end if
1023 : : end if
1024 : 15713133 : d2teold=d2te
1025 : :
1026 : : !write(msg, *)"iline, resid, skipme", iline, resid, skipme; call wrtout(std_out, msg)
1027 : :
1028 : : !DEBUG Keep this debugging feature !
1029 : : !call sqnorm_g(dotr,istwf_k,npw1*nspinor,direc,me_g0,comm_fft)
1030 : : !write(std_out,*)' dfpt_cgwf: before precon, direc**2=',dotr
1031 : : !if (gen_eigenpb) then
1032 : : !call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,1,cwaveq,&
1033 : : !& gscq(:,1+npw1*nspinor*(u1_band-1)+igscq:npw1*nspinor*u1_band+igscq),me_g0,mpi_enreg%comm_spinorfft)
1034 : : !else
1035 : : !call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwaveq,me_g0,comm_fft)
1036 : : !end if
1037 : : !write(std_out,*)' dfpt_cgwf: before precon, cwaveq**2=',dotr
1038 : : !ENDDEBUG
1039 : :
1040 : : ! ======================================================================
1041 : : ! ======== PRECONDITION THE STEEPEST DESCENT DIRECTION =================
1042 : : ! ======================================================================
1043 : :
1044 : : ! If wfoptalg>=10, the precondition matrix is kept constant
1045 : : ! during iteration; otherwise it is recomputed
1046 [ + + + + ]: 15713133 : if (wfoptalg<10.or.iline==1) then
1047 : : !print *, 'cwaveq cg_precon ', cwaveq; print *, 'direc cg_precon ', direc
1048 [ - + ]: 14446569 : call cg_precon(cwaveq,zero,istwf_k,kinpw1,npw1,nspinor,me_g0,0,pcon,direc,mpi_enreg%comm_fft)
1049 : : else
1050 [ + + ]: 2588504 : do ispinor=1,nspinor
1051 : 1321940 : igs=(ispinor-1)*npw1
1052 : : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(igs,npw,direc,pcon)
1053 [ + + ]: 185259575 : do ipw=1+igs,npw1+igs
1054 [ + + ]: 549335153 : direc(1:2,ipw)=direc(1:2,ipw)*pcon(ipw-igs)
1055 : : end do
1056 : : end do
1057 : : end if
1058 : :
1059 : : !DEBUG Keep this debugging feature !
1060 : : !call sqnorm_g(dotr,istwf_k,npw1*nspinor,direc,me_g0,comm_fft)
1061 : : !write(std_out,*)' dfpt_cgwf: after precon, direc**2=',dotr
1062 : : !ENDDEBUG
1063 : :
1064 : : ! ======================================================================
1065 : : ! ======= PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION ==============
1066 : : ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
1067 : : ! ======================================================================
1068 : :
1069 : : ! Projecting again out all bands:
1070 : : ! -For the simple eigenPb, gscq is used as dummy argument
1071 [ + + ]: 137842034 : do iband = 1, nband
1072 [ + + ]: 122128901 : if (bands_treated_now(iband) == 0) cycle
1073 : :
1074 [ + + + - : 13440829125 : if (rank_band(iband) == me_band) work = direc
+ - - + -
- + + +
+ ]
1075 : 15738153 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
1076 : :
1077 : : #ifdef HAVE_OPENMP_OFFLOAD
1078 : : !$OMP TARGET UPDATE TO (work) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
1079 : : #endif
1080 : : call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
1081 : 15738153 : gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft,gpu_option=gs_hamkq%gpu_option)
1082 : : #ifdef HAVE_OPENMP_OFFLOAD
1083 : : !$OMP TARGET UPDATE FROM (work,scprod) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
1084 : : #endif
1085 : :
1086 : 15738153 : call xmpi_sum_master(work, rank_band(iband), mpi_enreg%comm_band, ierr)
1087 : :
1088 : : ! save this for me_band only
1089 : : !TODO: make this a blas call? zaxpy
1090 [ + + ]: 47189439 : if (rank_band(iband) == me_band) then
1091 [ + - + - : 13440804105 : direc = work - (mpi_enreg%nproc_band-1)*direc
- + - - +
+ + + ]
1092 : : end if
1093 : : end do
1094 : :
1095 : : !DEBUG Keep this debugging feature !
1096 : : !call sqnorm_g(dotr,istwf_k,npw1*nspinor,direc,me_g0,comm_fft)
1097 : : !write(std_out,*)' dfpt_cgwf: after projbd, direc**2=',dotr
1098 : :
1099 : : ! ======================================================================
1100 : : ! ================= COMPUTE THE CONJUGATE-GRADIENT =====================
1101 : : ! ======================================================================
1102 : :
1103 : : ! get dot of direction vector with residual vector
1104 : 15713133 : call dotprod_g(dotgg,doti,istwf_k,npw1*nspinor,1,direc,gresid,me_g0,mpi_enreg%comm_spinorfft)
1105 : :
1106 [ + + ]: 15713133 : if (iline==1) then
1107 : : ! At first iteration, gamma is set to zero
1108 : 4608213 : gamma=zero
1109 : 4608213 : dotgp=dotgg
1110 : 4608213 : call cg_zcopy(npw1*nspinor,direc,conjgr)
1111 : : else
1112 : : ! At next iterations, h = g + gamma * h
1113 : 11104920 : gamma=dotgg/dotgp
1114 : 11104920 : dotgp=dotgg
1115 [ - + ]: 11104920 : if (prtvol==-level.or.prtvol==-19)then
1116 : 0 : write(msg,'(a,2es16.6)') 'dfpt_cgwf: dotgg,gamma = ',dotgg,gamma
1117 : 0 : call wrtout(std_out,msg)
1118 : : end if
1119 : : !$OMP PARALLEL DO
1120 [ + + ]: 3431347901 : do ipw=1,npw1*nspinor
1121 [ + + ]: 10271833863 : conjgr(1:2,ipw)=direc(1:2,ipw)+gamma*conjgr(1:2,ipw)
1122 : : end do
1123 [ + - ]: 11104920 : if (prtvol==-level.or.prtvol==-19) call wrtout(std_out,'dfpt_cgwf: conjugate direction has been found')
1124 : : end if
1125 : :
1126 : : ! ======================================================================
1127 : : ! ===== COMPUTE CONTRIBUTIONS TO 1ST AND 2ND DERIVATIVES OF ENERGY =====
1128 : : ! ======================================================================
1129 : : ! ...along the search direction
1130 : :
1131 : : ! Compute dedt, Eq.(29) of of PRB55, 10337 (1997) [[cite:Gonze1997]],
1132 : : ! with an additional factor of 2 for the difference between E(2) and the 2DTE
1133 : : dedt = zero
1134 : 15713133 : call dotprod_g(dedt,doti,istwf_k,npw1*nspinor,1,conjgr,gresid,me_g0,mpi_enreg%comm_spinorfft)
1135 : 15713133 : dedt=-two*two*dedt
1136 : :
1137 [ + - - + : 15713133 : if((prtvol==-level.or.prtvol==-19.or.prtvol==-20).and.dedt-tol14>0) then
- - ]
1138 : 0 : call wrtout(std_out, ' DFPT_CGWF WARNING: dedt > 0')
1139 : : end if
1140 [ - + + - : 47139399 : ABI_MALLOC(gvnlx_direc,(2,npw1*nspinor))
+ - + - -
+ - + -
+ ]
1141 [ - + + - : 31426266 : ABI_MALLOC(gh_direc,(2,npw1*nspinor))
- + - + -
+ ]
1142 [ + + ]: 15713133 : if (gen_eigenpb) then
1143 [ - + + - : 3728800 : ABI_MALLOC(sconjgr,(2,npw1*nspinor))
- + - + -
+ ]
1144 : : else
1145 [ - + - + ]: 13848733 : ABI_MALLOC(sconjgr,(0,0))
1146 : : end if
1147 [ + + ]: 15713133 : sij_opt=0;if (gen_eigenpb) sij_opt=1
1148 : 15713133 : cpopt=-1+usepaw
1149 : : call getghc(cpopt,conjgr,conjgrprj,gh_direc,sconjgr,gs_hamkq,gvnlx_direc,&
1150 : 15713133 : eshift,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
1151 : :
1152 : : ! ghc also includes the eigenvalue shift
1153 [ + + ]: 15713133 : if (gen_eigenpb) then
1154 : : !$OMP PARALLEL DO
1155 [ + + ]: 273975009 : do ipw=1,npw1*nspinor
1156 [ + + ]: 818196227 : gh_direc(1:2,ipw)=gh_direc(1:2,ipw)-eshift*sconjgr(1:2,ipw)
1157 : : end do
1158 : : else
1159 : : !$OMP PARALLEL DO
1160 [ + + ]: 4211530737 : do ipw=1,npw1*nspinor
1161 : : ! gh_direc(1:2,ipw)=gh_direc(1:2,ipw)-eshift*conjgr(1:2,ipw)
1162 : 4197682004 : gh_direc(1,ipw)=gh_direc(1,ipw)-eshift*conjgr(1,ipw)+eta_*conjgr(2,ipw)
1163 : 4211530737 : gh_direc(2,ipw)=gh_direc(2,ipw)-eshift*conjgr(2,ipw)-eta_*conjgr(1,ipw)
1164 : : end do
1165 : : end if
1166 : :
1167 : : ! compute d2edt2, Eq.(30) of of PRB55, 10337 (1997) [[cite:Gonze1997]],
1168 : : ! with an additional factor of 2 for the difference
1169 : : ! between E(2) and the 2DTE, and neglect of local fields (SC terms)
1170 : : d2edt2 = zero
1171 : 15713133 : call dotprod_g(d2edt2,doti,istwf_k,npw1*nspinor,1,conjgr,gh_direc,me_g0,mpi_enreg%comm_spinorfft)
1172 : 15713133 : d2edt2=two*two*d2edt2
1173 : :
1174 [ - + ]: 15713133 : if(prtvol==-level.or.prtvol==-19)then
1175 : 0 : write(msg,'(a,3es14.6)') 'dfpt_cgwf: dedt,d2edt2,resid=',dedt,d2edt2,resid; call wrtout(std_out,msg)
1176 : : end if
1177 : :
1178 : : ! ======================================================================
1179 : : ! ======= COMPUTE MIXING FACTOR - CHECK FOR CONVERGENCE ===============
1180 : : ! ======================================================================
1181 : : ! see Eq.(31) of PRB55, 10337 (1997) [[cite:Gonze1997]]
1182 : : !
1183 [ + + + + ]: 15713133 : if(d2edt2<-tol_restart .and. skipme == 0)then
1184 : : ! This may happen when the eigenvalue eig_mk(0) is higher than
1185 : : ! the lowest non-treated eig_mk+q(0). The solution adopted here
1186 : : ! is very crude, and rely upon the fact that occupancies of such
1187 : : ! levels should be smaller and smaller with increasing nband, so that
1188 : : ! a convergence study will give the right result.
1189 : 1311 : theta=zero
1190 : :
1191 [ + + + + ]: 2863515 : cwavef=zero
1192 [ + + + + ]: 2863515 : ghc =zero
1193 [ + + + + ]: 2863515 : gvnlxc =zero
1194 [ - + - - : 1311 : if (gen_eigenpb) gsc=zero
- - ]
1195 : 0 : if (usepaw==1) call pawcprj_set_zero(cwaveprj)
1196 : :
1197 : : ! A negative residual will be the signal of this problem ...
1198 : : !write(msg,'(a,3es16.6)') ' dfpt_cgwf: dedt,d2edt2,resid=',dedt,d2edt2,resid; call wrtout(std_out, msg)
1199 : : !write(msg, "(a, es16.6)")" eig0_k(n_k):", eig0_k(u1_band); call wrtout(std_out, msg)
1200 : : !write(msg, "(a, es16.6)")" eig0_kq(nband):", eig0_kq(nband); call wrtout(std_out, msg)
1201 : : !write(msg, "(a, es16.6)")" prev resid:", resid; call wrtout(std_out, msg)
1202 : :
1203 : 1311 : resid=-two
1204 [ + + + - ]: 1311 : if (prtvol > 0 .and. u1_band_ > 0) then
1205 : 225 : call wrtout(std_out,' dfpt_cgwf: problem of minimisation (likely metallic), set resid to -2')
1206 : : end if
1207 : :
1208 [ + + ]: 15711822 : else if (d2edt2 > 1.d-40) then
1209 : : ! Here, the value of theta that gives the minimum
1210 : 15711500 : theta=-dedt/d2edt2
1211 : : !write(std_out,*)' dfpt_cgwf: dedt,d2edt2=',dedt,d2edt2
1212 : :
1213 : : else
1214 [ + - ]: 322 : if (u1_band_ > 0) call wrtout(std_out, " DFPT_CGWF WARNING: d2edt2 is zero, skipping update")
1215 : 322 : theta=zero
1216 : : end if
1217 : :
1218 : : ! Check that result is above machine precision
1219 [ + + ]: 15713133 : if (one+theta==one) then
1220 [ + + + - ]: 1633 : if (prtvol > 0 .and. u1_band_ > 0) then
1221 : 235 : write(msg, '(a,es16.4)' )' dfpt_cgwf: converged with theta= ',theta
1222 : 235 : call wrtout(std_out,msg)
1223 : : end if
1224 : 1633 : nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
1225 : : skipme = 1
1226 : : end if
1227 : :
1228 : : ! ======================================================================
1229 : : ! ================ GENERATE NEW |wf>, H|wf>, Vnl|Wf ... ================
1230 : : ! ======================================================================
1231 : :
1232 [ + + ]: 15711500 : if (skipme == 0) then
1233 [ + + ]: 43834401 : call cg_zaxpy(npw1*nspinor, [theta, zero], conjgr,cwavef)
1234 : : ! Filter the wavefunctions for large modified kinetic energy (see routine mkkin.f)
1235 [ + + ]: 29742680 : do ispinor=1,nspinor
1236 : 15131213 : ipws=(ispinor-1)*npw1
1237 : : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef,kinpw1,ipws,npw1)
1238 [ + + ]: 4142163449 : do ipw=1+ipws,npw1+ipws
1239 [ - + ]: 4127551982 : if(kinpw1(ipw-ipws)>huge(zero)*1.d-11)then
1240 [ # # ]: 0 : cwavef(1:2,ipw)=zero
1241 : : end if
1242 : : end do
1243 : : end do
1244 : :
1245 [ + + ]: 43834401 : call cg_zaxpy(npw1*nspinor, [theta, zero], gh_direc,ghc)
1246 [ + + ]: 43834401 : call cg_zaxpy(npw1*nspinor, [theta, zero], gvnlx_direc,gvnlxc)
1247 : :
1248 [ + + + + ]: 18100951 : if (gen_eigenpb) call cg_zaxpy(npw1*nspinor, [theta, zero], sconjgr, gsc)
1249 : 1744742 : if (usepaw==1) call pawcprj_axpby(theta,one,conjgrprj,cwaveprj)
1250 : : end if
1251 : :
1252 [ - + ]: 15713133 : ABI_FREE(gh_direc)
1253 [ - + ]: 15713133 : ABI_FREE(gvnlx_direc)
1254 [ - + ]: 15713133 : ABI_FREE(sconjgr)
1255 : :
1256 : : ! ======================================================================
1257 : : ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
1258 : : ! ======================================================================
1259 : :
1260 [ + + ]: 15713133 : if (usetolrde__ /= 0) then
1261 : : ! Check reduction in trial energy deltae, Eq.(28) of PRB55, 10337 (1997) [[cite:Gonze1997]]
1262 : 13287034 : deltae=half*d2edt2*theta**2+theta*dedt
1263 : :
1264 [ + + ]: 13287034 : if (iline==1) then
1265 : 4444071 : deold=deltae
1266 : : ! The extra factor of two should be removed !
1267 [ + + + + ]: 8842963 : else if (abs(deltae)<tolrde*two*abs(deold) .and. iline/=nline) then
1268 [ + + - + ]: 1981264 : if(prtvol>=10.or.prtvol==-level.or.prtvol==-19)then
1269 : : write(msg, '(a,i0,1x,a,1p,e12.4,a,e12.4,a)' ) &
1270 : 82900 : ' dfpt_cgwf: line',iline,' deltae=',deltae,' < tolrde*',deold,' =>skip lines'
1271 : 82900 : call wrtout(std_out,msg)
1272 : : end if
1273 : 1981264 : nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
1274 : 1981264 : skipme = 1
1275 : : end if
1276 : : end if
1277 : :
1278 : : ! if all bands are skippable, we can exit the iline loop for good.
1279 : : ! otherwise, all procs are needed for the projbd and other operations,
1280 : : ! even if the present band will not be updated
1281 [ + + ]: 137842034 : bands_skipped_now = 0
1282 : 15713133 : bands_skipped_now(u1_band) = skipme
1283 [ - + ]: 15713133 : if (u1_band_ < 0) bands_skipped_now(u1_band) = 0 ! Handle negative index
1284 : 15713133 : call xmpi_sum(bands_skipped_now,mpi_enreg%comm_band,ierr)
1285 : :
1286 : : ! bands_skipped_now = bands_skipped_now - bands_treated_now
1287 [ + + + + ]: 137842034 : if (sum(abs(bands_skipped_now - bands_treated_now)) == 0) exit
1288 : :
1289 : 17308209 : nlines_done = nlines_done + 1
1290 : : end do ! iline
1291 : : ABI_NVTX_END_RANGE()
1292 : :
1293 : : !--------------------------------------------------------------------------
1294 : : ! DEBUG
1295 : : !--------------------------------------------------------------------------
1296 : : ! Check that final cwavef (Psi^(1)) satisfies the orthogonality condition
1297 [ + - ]: 4643781 : if (prtvol==-level.or.prtvol==-19) then
1298 [ # # ]: 0 : sij_opt=0 ; usevnl=1 ; optlocal=1 ; optnl=2 ; if (gen_eigenpb) sij_opt=1
1299 : : !ABI_MALLOC(work1,(2,npw1*nspinor))
1300 [ # # # # : 0 : ABI_MALLOC(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
# # # # #
# # # #
# ]
1301 [ # # # # : 0 : ABI_MALLOC(work2,(2,npw1*nspinor*sij_opt))
# # # # #
# # # #
# ]
1302 : 0 : iband_me = 0
1303 [ # # ]: 0 : do iband=1,nband
1304 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
1305 [ # # ]: 0 : if (rank_band(iband) == me_band) then
1306 : 0 : iband_me = iband_me+1
1307 [ # # ]: 0 : if (gen_eigenpb) then
1308 [ # # # # ]: 0 : work(:,:)=gscq(:,1+npw1*nspinor*(iband-1)+igscq:npw1*nspinor*iband+igscq)
1309 : : else
1310 [ # # # # ]: 0 : work(:,:)=cgq(:,1+npw1*nspinor*(iband-1)+icgq:npw1*nspinor*iband+icgq)
1311 : : end if
1312 : : end if
1313 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
1314 : :
1315 : : ! Compute: <Psi^(0)_i,k+q|Psi^(1)_j,k,q>
1316 : 0 : call dotprod_g(prod1,prod2,istwf_k,npw1*nspinor,2,work,cwavef,me_g0,mpi_enreg%comm_spinorfft)
1317 : :
1318 [ # # # # ]: 0 : if (ipert/=natom+10.and.ipert/=natom+11) then
1319 [ # # ]: 0 : if (gen_eigenpb) then
1320 : : ! note: hardcoded ndat=1
1321 : : call getgh1c(berryopt,cwave0,cwaveprj0,work1,gberry,work2,gs_hamkq,gvnlx1_saved,idir,ipert,(/eshift/),&
1322 [ # # ]: 0 : mpi_enreg,1,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
1323 : :
1324 [ # # ]: 0 : if (rank_band(iband) == me_band) then
1325 [ # # # # ]: 0 : work(:,:)=cgq(:,1+npw1*nspinor*(iband_me-1)+icgq:npw1*nspinor*iband_me+icgq)
1326 : : end if
1327 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
1328 : 0 : call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,work,work2,me_g0,mpi_enreg%comm_spinorfft)
1329 : : else
1330 : 0 : dotr=zero; doti=zero
1331 : : end if
1332 : 0 : dotr=prod1+half*dotr
1333 : 0 : doti=prod2+half*doti
1334 [ # # ]: 0 : else if(prtvol==-19) then ! 2nd order case
1335 : 0 : dotr=prod1+half*rf2%amn(1,iband+(u1_band-1)*nband)
1336 : 0 : doti=prod2+half*rf2%amn(2,iband+(u1_band-1)*nband)
1337 : : else
1338 : 0 : write(msg,'(a)') 'CGWF3_WARNING : Use prtvol=-19 to test orthogonality for ipert=natom+10 or +11'
1339 : 0 : call wrtout(std_out,msg)
1340 : : end if
1341 : 0 : dotr=sqrt(dotr**2+doti**2)
1342 [ # # ]: 0 : if(dotr>tol10) then
1343 : : !if (gen_eigenpb) then
1344 : : ! write(msg,'(2a,i3,a,2es22.15)') 'CGWF3_WARNING : <Psi^(1)_i,k,q|S^(0)|Psi^(0)_j,k+q>',&
1345 : : ! '+ 1/2<Psi^(0)_i,k|S^(1)|Psi^(0)_j,k+q>, for j= ',iband,' is ',dotr,doti
1346 : : !else
1347 : 0 : write(msg,'(a,i0,a,es22.15)') 'CGWF3_WARNING : |<Psi^(0)_i,k+q|Psi^(1)_j,k,q>+amn(i,j)/2|, for j= ',iband,' is ',dotr
1348 : : !end if
1349 : 0 : call wrtout(std_out,msg)
1350 : : end if
1351 : : end do
1352 [ # # ]: 0 : ABI_FREE(work1)
1353 [ # # ]: 0 : ABI_FREE(work2)
1354 [ # # # # ]: 0 : if (ipert/=natom+10.and.ipert/=natom+11) then
1355 [ # # ]: 0 : ABI_FREE(gvnlx1_saved)
1356 : : end if
1357 : : end if ! prtvol==-level.or.prtvol==-19
1358 : :
1359 [ - + ]: 4643781 : if (prtvol==-level.or.prtvol==-19)then
1360 : : ! Check that final cwavef Psi^(1) is Pc.Psi^(1)+delta_Psi^(1)
1361 [ # # # # : 0 : ABI_MALLOC(cwwork,(2,npw1*nspinor))
# # # # #
# # # #
# ]
1362 : : ! -Apply Pc to Psi^(1)
1363 [ # # ]: 0 : do iband = 1, nband
1364 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
1365 [ # # # # : 0 : if (rank_band(iband) == me_band) cwwork=cwavef
# # # # #
# # # #
# ]
1366 : :
1367 : 0 : call xmpi_bcast(cwwork,rank_band(iband),mpi_enreg%comm_band,ierr)
1368 : :
1369 [ # # ]: 0 : if(gen_eigenpb)then
1370 : : call projbd(cgq,cwwork,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
1371 : 0 : cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1372 : : else
1373 : : call projbd(cgq,cwwork,-1,icgq,0,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
1374 : 0 : dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1375 : : end if
1376 : :
1377 : 0 : call xmpi_sum(cwwork,mpi_enreg%comm_band,ierr)
1378 : :
1379 : : ! save this for me_band only
1380 : : !TODO: make this a blas call? zaxpy
1381 [ # # ]: 0 : if (rank_band(iband) == me_band) then
1382 [ # # # # : 0 : cwwork = cwwork - (mpi_enreg%nproc_band-1)*cwavef
# # # # #
# # # ]
1383 : : end if
1384 : : end do
1385 : :
1386 : : ! -Add delta_Psi^(1)
1387 [ # # # # : 0 : if (usedcwavef>0) cwwork=cwwork+dcwavef
# # # # #
# # # #
# ]
1388 [ # # # # : 0 : if(ipert==natom+10.or.ipert==natom+11) cwwork=cwwork+rf2%dcwavef(:,1+dc_shift_band:npw1*nspinor+dc_shift_band)
# # # # #
# ]
1389 : : ! -Compare to Psi^(1)
1390 [ # # # # : 0 : cwwork=cwwork-cwavef
# # # # #
# # # ]
1391 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwwork,me_g0,comm_fft)
1392 [ # # ]: 0 : ABI_FREE(cwwork)
1393 [ # # ]: 0 : if(sqrt(dotr)>tol10) then
1394 : : !if (gen_eigenpb) then
1395 : : ! write(msg,'(a,i3,a,es22.15)') &
1396 : : ! 'CGWF3_WARNING : |(Pc.Psi^(1)_i,k,q + delta_Psi^(1)_i,k) - Psi^(1)_i,k,q|^2 (band ',u1_band,')=',dotr
1397 : : !else
1398 : 0 : write(msg,'(a,es22.15)') 'CGWF3_WARNING : |(Pc.Psi^(1)_i,k,q + delta_Psi^(1)_i,k) - Psi^(1)_i,k,q| = ',sqrt(dotr)
1399 : : !end if
1400 : 0 : call wrtout(std_out,msg)
1401 : : end if
1402 : : end if ! prtvol==-level.or.prtvol==-19
1403 : :
1404 [ + - - + ]: 4643781 : if(prtvol==-level.or.prtvol==-19.or.prtvol==-20)then
1405 : : ! Check that final cwavef (Psi^(1)) solves the Sternheimer equation
1406 [ # # # # : 0 : ABI_MALLOC(cwwork,(2,npw1*nspinor))
# # # # #
# # # #
# ]
1407 : : ! -Apply Pc to Psi^(1)
1408 [ # # ]: 0 : do iband = 1, nband
1409 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
1410 : :
1411 [ # # # # : 0 : if (rank_band(iband) == me_band) cwwork=cwavef
# # # # #
# # # #
# ]
1412 : 0 : call xmpi_bcast(cwwork,rank_band(iband),mpi_enreg%comm_band,ierr)
1413 : :
1414 [ # # ]: 0 : if(gen_eigenpb)then
1415 : : call projbd(cgq,cwwork,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
1416 : 0 : gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1417 : : else
1418 : : call projbd(cgq,cwwork,-1,icgq,0,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
1419 : 0 : dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1420 : : end if
1421 : :
1422 : 0 : call xmpi_sum_master(cwwork, rank_band(iband), mpi_enreg%comm_band, ierr)
1423 : :
1424 : : ! save this for me_band only
1425 : : !TODO: make this a blas call? zaxpy
1426 [ # # ]: 0 : if (rank_band(iband) == me_band) then
1427 [ # # # # : 0 : cwwork = cwwork - (mpi_enreg%nproc_band-1)*cwavef
# # # # #
# # # ]
1428 : : end if
1429 : : end do
1430 : :
1431 : : ! - Apply H^(0)-E.S^(0)
1432 [ # # ]: 0 : sij_opt=0;if (gen_eigenpb) sij_opt=1
1433 : 0 : cpopt=-1
1434 [ # # # # : 0 : ABI_MALLOC(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
# # # # #
# # # #
# ]
1435 [ # # # # : 0 : ABI_MALLOC(work2,(2,npw1*nspinor))
# # # # #
# # # #
# ]
1436 : : call getghc(cpopt,cwwork,conjgrprj,work,work1,gs_hamkq,work2,eshift,&
1437 : 0 : mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
1438 [ # # ]: 0 : if (gen_eigenpb) then
1439 [ # # # # : 0 : cwwork=work-eshift*work1
# # # # #
# # # ]
1440 : : else
1441 [ # # # # : 0 : cwwork=work-eshift*cwwork
# # # # #
# # # ]
1442 : : end if
1443 [ # # ]: 0 : ABI_FREE(work1)
1444 [ # # ]: 0 : ABI_FREE(work2)
1445 : :
1446 : : ! The following is not mandatory, as Pc has been already applied to Psi^(1)
1447 : : ! and Pc^* H^(0) Pc = Pc^* H^(0) = H^(0) Pc (same for S^(0)).
1448 : : ! However, in PAW, to apply Pc^* here seems to reduce the numerical error
1449 : : ! -Apply Pc^*
1450 [ # # ]: 0 : do iband = 1, nband
1451 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
1452 [ # # # # : 0 : if (rank_band(iband) == me_band) work=cwwork
# # # # #
# # # # #
# # # # ]
1453 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
1454 : :
1455 [ # # ]: 0 : if(gen_eigenpb)then
1456 : : call projbd(gscq, work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
1457 : 0 : cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1458 : : else
1459 : : call projbd(cgq,work,-1,icgq,0,istwf_k,mcgq,mgscq,nband_me,npw1,nspinor,&
1460 : 0 : dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1461 : : end if
1462 : :
1463 : 0 : call xmpi_sum(work,mpi_enreg%comm_band,ierr)
1464 : :
1465 : : ! save this for me_band only
1466 : : !TODO: make this a blas call? zaxpy
1467 [ # # ]: 0 : if (rank_band(iband) == me_band) then
1468 [ # # # # : 0 : cwwork = cwwork - (mpi_enreg%nproc_band-1)*work
# # ]
1469 : : end if
1470 : : end do
1471 : :
1472 : : ! - Add Pc^*(H^(1)-E.S^(1)).Psi^(0)
1473 [ # # # # : 0 : cwwork=cwwork+gh1c
# # ]
1474 : 0 : call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwwork,me_g0,comm_fft)
1475 [ # # ]: 0 : ABI_FREE(cwwork)
1476 : : write(msg,'(a,i0,a,es22.15,2a,i4)') &
1477 : 0 : '*** CGWF3 Sternheimer equation test for band ',u1_band,'=',sqrt(dotr),ch10,&
1478 : 0 : 'It should go to zero for large nline : nlines_done = ',nlines_done
1479 : 0 : call wrtout(std_out,msg)
1480 : : end if ! prtvol==-level.or.prtvol==-19.or.prtvol==-20
1481 : :
1482 [ + - - + ]: 4643781 : if(prtvol==-level.or.prtvol==-19.or.prtvol==-20)then
1483 : : ! Check that < Psi^(0) | ( H^(0)-eps^(0) S^(0) ) | Psi^(1) > is in agreement with eig^(1)
1484 [ # # # # : 0 : ABI_MALLOC(cwwork,(2,npw1*nspinor))
# # # # #
# # # #
# ]
1485 [ # # # # : 0 : cwwork=cwavef
# # # # #
# ]
1486 : : ! - Apply H^(0)-E.S^(0)
1487 [ # # ]: 0 : sij_opt=0;if (gen_eigenpb) sij_opt=1
1488 : 0 : cpopt=-1
1489 [ # # # # : 0 : ABI_MALLOC(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
# # # # #
# # # #
# ]
1490 [ # # # # : 0 : ABI_MALLOC(work2,(2,npw1*nspinor))
# # # # #
# ]
1491 : : call getghc(cpopt,cwwork,conjgrprj,work,work1,gs_hamkq,work2,eshift,&
1492 : 0 : mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
1493 [ # # ]: 0 : if (gen_eigenpb) then
1494 [ # # # # : 0 : cwwork=work-eshift*work1
# # # # #
# # # ]
1495 : : else
1496 [ # # # # : 0 : cwwork=work-eshift*cwwork
# # # # #
# # # ]
1497 : : end if
1498 [ # # ]: 0 : ABI_FREE(work1)
1499 [ # # ]: 0 : ABI_FREE(work2)
1500 [ # # # # : 0 : cwwork=cwwork+gh1c_n
# # # # #
# # # ]
1501 : 0 : jband=(u1_band-1)*2*nband
1502 : 0 : iband_me = 0
1503 [ # # ]: 0 : do iband=1,nband
1504 [ # # ]: 0 : if (bands_treated_now(iband) == 0) cycle
1505 [ # # ]: 0 : if (rank_band(iband) == me_band) then
1506 : 0 : iband_me = iband_me+1
1507 [ # # # # ]: 0 : work(:,:)=cgq(:,1+npw1*nspinor*(iband_me-1)+icgq:npw1*nspinor*iband_me+icgq)
1508 : : end if
1509 : 0 : call xmpi_bcast(work,rank_band(iband),mpi_enreg%comm_band,ierr)
1510 : :
1511 : 0 : call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,work,cwwork,me_g0,mpi_enreg%comm_spinorfft)
1512 : 0 : dotr = dotr - eig1_k(2*iband-1+jband)
1513 : 0 : doti = doti - eig1_k(2*iband +jband)
1514 : 0 : dotr = sqrt(dotr**2+doti**2)
1515 [ # # ]: 0 : if (dotr > tol8) then
1516 : : write(msg,'(2(a,i0),a,es22.15)') &
1517 : 0 : 'CGWF3_WARNING < Psi^(0) | ( H^(0)-eps^(0) S^(0) ) | Psi^(1) > for i=',iband,' j=',u1_band,&
1518 : 0 : ' : ',sqrt(dotr**2+doti**2)
1519 : 0 : call wrtout(std_out,msg)
1520 : : end if
1521 : : end do
1522 : : !write(std_out,'(a)') '< Psi^(0) | ( H^(0)-eps^(0) S^(0) ) | Psi^(1) > is done.'
1523 [ # # ]: 0 : ABI_FREE(cwwork)
1524 : : end if ! prtvol==-level.or.prtvol==-19.or.prtvol==-20
1525 : : !--------------------------------------------------------------------------
1526 : : ! END DEBUG
1527 : : !--------------------------------------------------------------------------
1528 : :
1529 [ - + ]: 4643781 : ABI_SFREE(gh_direc)
1530 [ - + ]: 4643781 : ABI_SFREE(gvnlx_direc)
1531 : :
1532 [ - + ]: 4643781 : ABI_FREE(conjgr)
1533 [ - + ]: 4643781 : ABI_FREE(cwaveq)
1534 [ - + ]: 4643781 : ABI_FREE(direc)
1535 [ - + ]: 4643781 : ABI_FREE(gresid)
1536 [ + + ]: 4643781 : if (usepaw==1) call pawcprj_free(conjgrprj)
1537 [ + - + + : 5995269 : ABI_FREE(conjgrprj)
- + - + -
+ ]
1538 : :
1539 [ + + ]: 4643781 : if(u1_band>max(1,nband-nbdbuf))then
1540 : : ! A small negative residual will be associated with these
1541 : : ! in the present algorithm all bands need to be in the loops over cgq etc... for the parallelization
1542 : 398682 : resid=-0.1_dp
1543 : : end if
1544 : :
1545 : : ! At the end of the treatment of a set of bands, write the number of one-way 3D ffts skipped
1546 [ + + + + ]: 4643781 : if (xmpi_paral==1 .and. u1_band==nband .and. prtvol>=10) then
1547 : 51182 : write(msg,'(a,i0)')' dfpt_cgwf: number of one-way 3D ffts skipped in cgwf3 until now: ',nskip
1548 : 51182 : call wrtout(std_out,msg)
1549 : : end if
1550 : :
1551 : : #ifdef HAVE_OPENMP_OFFLOAD
1552 : : !$OMP TARGET EXIT DATA MAP(delete:work,scprod) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
1553 : : #endif
1554 [ - + ]: 4643781 : ABI_FREE(work)
1555 [ - + ]: 4643781 : ABI_FREE(gh1c)
1556 [ - + ]: 4643781 : ABI_FREE(pcon)
1557 [ - + ]: 4643781 : ABI_FREE(scprod)
1558 [ - + ]: 4643781 : ABI_FREE(gberry)
1559 : :
1560 : 4643781 : call timab(122,2,tsec)
1561 : :
1562 : : ABI_NVTX_END_RANGE()
1563 : :
1564 : : DBG_EXIT("COLL")
1565 : :
1566 [ - + - + : 13931343 : end subroutine dfpt_cgwf
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - + ]
1567 : : !!***
1568 : :
1569 : : !!****f* ABINIT/full_active_wf1
1570 : : !!
1571 : : !! NAME
1572 : : !! full_active_wf1
1573 : : !!
1574 : : !! FUNCTION
1575 : : !! Restore the full "active space" contribution to the 1st-order wavefunctions.
1576 : : !! The 1st-order WF corrected in this way will no longer be orthogonal to the other occupied states.
1577 : : !! This routine will be only used in a non self-consistent calculation of the
1578 : : !! 1st-order WF for post-processing purposes. Therefore, it does not compute
1579 : : !! the contribution of the 2DTE coming from the change of occupations.
1580 : : !!
1581 : : !! INPUTS
1582 : : !! cgq(2,mcgq)=planewave coefficients of wavefunctions at k+q
1583 : : !! cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
1584 : : !! cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
1585 : : !! cwaveprj(natom,nspinor)= 1st-order wave-function before correction
1586 : : !! projected on NL projectors (PAW)
1587 : : !! cycle_bands(nband)=array of logicals for bands we have on this cpu
1588 : : !! eig1(2*nband**2)=first-order eigenvalues (hartree)
1589 : : !! fermie1=derivative of fermi energy wrt (strain) perturbation
1590 : : !! eig0nk=energy of the band at k being corrected
1591 : : !! eig0_kq(nband)=energies of the bands at k+q
1592 : : !! elph2_imagden=imaginary parameter to broaden the energy denominators
1593 : : !! iband=index of current band
1594 : : !! ibgq=shift to be applied on the location of data in the array cprjq
1595 : : !! icgq=shift to be applied on the location of data in the array cgq
1596 : : !! mcgq=second dimension of the cgq array
1597 : : !! mcprjq=second dimension of the cprjq array
1598 : : !! mpi_enreg=information about MPI parallelization
1599 : : !! natom=number of atoms in cell
1600 : : !! nband=number of bands
1601 : : !! npw1=number of plane waves at this k+q point
1602 : : !! nspinor=number of spinorial components of the wavefunctions
1603 : : !! timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
1604 : : !! usepaw=flag for PAW
1605 : : !!
1606 : : !! OUTPUT
1607 : : !! cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
1608 : : !! cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
1609 : : !! projected on NL projectors (PAW)
1610 : : !!
1611 : : !! SOURCE
1612 : :
1613 : 88145 : subroutine full_active_wf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,eig1,&
1614 : 88145 : fermie1,eig0nk,eig0_kq,elph2_imagden,&
1615 : : iband,ibgq,icgq,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,&
1616 : : nspinor,timcount,usepaw)
1617 : :
1618 : : !Arguments ------------------------------------
1619 : : !scalars
1620 : : integer,intent(in) :: iband,ibgq,icgq,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw
1621 : : real(dp),intent(in) :: fermie1, eig0nk, elph2_imagden
1622 : : type(MPI_type),intent(in) :: mpi_enreg
1623 : : !arrays
1624 : : logical,intent(in) :: cycle_bands(nband)
1625 : : real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor)
1626 : : real(dp),intent(in) :: eig0_kq(nband), eig1(2*nband**2)
1627 : : real(dp),intent(out) :: cwave1(2,npw1*nspinor)
1628 : : type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw)
1629 : : type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw)
1630 : :
1631 : : !Local variables-------------------------------
1632 : : !scalars
1633 : : integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii,ibandkq_me, ierr
1634 : : real(dp) :: facti,factr,eta,delta_E,inv_delta_E,gkkr
1635 : : !arrays
1636 : : real(dp) :: tsec(2)
1637 : : ! *********************************************************************
1638 : :
1639 : : DBG_ENTER("COLL")
1640 : :
1641 : 88145 : call timab(214+timcount,1,tsec)
1642 : :
1643 : : ! At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf).
1644 : : ! Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 .
1645 : :
1646 : : ! New logic 11/11/2019: accumulate correction in cwave1 and cwaveprj1,
1647 : : ! then add it to cwavef at the end with a modified blas call
1648 [ + + + + ]: 323369354 : cwave1 = zero
1649 : :
1650 [ - + ]: 88145 : if (usepaw==1) call pawcprj_set_zero(cwaveprj1)
1651 : :
1652 : 88145 : eta = elph2_imagden
1653 : :
1654 : : ! Loop over WF at k+q subspace
1655 : 88145 : ibandkq_me = 0
1656 [ + + ]: 1030081 : do ibandkq=1,nband
1657 : :
1658 : : !TODO MJV: here we have an issue - the cgq are no longer present for all bands!
1659 : : ! we only have diagonal terms for iband iband1 and ibandq in same set of bands
1660 : : ! 1) filter with distrb
1661 [ + - ]: 941936 : if (cycle_bands(ibandkq)) cycle
1662 : 941936 : ibandkq_me = ibandkq_me + 1
1663 : :
1664 : : ! 2) get contributions for correction factors of cgq from bands present on this cpu
1665 : 941936 : delta_E = eig0nk - eig0_kq(ibandkq)
1666 : 941936 : inv_delta_E = delta_E / ( delta_E ** 2 + eta ** 2)
1667 : :
1668 : 941936 : index_eig1=2*ibandkq-1+(iband-1)*2*nband
1669 : 941936 : index_cgq=npw1*nspinor*(ibandkq_me-1)+icgq
1670 : :
1671 [ + + ]: 941936 : if(ibandkq==iband) then
1672 : 88145 : gkkr = eig1(index_eig1) - fermie1
1673 : : else
1674 : 853791 : gkkr = eig1(index_eig1)
1675 : : end if
1676 : 941936 : factr = inv_delta_E * gkkr
1677 : 941936 : facti = inv_delta_E * eig1(index_eig1+1)
1678 : : !print *, "ibandkq, iband, factr, facti:", ibandkq, iband, factr, facti
1679 : :
1680 : : ! Apply correction to 1st-order WF
1681 : : !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwave1,facti,factr,index_cgq,npw1,nspinor)
1682 [ + + ]: 1137223944 : do ii=1,npw1*nspinor
1683 : 1136282008 : cwave1(1,ii)=cwave1(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq))
1684 : 1137223944 : cwave1(2,ii)=cwave1(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq))
1685 : : end do
1686 : :
1687 : : ! In the PAW case, also apply correction to projected WF
1688 [ + - ]: 1030081 : if (usepaw==1) then
1689 : 0 : index_cprjq=nspinor*(ibandkq_me-1)+ibgq
1690 [ # # ]: 0 : call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1)
1691 : : end if
1692 : :
1693 : : end do ! Loop over k+q subspace
1694 : :
1695 : : ! 3) reduce over bands to get all contributions to correction. need MPI reduce over band communicator only
1696 : 88145 : call xmpi_sum(cwave1,mpi_enreg%comm_band,ierr)
1697 [ - + ]: 88145 : if (usepaw == 1) call pawcprj_mpi_sum(cwaveprj1, mpi_enreg%comm_band, ierr)
1698 : :
1699 : : ! 4) add correction to the cwave1
1700 : : ! Now add on input WF into output WF
1701 : 88145 : call cg_zaxpy(npw1*nspinor,(/one,zero/),cwavef,cwave1)
1702 : :
1703 : : ! Idem for cprj
1704 [ - + ]: 88145 : if (usepaw==1) call pawcprj_zaxpby((/one,zero/),(/one,zero/),cwaveprj,cwaveprj1)
1705 : :
1706 : 88145 : call timab(214+timcount,2,tsec)
1707 : :
1708 : : DBG_EXIT("COLL")
1709 : :
1710 : 88145 : end subroutine full_active_wf1
1711 : : !!***
1712 : :
1713 : : !!****f* m_dfpt_cgwf/stern_init
1714 : : !! NAME
1715 : : !! stern_init
1716 : : !!
1717 : : !! FUNCTION
1718 : : !! Initialize the object.
1719 : : !!
1720 : : !! INPUT
1721 : : !! dtset
1722 : : !! npw_k
1723 : : !! npw_kq
1724 : : !! nspinor
1725 : : !! nband
1726 : : !! nband_me
1727 : : !! fermie1_idir_ipert
1728 : : !! use_cache
1729 : : !! work_ngfft
1730 : : !! mpi_enreg
1731 : : !! comm_band
1732 : : !!
1733 : : !! OUTPUT
1734 : : !!
1735 : : !! SOURCE
1736 : :
1737 [ + - ]: 2040 : subroutine stern_init(stern, dtset, npw_k, npw_kq, nspinor, nband, nband_me, fermie1_idir_ipert, &
1738 : : use_cache, work_ngfft, mpi_enreg, comm_band)
1739 : :
1740 : : !Arguments ------------------------------------
1741 : : class(stern_t),intent(out) :: stern
1742 : : type(dataset_type),target,intent(in) :: dtset
1743 : : integer,intent(in) :: npw_k, npw_kq, nspinor, nband, nband_me, comm_band
1744 : : real(dp),intent(in) :: fermie1_idir_ipert(3, dtset%natom)
1745 : : integer,intent(in) :: work_ngfft(18)
1746 : : logical,intent(in) :: use_cache
1747 : : type(mpi_type),intent(in) :: mpi_enreg
1748 : :
1749 : : !Local variables ------------------------------
1750 : : integer :: natom, usepaw
1751 : : ! *************************************************************************
1752 : :
1753 : 2040 : natom = dtset%natom; usepaw = dtset%usepaw
1754 : :
1755 : 2040 : stern%npw_k = npw_k; stern%npw_kq = npw_kq; stern%nspinor = nspinor; stern%nband = nband; stern%dtset => dtset
1756 : 2040 : stern%usepaw = usepaw; stern%nband_me = nband_me; stern%usedcwavef = 0
1757 [ + + ]: 38760 : stern%use_cache = use_cache; stern%work_ngfft = work_ngfft
1758 [ + + ]: 2040 : if (stern%use_cache) then
1759 [ - + + - : 2895 : ABI_MALLOC(stern%work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1760 : : end if
1761 : :
1762 [ - + + - : 6120 : ABI_MALLOC(stern%fermie1_idir_ipert, (3, natom))
+ - + - -
+ - + -
+ ]
1763 [ + - - + : 20400 : stern%fermie1_idir_ipert = fermie1_idir_ipert
- - + + +
+ ]
1764 : :
1765 : 2040 : call copy_mpi_enreg(mpi_enreg, stern%mpi_enreg)
1766 : 2040 : stern%mpi_enreg%comm_band = comm_band
1767 : 2040 : stern%mpi_enreg%me_band = xmpi_comm_rank(comm_band)
1768 : 2040 : stern%mpi_enreg%nproc_band = xmpi_comm_size(comm_band)
1769 : 2040 : stern%has_band_para = stern%mpi_enreg%nproc_band /= 1
1770 : :
1771 [ - + + - : 1042166 : ABI_CALLOC(stern%eig1_k, (2, nband, nband))
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
1772 [ + - - - : 6120 : ABI_MALLOC(stern%dcwavef, (2, npw_kq*nspinor*stern%usedcwavef))
+ - - + -
+ - + -
+ ]
1773 [ - + + - : 6120 : ABI_MALLOC(stern%gh1c_n, (2, npw_kq*nspinor))
+ - + - -
+ - + -
+ ]
1774 [ - + + - : 4080 : ABI_MALLOC(stern%ghc, (2, npw_kq*nspinor))
- + - + -
+ ]
1775 [ - + + - : 4080 : ABI_MALLOC(stern%gsc, (2, npw_kq*nspinor))
- + - + -
+ ]
1776 [ - + + - : 4080 : ABI_MALLOC(stern%gvnlxc, (2, npw_kq*nspinor))
- + - + -
+ ]
1777 [ - + + - : 4080 : ABI_MALLOC(stern%gvnlx1, (2, npw_kq*nspinor))
- + - + -
+ ]
1778 : :
1779 : : ! TODO: to distribute cgq and kets memory, use mband_mem per core in band comm, but coordinate everyone with
1780 : : ! the following array (as opposed to the distribution of cg1 which is done in the normal dfpt calls
1781 [ - + + - : 6120 : ABI_MALLOC(stern%bands_treated_now, (nband))
+ - + - -
+ - + -
+ ]
1782 [ - + + - : 4080 : ABI_MALLOC (stern%rank_band, (nband))
- + - + -
+ ]
1783 [ + + ]: 26321 : stern%rank_band = 0
1784 : :
1785 [ + + ]: 2040 : stern%nline_in = min(100, npw_kq); if (dtset%nline > stern%nline_in) stern%nline_in = min(dtset%nline, npw_kq)
1786 [ - + ]: 2040 : if (dtset%tolwfr <= tiny(one)) then
1787 : 0 : ABI_ERROR("tolwfr must be specified in input when solving the Sterheimer equation non-self-consistently")
1788 : : end if
1789 : :
1790 [ - + + - : 8160 : ABI_MALLOC(stern%cgq, (2, npw_kq * nspinor, stern%nband_me))
- + + - +
- + - - +
- + - + ]
1791 [ - + + - : 8160 : ABI_MALLOC(stern%gscq, (2, npw_kq * nspinor, nband_me*usepaw))
+ - - - +
- - + - +
- + - + ]
1792 : :
1793 : 2040 : stern%mcgq = stern%npw_kq * stern%nspinor * stern%nband_me
1794 : 2040 : stern%mgscq = stern%npw_kq * stern%nspinor * stern%nband_me * usepaw
1795 : :
1796 : : ! PAW wave functions at k+q projected with non-local projectors
1797 : 2040 : stern%mcprjq = 0
1798 [ - + + - : 4080 : ABI_MALLOC(stern%cprjq, (natom, stern%mcprjq))
- + - + -
+ ]
1799 [ - + + - : 8160 : ABI_MALLOC(stern%cwaveprj1, (natom, nspinor*usepaw))
+ - - - +
- - + - +
- + - + -
+ - - ]
1800 : :
1801 : 2040 : end subroutine stern_init
1802 : : !!***
1803 : :
1804 : : !!****f* m_dfpt_cgwf/stern_solve
1805 : : !! NAME
1806 : : !! stern_solve
1807 : : !!
1808 : : !! FUNCTION
1809 : : !! Solves the NSCF Sternheimer equation. Simplified wrapper around dfpt_cgwf.
1810 : : !!
1811 : : !! INPUT
1812 : : !!
1813 : : !! OUTPUT
1814 : : !!
1815 : : !! SOURCE
1816 : :
1817 : 104430 : subroutine stern_solve(stern, u1_band, band_me, idir, ipert, qpt, gs_hamkq, rf_hamkq, eig0_k, eig0_kq, &
1818 : 104430 : ug0_nk, cprj0_nk, & ! in
1819 : 104430 : ug1_nkq, cprj1_nkq, err_msg, ierr, & ! out
1820 : : full_cg1, full_ur1, init_mode) ! optional
1821 : :
1822 : : !Arguments ------------------------------------
1823 : : class(stern_t),target,intent(inout) :: stern
1824 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
1825 : : type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
1826 : : integer,intent(in) :: u1_band, band_me, idir, ipert
1827 : : character(len=*),optional,intent(in) :: init_mode
1828 : : !arrays
1829 : : real(dp),intent(in) :: qpt(3), eig0_k(stern%nband), eig0_kq(stern%nband)
1830 : : real(dp),intent(inout) :: ug0_nk(2, stern%npw_k*stern%nspinor), ug1_nkq(2, stern%npw_kq*stern%nspinor)
1831 : : type(pawcprj_type),intent(inout) :: cprj0_nk(gs_hamkq%natom, stern%nspinor*gs_hamkq%usecprj)
1832 : : type(pawcprj_type),intent(inout) :: cprj1_nkq(gs_hamkq%natom, stern%nspinor)
1833 : : integer,intent(out) :: ierr
1834 : : character(len=*),intent(out) :: err_msg
1835 : : real(dp),optional,target,intent(out) :: full_cg1(2, stern%npw_kq*stern%nspinor)
1836 : : complex(gwp),optional,intent(out) :: full_ur1(gs_hamkq%nfft*stern%nspinor)
1837 : :
1838 : : !Local variables ------------------------------
1839 : : !scalars
1840 : : integer,parameter :: berryopt0 = 0, igscq0 = 0, icgq0 = 0, ibgq0 = 0, quit0 = 0, istwfk1 = 1, ndat1 = 1, timcount0 = 0
1841 : : integer :: opt_gvnlx1, grad_berry_size_mpw1, iband, gpu_option
1842 : : real(dp) :: out_resid, fermie1, eig0nk !, dotr
1843 : : logical :: map_cgq, map_vlocal
1844 : : character(len=500) :: init_mode__
1845 : 104430 : type(rf2_t) :: rf2
1846 : : !arrays
1847 : 104430 : real(dp),allocatable :: grad_berry(:,:)
1848 : 104430 : real(dp), contiguous, pointer :: cgq_ptr(:,:,:), vlocal_ptr(:,:,:,:) !, work_ptr(:,:,:,:), gscq_ptr(:,:,:)
1849 : : complex(gwp),allocatable :: cwork_sp(:)
1850 : 104430 : logical :: cycle_bands(stern%nband)
1851 : : #ifdef HAVE_GW_DPC
1852 : 104430 : complex(dp),pointer :: full_ug1_dp_ptr(:)
1853 : : #endif
1854 : : ! *************************************************************************
1855 : :
1856 : : !if (any(abs(gs_hamkq%kpt_k + qpt - gs_hamkq%kpt_kp) > tol12)) then
1857 : : ! ABI_ERROR(sjoin(ktoa(gs_hamkq%kpt_k + qpt), " != ", ktoa(gs_hamkq%kpt_kp)))
1858 : : !end if
1859 : :
1860 : : ! TODO: grad_berry is problematic because in dfpt_cgwf, the array is declared with
1861 : : !
1862 : : ! real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband)
1863 : : !
1864 : : ! and
1865 : : !
1866 : : ! npw1_k = number of plane waves at this k+q point
1867 : : !
1868 : : ! So in principle we should allocate lot of memory to avoid bound checking error!
1869 : : ! For the time being use mpw1 = 0 because mpw1 is not used in this call to dfpt_cgwf
1870 : : ! still it's clear that the treatment of this array must be completely refactored in the DFPT code.
1871 : : !
1872 : 104430 : opt_gvnlx1 = 0 ! gvnlx1 is output
1873 : 104430 : stern%nlines_done = 0
1874 : 104430 : grad_berry_size_mpw1 = 0
1875 : 0 : ABI_MALLOC(grad_berry, (2, stern%nspinor*(berryopt0/4)))
1876 : :
1877 : : !if (psps%usepaw==1) mcprjq = stern%nspinor*mband_mem*mkqmem*nsppol*usecprj
1878 : :
1879 : 104430 : gpu_option = stern%dtset%gpu_option
1880 [ + + - + ]: 104430 : init_mode__ = "None"; if (present(init_mode)) init_mode__ = init_mode
1881 : :
1882 : 104430 : select case (init_mode__)
1883 : : case ("None")
1884 : : ! Init entry in cg1s_kq, either from cache or with zeros.
1885 [ + + ]: 104430 : if (stern%use_cache) then
1886 [ + + + + ]: 38664072 : ug1_nkq = zero
1887 : : !u1c_ib_k = stern%u1c%find_band(band_ks)
1888 : : !if (u1c_ib_k /= -1) then
1889 : : ! call cgtk_change_gsphere(stern%nspinor, &
1890 : : ! stern%u1c%prev_npw_kq, istwfk1, stern%u1c%prev_kg_kq, stern%u1c%prev_cg1s_kq(1,1,ipc,u1c_ib_k), &
1891 : : ! stern%npw_kq, istwfk1, kg_kq, cg1s_kq(1,1,ipc,ib_k), stern%work_ngfft, stern%work)
1892 : : !else
1893 : : ! cg1s_kq(:,:,ipc,ib_k) = zero
1894 : : !end if
1895 : : else
1896 [ + + + + ]: 332154372 : ug1_nkq = zero
1897 : : end if
1898 : :
1899 : : case ("input")
1900 : 0 : continue
1901 : :
1902 : : !case ("random")
1903 : : !call cg_randomize(istwfk1, stern%npw_kq, stern%nspinor, ndat1, stern%mpi_enreg%me_g0, ug1_nkq)
1904 : : !call projbd(cgq, work, -1, icgq, igscq, istwfk1, mcgq, mgscq, nband_me, npw1, nspinor, &
1905 : : ! gscq, scprod, 0, tim_projbd, useoverlap, me_g0, comm_fft, gpu_option=gs_hamkq%gpu_option)
1906 : :
1907 : : case default
1908 [ + - - ]: 104430 : ABI_ERROR(sjoin("Invalid init_mode:", init_mode__))
1909 : : end select
1910 : :
1911 : 104430 : cgq_ptr => stern%cgq
1912 : 104430 : vlocal_ptr => gs_hamkq%vlocal
1913 : :
1914 [ - + ]: 104430 : if (gpu_option == ABI_GPU_OPENMP) then
1915 : : ! Upload cgq array to GPU
1916 : 0 : map_cgq = .not. (xomp_target_is_present(c_loc(cgq_ptr)))
1917 : 0 : map_vlocal = .not. (xomp_target_is_present(c_loc(vlocal_ptr)))
1918 : : #ifdef HAVE_OPENMP_OFFLOAD
1919 : : !$OMP TARGET ENTER DATA MAP(to:cgq_ptr) IF (map_cgq)
1920 : : !$OMP TARGET ENTER DATA MAP(to:vlocal_ptr) IF (map_vlocal)
1921 : : #endif
1922 : : end if
1923 : :
1924 : : !print *, "before dfpt_cgwf
1925 : : call dfpt_cgwf(u1_band, band_me, stern%rank_band, stern%bands_treated_now, berryopt0, &
1926 : : stern%cgq, ug1_nkq, ug0_nk, & ! Important stuff
1927 : : cprj1_nkq, cprj0_nk, rf2, stern%dcwavef, &
1928 : : eig0_k, eig0_kq, stern%eig1_k, &
1929 : : stern%ghc, stern%gh1c_n, grad_berry, stern%gsc, stern%gscq, &
1930 : : gs_hamkq, stern%gvnlxc, stern%gvnlx1, icgq0, idir, ipert, igscq0, &
1931 : : stern%mcgq, stern%mgscq, stern%mpi_enreg, grad_berry_size_mpw1, stern%dtset%natom, stern%nband, stern%nband_me, &
1932 : : stern%dtset%nbdbuf, stern%nline_in, stern%npw_k, stern%npw_kq, stern%nspinor, &
1933 : : opt_gvnlx1, stern%dtset%prtvol, quit0, out_resid, rf_hamkq, stern%dtset%dfpt_sciss, -one, stern%dtset%tolwfr, &
1934 : 104430 : stern%usedcwavef, stern%dtset%wfoptalg, stern%nlines_done, usetolrde=0)
1935 : : !print *, "after dfpt_cgwf
1936 : :
1937 [ - + ]: 104430 : ABI_FREE(grad_berry)
1938 : :
1939 : : if (gpu_option == ABI_GPU_OPENMP) then
1940 : : if (map_vlocal) then
1941 : : #ifdef HAVE_OPENMP_OFFLOAD
1942 : : !$OMP TARGET EXIT DATA MAP(delete:vlocal_ptr)
1943 : : #endif
1944 : : end if
1945 : : end if
1946 : :
1947 : : if (stern%use_cache) then
1948 : : ! Store |Psi_1> to init Sternheimer solver for the next q-point.
1949 : : ABI_UNUSED(qpt(1))
1950 : : !call stern%u1c%store(qpt, stern%npw_kq, stern%nspinor, natom3, bstart_ks, nbcalc_ks, kg_kq, cg1s_kq)
1951 : : end if
1952 : :
1953 : : ! Handle possible convergence error.
1954 [ + - ]: 104430 : err_msg = ""; ierr = 0
1955 [ + - ]: 104430 : if (u1_band > 0) then
1956 [ - + ]: 104430 : if (out_resid > stern%dtset%tolwfr) then
1957 : : write(err_msg, "(2(a,i0),a, 2(a,es13.5), 2a,i0,a)") &
1958 : 0 : " Sternheimer didn't convergence for band: ", u1_band, "/", stern%nband, ch10, &
1959 : 0 : " resid:", out_resid, " >= tolwfr: ", stern%dtset%tolwfr, ch10, &
1960 : 0 : " after nline: ", stern%nlines_done, " iterations. Increase nline and/or tolwfr."
1961 : 0 : ierr = 1
1962 : :
1963 [ + + ]: 104430 : else if (out_resid < zero) then
1964 : : ! This may happen when the eigenvalue eig_mk(0) is higher than the lowest non-treated eig_mk+q(0). See dfpt_cgwf.
1965 : : write(err_msg, "(2(a,i0),a, (a,es13.5), a,i0)") &
1966 : 364 : " Sternheimer solver didn't convergence for band: ", u1_band, "/", stern%nband, ch10, &
1967 : 728 : " resid:", out_resid, ", after nline: ", stern%nlines_done
1968 : : ! write(std_out, *)iband, eig0_k(iband), eig0_kq(iband), eig0_k(iband) - eig0_kq(iband)
1969 : 364 : ierr = -2
1970 : : end if
1971 : :
1972 : : !if (my_rank == master .and. (enough_stern <= 5 .or. stern%dtset%prtvol > 10)) then
1973 : : ! write(msg, "(2(a,es13.5),a,i0)") &
1974 : : ! " Sternheimer converged with resid: ", out_resid, " <= tolwfr: ", dtset%tolwfr, &
1975 : : ! " after nlines_done: ", nlines_done
1976 : : ! enough_stern = enough_stern + 1
1977 : : !end if
1978 : : end if ! u1_band > 0
1979 : :
1980 [ + + ]: 105522 : if (ierr /= 0) return
1981 : :
1982 [ + + ]: 104066 : if (present(full_ur1)) then
1983 [ - + ]: 87572 : ABI_CHECK(present(full_cg1), "full_ur1 requires full_cg1")
1984 : : end if
1985 : :
1986 : : ! At this stage, the 1st order function ug1_nkq is orthogonal to cgq (unlike
1987 : : ! when it is input to dfpt_cgwf). Here, restore the "active space" content
1988 : : ! of the first-order wavefunction, to give cwave1.
1989 : : ! PAW: note that dcwavef (1st-order change of WF due to overlap change)
1990 : : ! remains in the subspace orthogonal to cgq
1991 : :
1992 : : ! A similar section of code is present in m_dfpt_vtowfk when if (dtset%prtfull1wf > 0).
1993 : : ! Also, note that in metals one should include fermie1
1994 : : ! See need_fermie1 in m_dfpt_scfcv
1995 : :
1996 [ + + ]: 104066 : if (present(full_cg1)) then
1997 : : ! =====================================
1998 : : ! Compute full first order wavefunction
1999 : : ! =====================================
2000 : :
2001 : : ! NOTE: Assuming all bands at k+q are on this cpu.
2002 [ + + ]: 1015756 : cycle_bands(:) = .False.
2003 [ - + ]: 87572 : ABI_CHECK_IGEQ(u1_band, 1, "u1_band")
2004 : 87572 : eig0nk = eig0_k(u1_band)
2005 [ + + + + ]: 350288 : fermie1 = zero; if (sum(qpt**2) < tol14) fermie1 = stern%fermie1_idir_ipert(idir, ipert)
2006 : :
2007 : 87572 : iband = u1_band
2008 : : call full_active_wf1(stern%cgq, stern%cprjq, ug1_nkq, full_cg1, cprj1_nkq, stern%cwaveprj1, cycle_bands, stern%eig1_k, fermie1, &
2009 : : eig0nk, eig0_kq, stern%dtset%elph2_imagden, iband, ibgq0, icgq0, stern%mcgq, stern%mcprjq, stern%mpi_enreg, &
2010 : 87572 : stern%dtset%natom, stern%nband, stern%npw_kq, stern%nspinor, timcount0, gs_hamkq%usepaw)
2011 : :
2012 [ + - ]: 87572 : if (present(full_ur1)) then
2013 : : ! Note the use use of _kp pointers in gs_hamkq as full_ug1 is given on the k+q g-sphere.
2014 : : #ifdef HAVE_GW_DPC
2015 : : ! we are using double precision -> cast full_cg1 to dp complex pointer.
2016 [ + + ]: 175144 : call c_f_pointer(c_loc(full_cg1), full_ug1_dp_ptr, [stern%npw_kq * stern%nspinor])
2017 : : call fft_ug(stern%npw_kq, gs_hamkq%nfft, stern%nspinor, ndat1, gs_hamkq%mgfft, gs_hamkq%ngfft, &
2018 [ + - - - : 87572 : istwfk1, gs_hamkq%kg_kp, gs_hamkq%gbound_kp, full_ug1_dp_ptr, full_ur1)
- - - - +
- - - - -
- - - + -
+ - + ]
2019 : : #else
2020 : : ! Transfer cgs_kbz from dp to sp and perform FFT in single precision.
2021 : : ABI_MALLOC(cwork_sp, (stern%npw_kq*stern%nspinor))
2022 : : cwork_sp(:) = full_cg1(1,:) + j_sp * full_cg1(2,:)
2023 : : call fft_ug(stern%npw_kq, gs_hamkq%nfft, stern%nspinor, ndat1, gs_hamkq%mgfft, gs_hamkq%ngfft, &
2024 : : istwfk1, gs_hamkq%kg_kp, gs_hamkq%gbound_kp, cwork_sp, full_ur1)
2025 : : ABI_FREE(cwork_sp)
2026 : : #endif
2027 : : end if
2028 : : end if
2029 : :
2030 [ - + + - : 209588 : end subroutine stern_solve
- + - + -
+ + - +
- ]
2031 : : !!***
2032 : :
2033 : : !!****f* m_dfpt_cgwf/stern_free
2034 : : !! NAME
2035 : : !! stern_free
2036 : : !!
2037 : : !! FUNCTION
2038 : : !! Free dynamic memory
2039 : : !!
2040 : : !! SOURCE
2041 : :
2042 : 5313 : subroutine stern_free(stern)
2043 : :
2044 : : !Arguments ------------------------------------
2045 : : class(stern_t),target,intent(inout) :: stern
2046 : :
2047 : : !Local variables ------------------------------
2048 : 5313 : real(dp), contiguous, pointer :: cgq_ptr(:,:,:) !, work_ptr(:,:,:,:), gscq_ptr(:,:,:)
2049 : : !************************************************************************
2050 : :
2051 : : ! integer
2052 : 2040 : ABI_SFREE(stern%bands_treated_now)
2053 [ + + ]: 5313 : ABI_SFREE(stern%rank_band)
2054 : :
2055 : : ! real
2056 [ + + ]: 5313 : ABI_SFREE(stern%fermie1_idir_ipert)
2057 [ + + ]: 5313 : ABI_SFREE(stern%eig1_k)
2058 [ + + ]: 5313 : ABI_SFREE(stern%dcwavef)
2059 [ + + ]: 5313 : ABI_SFREE(stern%gh1c_n)
2060 [ + + ]: 5313 : ABI_SFREE(stern%ghc)
2061 [ + + ]: 5313 : ABI_SFREE(stern%gsc)
2062 [ + + ]: 5313 : ABI_SFREE(stern%gvnlxc)
2063 : :
2064 [ + + ]: 5313 : ABI_SFREE(stern%gscq)
2065 [ + + ]: 5313 : ABI_SFREE(stern%gvnlx1)
2066 [ + + ]: 5313 : ABI_SFREE(stern%work)
2067 : :
2068 : : !call stern%u1c%free()
2069 : 5313 : call destroy_mpi_enreg(stern%mpi_enreg)
2070 : : !call xmpi_comm_free(stern%mpi_enreg%comm_band)
2071 : :
2072 : : !if (stern%usepaw == 1) then
2073 : : ! call pawcprj_free(stern%cprjq)
2074 : : ! call pawcprj_free(stern%cwaveprj1)
2075 : : !end if
2076 [ + + - + : 5313 : ABI_SFREE(stern%cprjq)
- - - - -
+ ]
2077 [ + + - + : 5313 : ABI_SFREE(stern%cwaveprj1)
- - - - -
+ ]
2078 : :
2079 : 5313 : cgq_ptr => stern%cgq
2080 : : #ifdef HAVE_OPENMP_OFFLOAD
2081 : : ! Free array on the GPU
2082 : : if (xomp_target_is_present(c_loc(cgq_ptr))) then
2083 : : !$OMP TARGET EXIT DATA MAP(delete:cgq_ptr)
2084 : : end if
2085 : : #endif
2086 [ + + ]: 5313 : ABI_SFREE(stern%cgq)
2087 : :
2088 [ + + ]: 5313 : end subroutine stern_free
2089 : : !!***
2090 : :
2091 [ - + - - : 4080 : end module m_dfpt_cgwf
+ + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ + - - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - + -
+ - - - -
- - - + ]
2092 : : !!***
|
