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 4629309 : subroutine dfpt_cgwf(u1_band_,band_me,rank_band,bands_treated_now,berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,&
249 4629309 : & rf2,dcwavef,&
250 4629309 : & eig0_k,eig0_kq,eig1_k,ghc,gh1c_n,grad_berry,gsc,gscq,&
251 4629309 : & 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 9258618 : integer :: skipme, bands_skipped_now(nband)
302 : character(len=500) :: msg
303 : !arrays
304 : real(dp) :: dummy(0,0),tsec(2)
305 9258618 : real(dp) :: eig1_k_loc(2,nband,nband)
306 9258618 : real(dp),allocatable :: conjgr(:,:),cwaveq(:,:),cwwork(:,:),direc(:,:)
307 4629309 : real(dp),allocatable :: gberry(:,:),gh1c(:,:),gh_direc(:,:),gresid(:,:),gvnlx1_saved(:,:)
308 4629309 : real(dp),allocatable :: gs1c(:,:),gvnlx_direc(:,:),pcon(:),sconjgr(:,:)
309 4629309 : real(dp),allocatable :: scprod(:,:),work(:,:),work1(:,:),work2(:,:)
310 4629309 : real(dp),pointer :: kinpw1(:)
311 4629309 : type(pawcprj_type),allocatable :: conjgrprj(:,:)
312 32405163 : type(pawcprj_type) :: cprj_dummy(1,1)
313 : ! *********************************************************************
314 :
315 : DBG_ENTER("COLL")
316 :
317 : ABI_NVTX_START_RANGE(NVTX_DFPT_CGWF)
318 :
319 4629309 : call timab(122,1,tsec)
320 :
321 : !==================================================================
322 : !========= LOCAL VARIABLES DEFINITIONS AND ALLOCATIONS ============
323 : !==================================================================
324 :
325 4629309 : 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 4629309 : nline = nline_in
338 4629309 : usetolrde__ = 1
339 4629309 : 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 4629309 : if (present(usetolrde)) usetolrde__ = usetolrde
345 :
346 4629309 : 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 4629309 : me_g0 = mpi_enreg%me_g0
353 4629309 : comm_fft = mpi_enreg%comm_fft
354 4629309 : me_band = mpi_enreg%me_band
355 4629309 : np_band = mpi_enreg%nproc_band
356 :
357 4629309 : skipme = 0
358 :
359 : ! if PAW, one has to solve a generalized eigenproblem
360 4629309 : usepaw=gs_hamkq%usepaw
361 4629309 : gen_eigenpb=(usepaw==1)
362 4629309 : useoverlap=0;if (gen_eigenpb) useoverlap=1
363 :
364 : ! Use scissor shift on 0-order eigenvalue
365 4629309 : eshift=eig0_k(u1_band)-dfpt_sciss
366 :
367 : ! Remove omega for a finite-frequency calculation
368 4629309 : omega_=zero ; if (present(omega)) omega_=omega
369 4629309 : eta_=zero ; if (present(eta)) eta_=eta
370 4629309 : eshift=eshift+omega_
371 :
372 : ! Additional initializations
373 4629309 : istwf_k=gs_hamkq%istwf_k
374 4629309 : optekin=0;if (wfoptalg>=10) optekin=1
375 4629309 : tol_restart=tol12;if (gen_eigenpb) tol_restart=tol8
376 4629309 : if (ipert == natom+10 .or. ipert == natom+11) tol_restart = tol7
377 : !tol_restart = tol7
378 :
379 4629309 : kinpw1 => gs_hamkq%kinpw_kp
380 :
381 : ! Memory allocations
382 13887927 : ABI_MALLOC(gh1c,(2,npw1*nspinor))
383 13887927 : ABI_MALLOC(pcon,(npw1))
384 13887927 : ABI_MALLOC(scprod,(2,nband_me))
385 :
386 4629309 : 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 4626205 : 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 4629309 : dc_shift_band=(u1_band-1)*npw1*nspinor
395 :
396 : ! this is used many times - no use de and re allocating
397 9258618 : 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 4629309 : 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 4629309 : 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 4569597 : if (gen_eigenpb) then
569 589604 : sij_opt=1
570 1179208 : ABI_MALLOC(gs1c,(2,npw1*nspinor))
571 : else
572 3979993 : ABI_MALLOC(gs1c,(0,0))
573 3979993 : sij_opt=0
574 : end if
575 4569597 : usevnl=1; optlocal=1; optnl=2
576 4569597 : 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 9139194 : mpi_enreg,1,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
582 :
583 4569597 : if (gen_eigenpb) then
584 589604 : if (ipert/=natom+2) then ! S^(1) is zero for ipert=natom+2
585 : !$OMP PARALLEL
586 : !$OMP DO
587 74758278 : do ipw=1,npw1*nspinor
588 223205386 : gh1c (1:2,ipw)=gh1c (1:2,ipw)-eshift*gs1c(1:2,ipw)
589 : end do
590 : !$OMP END DO NOWAIT
591 534724 : if (opt_gvnlx1/=1) then
592 : !$OMP DO
593 74758278 : do ipw=1,npw1*nspinor
594 223205386 : 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 589604 : 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 534724 : & 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 4629309 : 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 4629309 : 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 1716704467 : eig1_k_loc = zero
666 41160019 : do iband = 1, nband
667 36530710 : if (bands_treated_now(iband) == 0) cycle
668 3140849256 : if (rank_band(iband) == me_band) work = gh1c
669 4635357 : 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 4635357 : if(gen_eigenpb)then
675 : call projbd(gscq,work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
676 596612 : 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 4635357 : 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 4635357 : jband_me = 0
690 41311219 : do jband=1,nband
691 36675862 : if (rank_band(jband) /= me_band) cycle
692 36530710 : jband_me = jband_me + 1
693 114227487 : 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 18535380 : if (rank_band(iband) == me_band) then
699 3140843208 : gh1c = work - (mpi_enreg%nproc_band-1)*gh1c
700 : end if
701 :
702 : end do !iband
703 :
704 :
705 4629309 : 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 4569597 : if (gen_eigenpb) then
715 4626252 : do iband=1,nband
716 4036648 : if (bands_treated_now(iband) == 0) cycle
717 246639406 : 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 592772 : 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 592772 : indx_cgq=icgq
723 5295056 : do jband=1,nband
724 4112680 : if (rank_band(jband) /= me_band) cycle
725 :
726 4036648 : 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 4036648 : me_g0,mpi_enreg%comm_spinorfft)
729 4036648 : eig1_k_loc(1,jband,iband)=eig1_k_loc(1,jband,iband)-eshiftkq*dotr
730 4036648 : eig1_k_loc(2,jband,iband)=eig1_k_loc(2,jband,iband)-eshiftkq*doti
731 8149328 : 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 4569597 : 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 4569597 : call xmpi_sum(eig1_k_loc, mpi_enreg%comm_band, ierr)
742 :
743 : ! TODO: I think this is just a reshape
744 40807059 : do iband=1,nband
745 36237462 : if (bands_treated_now(iband) == 0) cycle
746 4575645 : band_off=(iband-1)*2*nband
747 45527856 : do jband=1,nband
748 36382614 : eig1_k(2*jband-1+band_off) = eig1_k_loc(1,jband,iband)
749 72620076 : 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 9402056 : do ispinor=1,nspinor
757 4772747 : ipws=(ispinor-1)*npw1
758 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef,kinpw1,ipws,npw1)
759 1053263586 : do ipw=1+ipws,npw1+ipws
760 1048634277 : 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 41160019 : do iband = 1, nband
770 36530710 : if (bands_treated_now(iband) == 0) cycle
771 :
772 3140849256 : if (rank_band(iband) == me_band) work = cwavef
773 4635357 : 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 4635357 : 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 4635357 : 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 13900023 : if (rank_band(iband) == me_band) then
789 3136213899 : cwavef = work - (mpi_enreg%nproc_band-1)*cwavef
790 : end if
791 : end do
792 :
793 :
794 4629309 : 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 4569597 : if (usepaw==1.and.usedcwavef>0) then
797 : !$OMP PARALLEL DO
798 74758278 : do ipw=1,npw1*nspinor
799 223205386 : 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 4629309 : if (u1_band>max(1,nband-nbdbuf))then
810 : ! Treat the case of buffer bands
811 276271836 : cwavef=zero
812 276271836 : ghc =zero
813 276271836 : gvnlxc =zero
814 19409754 : if (gen_eigenpb) gsc=zero
815 17532954 : if (usedcwavef==2) dcwavef=zero
816 395802 : if (usepaw==1) then
817 24052 : call pawcprj_set_zero(cwaveprj)
818 : end if
819 17532954 : if (usedcwavef==2) dcwavef=zero
820 : ! Number of one-way 3D ffts skipped
821 395802 : 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 395802 : 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 9258618 : ABI_MALLOC(conjgr,(2,npw1*nspinor))
835 9258618 : ABI_MALLOC(direc,(2,npw1*nspinor))
836 9258618 : ABI_MALLOC(gresid,(2,npw1*nspinor))
837 9258618 : ABI_MALLOC(cwaveq,(2,npw1*nspinor))
838 4629309 : if (usepaw==1) then
839 4287880 : ABI_MALLOC(conjgrprj,(natom,nspinor))
840 593444 : call pawcprj_alloc(conjgrprj,0,gs_hamkq%dimcprj)
841 : else
842 4035865 : ABI_MALLOC(conjgrprj,(0,0))
843 : end if
844 :
845 3136213899 : cwaveq(:,:)=cgq(:,1+npw1*nspinor*(band_me-1)+icgq:npw1*nspinor*band_me+icgq)
846 4629309 : dotgp=one
847 :
848 : ! Here apply H(0) at k+q to input orthogonalized 1st-order wfs
849 4629309 : sij_opt=0;if (gen_eigenpb) sij_opt=1
850 4629309 : cpopt=-1+usepaw
851 : call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamkq,gvnlxc,eshift,mpi_enreg,1,&
852 4629309 : prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
853 :
854 : ! ghc also includes the eigenvalue shift
855 4629309 : if (gen_eigenpb) then
856 1780332 : 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 4629309 : resid=zero
863 :
864 41160019 : bands_skipped_now = 0
865 4629309 : bands_skipped_now(u1_band) = skipme
866 4629309 : 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 17261117 : 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 15670759 : 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 4458550400 : do ipw=1,npw1*nspinor
897 13344329112 : 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 136535616 : do iband = 1, nband
915 120864857 : if (bands_treated_now(iband) == 0) cycle
916 13361244623 : if (rank_band(iband) == me_band) work = gresid
917 15694951 : 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 15694951 : if(gen_eigenpb)then
923 : call projbd(gscq,work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband_me,npw1,nspinor,&
924 1834698 : 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 15694951 : 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 47060661 : if (rank_band(iband) == me_band) then
938 13361220431 : gresid = work - (mpi_enreg%nproc_band-1)*gresid
939 : end if
940 : end do
941 :
942 15670759 : 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 15670759 : call dotprod_g(prod1,doti,istwf_k,npw1*nspinor,1,cwavef,gresid,me_g0,mpi_enreg%comm_spinorfft)
950 15670759 : call dotprod_g(prod2,doti,istwf_k,npw1*nspinor,1,cwavef,gh1c,me_g0,mpi_enreg%comm_spinorfft)
951 15670759 : 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 15670759 : u1h0me0u1=-prod1-prod2
966 :
967 : ! Some tolerance is allowed, to account for very small numerical inaccuracies and cancellations.
968 15670759 : 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 15670581 : if (skipme == 0) then
988 : ! Compute residual (squared) norm
989 15278749 : call sqnorm_g(resid,istwf_k,npw1*nspinor,gresid,me_g0,comm_fft)
990 15278749 : 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 15670759 : if (resid<tolwfr .and. skipme == 0) then
1000 705413 : 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 705413 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
1006 705413 : skipme = 1
1007 : end if
1008 :
1009 : ! If user require exiting the job, stop line minimisations
1010 15670759 : 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 15670759 : if (iline/=1) then
1019 11077018 : 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 15670759 : 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 15670759 : if (wfoptalg<10.or.iline==1) then
1047 : !print *, 'cwaveq cg_precon ', cwaveq; print *, 'direc cg_precon ', direc
1048 14432097 : call cg_precon(cwaveq,zero,istwf_k,kinpw1,npw1,nspinor,me_g0,0,pcon,direc,mpi_enreg%comm_fft)
1049 : else
1050 2532700 : do ispinor=1,nspinor
1051 1294038 : igs=(ispinor-1)*npw1
1052 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(igs,npw,direc,pcon)
1053 167273497 : do ipw=1+igs,npw1+igs
1054 495516429 : 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 136535616 : do iband = 1, nband
1072 120864857 : if (bands_treated_now(iband) == 0) cycle
1073 :
1074 13361244623 : if (rank_band(iband) == me_band) work = direc
1075 15694951 : 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 15694951 : 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 15694951 : 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 47060661 : if (rank_band(iband) == me_band) then
1091 13361220431 : 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 15670759 : call dotprod_g(dotgg,doti,istwf_k,npw1*nspinor,1,direc,gresid,me_g0,mpi_enreg%comm_spinorfft)
1105 :
1106 15670759 : if (iline==1) then
1107 : ! At first iteration, gamma is set to zero
1108 4593741 : gamma=zero
1109 4593741 : dotgp=dotgg
1110 4593741 : call cg_zcopy(npw1*nspinor,direc,conjgr)
1111 : else
1112 : ! At next iterations, h = g + gamma * h
1113 11077018 : gamma=dotgg/dotgp
1114 11077018 : dotgp=dotgg
1115 11077018 : 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 3413389725 : do ipw=1,npw1*nspinor
1121 10218015139 : conjgr(1:2,ipw)=direc(1:2,ipw)+gamma*conjgr(1:2,ipw)
1122 : end do
1123 11077018 : 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 15670759 : call dotprod_g(dedt,doti,istwf_k,npw1*nspinor,1,conjgr,gresid,me_g0,mpi_enreg%comm_spinorfft)
1135 15670759 : dedt=-two*two*dedt
1136 :
1137 15670759 : 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 47012277 : ABI_MALLOC(gvnlx_direc,(2,npw1*nspinor))
1141 31341518 : ABI_MALLOC(gh_direc,(2,npw1*nspinor))
1142 15670759 : if (gen_eigenpb) then
1143 3644052 : ABI_MALLOC(sconjgr,(2,npw1*nspinor))
1144 : else
1145 13848733 : ABI_MALLOC(sconjgr,(0,0))
1146 : end if
1147 15670759 : sij_opt=0;if (gen_eigenpb) sij_opt=1
1148 15670759 : cpopt=-1+usepaw
1149 : call getghc(cpopt,conjgr,conjgrprj,gh_direc,sconjgr,gs_hamkq,gvnlx_direc,&
1150 15670759 : eshift,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
1151 :
1152 : ! ghc also includes the eigenvalue shift
1153 15670759 : if (gen_eigenpb) then
1154 : !$OMP PARALLEL DO
1155 247432993 : do ipw=1,npw1*nspinor
1156 738654927 : 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 15670759 : call dotprod_g(d2edt2,doti,istwf_k,npw1*nspinor,1,conjgr,gh_direc,me_g0,mpi_enreg%comm_spinorfft)
1172 15670759 : d2edt2=two*two*d2edt2
1173 :
1174 15670759 : 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 15670759 : 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 15669448 : else if (d2edt2 > 1.d-40) then
1209 : ! Here, the value of theta that gives the minimum
1210 15669126 : 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 15670759 : 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 15669126 : if (skipme == 0) then
1233 43716075 : call cg_zaxpy(npw1*nspinor, [theta, zero], conjgr,cwavef)
1234 : ! Filter the wavefunctions for large modified kinetic energy (see routine mkkin.f)
1235 29663796 : do ispinor=1,nspinor
1236 15091771 : ipws=(ispinor-1)*npw1
1237 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef,kinpw1,ipws,npw1)
1238 4118130471 : do ipw=1+ipws,npw1+ipws
1239 4103558446 : 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 43716075 : call cg_zaxpy(npw1*nspinor, [theta, zero], gh_direc,ghc)
1246 43716075 : call cg_zaxpy(npw1*nspinor, [theta, zero], gvnlx_direc,gvnlxc)
1247 :
1248 17982625 : if (gen_eigenpb) call cg_zaxpy(npw1*nspinor, [theta, zero], sconjgr, gsc)
1249 1705300 : if (usepaw==1) call pawcprj_axpby(theta,one,conjgrprj,cwaveprj)
1250 : end if
1251 :
1252 15670759 : ABI_FREE(gh_direc)
1253 15670759 : ABI_FREE(gvnlx_direc)
1254 15670759 : ABI_FREE(sconjgr)
1255 :
1256 : ! ======================================================================
1257 : ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
1258 : ! ======================================================================
1259 :
1260 15670759 : if (usetolrde__ /= 0) then
1261 : ! Check reduction in trial energy deltae, Eq.(28) of PRB55, 10337 (1997) [[cite:Gonze1997]]
1262 13244660 : deltae=half*d2edt2*theta**2+theta*dedt
1263 :
1264 13244660 : if (iline==1) then
1265 4429599 : deold=deltae
1266 : ! The extra factor of two should be removed !
1267 8815061 : else if (abs(deltae)<tolrde*two*abs(deold) .and. iline/=nline) then
1268 1974379 : 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 1974379 : nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
1274 1974379 : 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 136535616 : bands_skipped_now = 0
1282 15670759 : bands_skipped_now(u1_band) = skipme
1283 15670759 : if (u1_band_ < 0) bands_skipped_now(u1_band) = 0 ! Handle negative index
1284 15670759 : call xmpi_sum(bands_skipped_now,mpi_enreg%comm_band,ierr)
1285 :
1286 : ! bands_skipped_now = bands_skipped_now - bands_treated_now
1287 136535616 : if (sum(abs(bands_skipped_now - bands_treated_now)) == 0) exit
1288 :
1289 17261117 : 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 4629309 : 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 4629309 : 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 4629309 : 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 4629309 : 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 4629309 : ABI_SFREE(gh_direc)
1530 4629309 : ABI_SFREE(gvnlx_direc)
1531 :
1532 4629309 : ABI_FREE(conjgr)
1533 4629309 : ABI_FREE(cwaveq)
1534 4629309 : ABI_FREE(direc)
1535 4629309 : ABI_FREE(gresid)
1536 4629309 : if (usepaw==1) call pawcprj_free(conjgrprj)
1537 5923125 : ABI_FREE(conjgrprj)
1538 :
1539 4629309 : 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 395802 : 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 4629309 : 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 4629309 : ABI_FREE(work)
1555 4629309 : ABI_FREE(gh1c)
1556 4629309 : ABI_FREE(pcon)
1557 4629309 : ABI_FREE(scprod)
1558 4629309 : ABI_FREE(gberry)
1559 :
1560 4629309 : call timab(122,2,tsec)
1561 :
1562 : ABI_NVTX_END_RANGE()
1563 :
1564 : DBG_EXIT("COLL")
1565 :
1566 13887927 : 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 : !!***
|
