Line data Source code
1 : !!****m* ABINIT/m_cgwf_cprj
2 : !! NAME
3 : !! m_cgwf_cprj
4 : !!
5 : !! FUNCTION
6 : !! Conjugate-gradient eigensolver.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (DCA, XG, GMR, MT, MVeithen, ISouza, JIniguez)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_cgwf_cprj
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_xmpi
27 : use m_abicore
28 : use m_cgtools
29 :
30 : use defs_abitypes, only : MPI_type
31 : use m_time, only : timab
32 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, &
33 : pawcprj_free,&
34 : pawcprj_axpby,pawcprj_zaxpby,pawcprj_projbd
35 : use m_hamiltonian, only : gs_hamiltonian_type
36 : use m_getchc, only : getchc,getcsc
37 : use m_getghc, only : getghc
38 : use m_nonlop, only : nonlop
39 : use m_paw_overlap, only : smatrix_k_paw
40 : use m_cgprj, only : getcprj
41 : use m_fft, only : fourwf
42 : use m_dtset, only : dataset_type
43 :
44 : implicit none
45 :
46 : private
47 : !!***
48 :
49 : public :: cgwf_cprj
50 : public :: mksubovl
51 : public :: cprj_update
52 : public :: cprj_update_oneband
53 : public :: cprj_check
54 : public :: get_cprj_id
55 :
56 : !!***
57 :
58 : contains
59 : !!***
60 :
61 : !!****f* m_cgwf/cgwf_cprj
62 : !! NAME
63 : !! cgwf_cprj
64 : !!
65 : !! FUNCTION
66 : !! Update all wavefunction |C>, non self-consistently.
67 : !! Uses a conjugate-gradient algorithm.
68 : !! This version is available for PAW only, as it needs the cprj in memory.
69 : !!
70 : !! INPUTS
71 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
72 : !! icg=shift to be applied on the location of data in the array cg
73 : !! mcg=second dimension of the cg array
74 : !! mpi_enreg=information about MPI parallelization
75 : !! nband=number of bands.
76 : !! nline=number of line minimizations per band.
77 : !! ortalg=governs the choice of the algorithm for orthogonalisation.
78 : !! prtvol=control print volume and debugging output
79 : !! tolrde=tolerance on the ratio of differences of energies (for the line minimisation)
80 : !! tolwfr=tolerance on largest wf residual
81 : !! wfoptalg=govern the choice of algorithm for wf optimisation
82 : !! (0, 1, 10 and 11 : in the present routine, usual CG algorithm ;
83 : !!
84 : !! OUTPUT
85 : !! resid(nband)=wf residual for new states=|(H-e)|C>|^2 (hartree^2)
86 : !!
87 : !! SIDE EFFECTS
88 : !! cg(2,mcg)
89 : !! at input =wavefunction <G|C band,k> coefficients for ALL bands
90 : !! at output same as input except that
91 : !! the current band, with number 'band' has been updated
92 : !! cprj_cwavef_bands=<p_i|c_n> coefficients for all bands n, they are updated when the WFs change.
93 : !! quit= if 1, proceeds to smooth ending of the job.
94 : !!
95 : !! NOTES
96 : !!
97 : !! SOURCE
98 :
99 0 : subroutine cgwf_cprj(cg,cprj_cwavef_bands,cprj_update_lvl,eig,&
100 0 : & gs_hamk,icg,mcg,mpi_enreg,nband,nline,ortalg,prtvol,quit,resid,subham,&
101 : & tolrde,tolwfr,wfoptalg)
102 : !Arguments ------------------------------------
103 : integer,intent(in) :: cprj_update_lvl,icg
104 : integer,intent(in) :: mcg,nband,nline,ortalg,prtvol
105 : integer,intent(in) :: wfoptalg,quit
106 : real(dp),intent(in) :: tolrde,tolwfr
107 : type(MPI_type),intent(in) :: mpi_enreg
108 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
109 : !arrays
110 : real(dp),intent(inout),target :: cg(2,mcg)
111 : real(dp), intent(inout) :: eig(:)
112 : real(dp),intent(out) :: resid(:),subham(:)
113 : type(pawcprj_type),intent(inout),target :: cprj_cwavef_bands(:,:)
114 :
115 : !Local variables-------------------------------
116 : integer,parameter :: level=113,tim_getghc=1,tim_projbd=1,type_calc=0,tim_getcprj=3
117 : integer,parameter :: useoverlap=0,tim_getcsc=3
118 : integer,save :: nskip=0
119 : integer :: cpopt,i1,i2,i3,iband,isubh,isubh0,jband,me_g0,igs
120 : integer :: iline,ipw,ispinor,istwf_k
121 : integer :: n4,n5,n6,natom,ncpgr,npw,nspinor
122 : integer :: optekin,sij_opt,wfopta10
123 : real(dp) :: chc,costh,deltae,deold,dhc,dhd,diff,dotgg,dotgp,doti,dotr,eval,gamma
124 : real(dp) :: lam0,lamold,root,sinth,sintn,swap,tan2th,xnorm,xnormd,xnormd_previous
125 : character(len=500) :: message
126 : !arrays
127 : real(dp) :: dot(2)
128 : real(dp) :: tsec(2)
129 0 : real(dp),allocatable :: conjgr(:,:),gvnlxc(:,:)
130 0 : real(dp), pointer :: cwavef(:,:),cwavef_left(:,:),cwavef_bands(:,:)
131 0 : real(dp), allocatable :: cwavef_r(:,:,:,:,:)
132 0 : real(dp),allocatable,target :: direc(:,:)
133 0 : real(dp),allocatable :: direc_tmp(:,:),pcon(:),scprod(:,:),scwavef_dum(:,:)
134 0 : real(dp),allocatable :: direc_r(:,:,:,:,:),scprod_csc(:)
135 : real(dp) :: z_tmp(2),z_tmp2(2)
136 0 : type(pawcprj_type),pointer :: cprj_cwavef(:,:),cprj_cwavef_left(:,:)
137 0 : type(pawcprj_type),allocatable :: cprj_direc(:,:),cprj_conjgr(:,:)
138 :
139 : ! *********************************************************************
140 :
141 : DBG_ENTER("COLL")
142 :
143 : !Starting the routine
144 0 : call timab(1300,1,tsec)
145 :
146 : !Some checks
147 : !Only wfoptalg==10 for now
148 0 : if (wfoptalg/=10) then
149 0 : ABI_ERROR("cgwf_cprj is implemented for wfoptalg==10 only")
150 : end if
151 :
152 : !======================================================================
153 : !========= LOCAL VARIABLES DEFINITIONS AND ALLOCATIONS ================
154 : !======================================================================
155 :
156 : !MPI data
157 0 : me_g0 = mpi_enreg%me_g0
158 :
159 : ! cprj are already in memory
160 0 : cpopt = 2
161 :
162 : !Initializations and allocations
163 0 : istwf_k=gs_hamk%istwf_k
164 0 : wfopta10=mod(wfoptalg,10)
165 :
166 0 : optekin=0;if (wfoptalg>=10) optekin=1
167 0 : natom=gs_hamk%natom
168 0 : npw=gs_hamk%npw_k
169 0 : nspinor=gs_hamk%nspinor
170 :
171 0 : ABI_MALLOC(pcon,(npw))
172 0 : ABI_MALLOC(conjgr,(2,npw*nspinor))
173 0 : ABI_MALLOC(scwavef_dum,(0,0))
174 0 : ABI_MALLOC(direc,(2,npw*nspinor))
175 0 : ABI_MALLOC(scprod,(2,nband))
176 0 : ABI_MALLOC(scprod_csc,(2*nband))
177 0 : ABI_MALLOC(direc_tmp,(2,npw*nspinor))
178 0 : ABI_MALLOC(gvnlxc,(0,0))
179 :
180 0 : ABI_MALLOC(cprj_direc ,(natom,nspinor))
181 0 : ABI_MALLOC(cprj_conjgr ,(natom,nspinor))
182 :
183 0 : n4=gs_hamk%ngfft(4);n5=gs_hamk%ngfft(5);n6=gs_hamk%ngfft(6)
184 0 : ABI_MALLOC(cwavef_r,(2,n4,n5,n6,nspinor))
185 0 : ABI_MALLOC(direc_r, (2,n4,n5,n6,nspinor))
186 :
187 0 : ncpgr = 0 ! no need of gradients here...
188 0 : call pawcprj_alloc(cprj_direc,ncpgr,gs_hamk%dimcprj)
189 0 : call pawcprj_alloc(cprj_conjgr,ncpgr,gs_hamk%dimcprj)
190 :
191 0 : cwavef_bands => cg(:,1+icg:nband*npw*nspinor+icg)
192 :
193 0 : if (cprj_update_lvl==-2)then
194 0 : call cprj_update(cg,cprj_cwavef_bands,gs_hamk,icg,nband,mpi_enreg,tim_getcprj)
195 : endif
196 :
197 0 : isubh=1
198 0 : isubh0=1
199 :
200 : ! Big iband loop
201 0 : do iband=1,nband
202 :
203 : ! ======================================================================
204 : ! ========== INITIALISATION OF MINIMIZATION ITERATIONS =================
205 : ! ======================================================================
206 :
207 0 : if (prtvol>=10) then ! Tell us what is going on:
208 0 : write(message, '(a,i6,2x,a,i3,a)' )' --- cgwf is called for band',iband,'for',nline,' lines'
209 0 : call wrtout(std_out,message,'PERS')
210 : end if
211 :
212 0 : dotgp=one
213 :
214 : ! Extraction of the vector that is iteratively updated
215 0 : cwavef => cwavef_bands(:,1+(iband-1)*npw*nspinor:iband*npw*nspinor)
216 0 : cprj_cwavef => cprj_cwavef_bands(:,nspinor*(iband-1)+1:nspinor*iband)
217 :
218 : ! Normalize incoming wf:
219 : call getcsc(dot,cpopt,cwavef,cwavef,cprj_cwavef,cprj_cwavef,&
220 0 : & gs_hamk,mpi_enreg,1,tim_getcsc)
221 0 : xnorm=one/sqrt(dot(1))
222 0 : z_tmp = (/xnorm,zero/)
223 : ! cwavef = xnorm * cwavef
224 0 : call timab(1305,1,tsec)
225 0 : call cg_zscal(npw*nspinor,z_tmp,cwavef)
226 0 : call timab(1305,2,tsec)
227 : ! cprj = xnorm * cprj
228 0 : call timab(1302,1,tsec)
229 0 : call pawcprj_axpby(zero,xnorm,cprj_cwavef,cprj_cwavef)
230 0 : call timab(1302,2,tsec)
231 :
232 : ! Compute wavefunction in real space
233 0 : call get_cwavefr(cwavef,cwavef_r,gs_hamk,mpi_enreg)
234 :
235 0 : if (prtvol==-level) then
236 0 : write(message,'(a,f26.14)')' cgwf: xnorm = ',xnorm
237 0 : call wrtout(std_out,message,'PERS')
238 : end if
239 :
240 : ! ======================================================================
241 : ! ====== BEGIN LOOP FOR A GIVEN BAND: MINIMIZATION ITERATIONS ==========
242 : ! ======================================================================
243 0 : if(nline/=0)then
244 0 : do iline=1,nline
245 :
246 : ! === COMPUTE THE RESIDUAL ===
247 : ! Compute lambda = <C|H|C>
248 0 : sij_opt=0
249 : call getchc(z_tmp,cpopt,cwavef,cwavef,cprj_cwavef,cprj_cwavef,cwavef_r,cwavef_r,&
250 0 : & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc)
251 0 : chc=z_tmp(1)
252 0 : lam0=chc
253 0 : eval=chc
254 0 : eig(iband)=chc
255 :
256 : ! Check that lam0 is decreasing on succeeding lines:
257 0 : if (iline==1) then
258 0 : lamold=lam0
259 : else
260 0 : if (lam0 > lamold+tol12) then
261 : write(message, '(a,i8,a,1p,e14.6,a1,3x,a,1p,e14.6,a1)')&
262 0 : & 'New trial energy at line ',iline,' = ',lam0,ch10,&
263 0 : & 'is higher than former =',lamold,ch10
264 0 : ABI_WARNING(message)
265 : end if
266 0 : lamold=lam0
267 : end if
268 :
269 : ! ======================================================================
270 : ! =========== COMPUTE THE STEEPEST DESCENT DIRECTION ===================
271 : ! ======================================================================
272 :
273 0 : sij_opt=-1
274 : call getghc(cpopt,cwavef,cprj_cwavef,direc,scwavef_dum,gs_hamk,gvnlxc,&
275 0 : & eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc,cwavef_r=cwavef_r)
276 :
277 : ! Compute residual (squared) norm
278 0 : call timab(1305,1,tsec)
279 0 : call sqnorm_g(resid(iband),istwf_k,npw*nspinor,direc,me_g0,mpi_enreg%comm_fft)
280 0 : call timab(1305,2,tsec)
281 :
282 0 : if (prtvol==-level) then
283 0 : write(message,'(a,i0,f26.14,es27.14e3)')' cgwf: iline,eval,resid = ',iline,eval,resid(iband)
284 0 : call wrtout(std_out,message,'PERS')
285 : end if
286 0 : xnormd=1/sqrt(resid(iband))
287 0 : z_tmp = (/xnormd,zero/)
288 0 : call cg_zscal(npw*nspinor,z_tmp,direc)
289 :
290 : ! ======================================================================
291 : ! ============== CHECK FOR CONVERGENCE CRITERIA ========================
292 : ! ======================================================================
293 :
294 : ! If residual sufficiently small stop line minimizations
295 0 : if (resid(iband)<tolwfr) then
296 0 : if (prtvol>=10) then
297 : write(message, '(a,i4,a,i2,a,es12.4)' ) &
298 0 : & ' cgwf: band ',iband,' converged after ',iline,' line minimizations: resid =',resid(iband)
299 0 : call wrtout(std_out,message,'PERS')
300 : end if
301 0 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
302 0 : exit ! Exit from the loop on iline
303 : end if
304 :
305 : ! If user require exiting the job, stop line minimisations
306 0 : if (quit==1) then
307 0 : write(message, '(a,i0)' )' cgwf: user require exiting => skip update of band ',iband
308 0 : call wrtout(std_out,message,'PERS')
309 :
310 0 : nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped
311 0 : exit ! Exit from the loop on iline
312 : end if
313 :
314 : ! =========== PROJECT THE STEEPEST DESCENT DIRECTION ===================
315 : ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
316 :
317 : ! The following projection over the subspace orthogonal to occupied bands
318 : ! is optional. It is a bit more accurate, but doubles the number of N^3 ops.
319 : ! It is done only if ortalg>=0.
320 :
321 : ! Project the steepest descent direction:
322 : ! direc(2,npw)=<G|H|Cnk> - \sum_{(i/=n)} <G|H|Cik>
323 0 : if(ortalg>=0)then
324 0 : call cprj_update_oneband(direc,cprj_direc,gs_hamk,mpi_enreg,tim_getcprj)
325 : call getcsc(scprod_csc,cpopt,direc,cwavef_bands,cprj_direc,cprj_cwavef_bands,&
326 0 : & gs_hamk,mpi_enreg,nband,tim_getcsc)
327 0 : scprod = reshape(scprod_csc,(/2,nband/))
328 : ! Note that the current band (|C_iband>) is not used here
329 : call projbd(cg,direc,iband,icg,icg,istwf_k,mcg,mcg,nband,npw,nspinor,&
330 0 : & direc,scprod,1,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
331 : end if
332 :
333 : ! For a generalized eigenpb, store the steepest descent direction
334 0 : direc_tmp=direc
335 :
336 : ! ======================================================================
337 : ! ======== PRECONDITION THE STEEPEST DESCENT DIRECTION =================
338 : ! ======================================================================
339 :
340 : ! If wfoptalg>=10, the precondition matrix is kept constant during iteration ; otherwise it is recomputed
341 0 : call timab(1305,1,tsec)
342 0 : if (wfoptalg<10.or.iline==1) then
343 0 : call cg_precon(cwavef,zero,istwf_k,gs_hamk%kinpw_k,npw,nspinor,me_g0,optekin,pcon,direc,mpi_enreg%comm_fft)
344 : else
345 0 : do ispinor=1,nspinor
346 0 : igs=(ispinor-1)*npw
347 : !$OMP PARALLEL DO
348 0 : do ipw=1+igs,npw+igs
349 0 : direc(1,ipw)=direc(1,ipw)*pcon(ipw-igs)
350 0 : direc(2,ipw)=direc(2,ipw)*pcon(ipw-igs)
351 : end do
352 : end do
353 : end if
354 0 : call timab(1305,2,tsec)
355 :
356 : ! ======= PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION ==============
357 : ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
358 0 : call cprj_update_oneband(direc,cprj_direc,gs_hamk,mpi_enreg,tim_getcprj)
359 : call getcsc(scprod_csc,cpopt,direc,cwavef_bands,cprj_direc,cprj_cwavef_bands,&
360 0 : & gs_hamk,mpi_enreg,nband,tim_getcsc)
361 : ! if (abs(xnorm-one)>tol10) then ! True if iline==1 and if input WFs are random
362 0 : if (iline==1) then
363 : ! We compensate the normalization of the current band
364 0 : scprod_csc(2*iband-1:2*iband) = scprod_csc(2*iband-1:2*iband)/xnorm
365 : end if
366 : ! Projecting again out all bands (not normalized).
367 0 : scprod = reshape(scprod_csc,(/2,nband/))
368 : call projbd(cg,direc,-1,icg,icg,istwf_k,mcg,mcg,nband,npw,nspinor,&
369 0 : & direc,scprod,1,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
370 0 : if (iline==1) then
371 : ! Again we have to compensate the normalization of the current band.
372 : ! Indeed, by calling projbd we compute:
373 : ! |direc'_i> = |direc_i> - \sum_j <c'_j|S|direc_i>|c'_j>
374 : ! = |direc_i> - \sum_{j/=i} <c_j|S|direc_i>|c_j> - <c'_i|S|direc_i>|c'_i>
375 : ! where |c'_j> = |c_j> for j/=i and |c'_i> = xnorm.|c_i>
376 : ! As we compensated "scprod" before the call of projbd we actually computed:
377 : ! |direc'_i> = |direc_i> - \sum_{j/=i} <c_j|S|direc_i>|c_j> - <c_i|S|direc_i>|c'_i>
378 : ! = |direc_i> - \sum_{j/=i} <c_j|S|direc_i>|c_j> - xnorm.<c_i|S|direc_i>|c_i>
379 : ! The correct projected direction should be:
380 : ! |projdirec_i> = |direc_i> - \sum_j <c_j|S|direc_i>|c_j>
381 : ! = |direc_i> - \sum_{j/=i} <c_j|S|direc_i>|c_j> - <c_i|S|direc_i>|c_i>
382 : ! So:
383 : ! |projdirec_i> = |direc'_i> - <c_i|S|direc_i>|c_i> + xnorm.<c_i|S|direc_i>|c_i>
384 : ! = |direc'_i> - (1-xnorm) <c_i|S|direc_i>|c_i>
385 : ! = |direc'_i> - (1-xnorm)/xnorm <c_i|S|direc_i>|c'_i>
386 : !
387 0 : z_tmp = -scprod_csc(2*iband-1:2*iband)*(one-xnorm)/xnorm
388 0 : call timab(1305,1,tsec)
389 0 : do ispinor=1,nspinor
390 0 : igs=(ispinor-1)*npw
391 0 : do ipw=1+igs,npw+igs
392 0 : direc(1,ipw)=direc(1,ipw) + z_tmp(1)*cwavef(1,ipw) - z_tmp(2)*cwavef(2,ipw)
393 0 : direc(2,ipw)=direc(2,ipw) + z_tmp(1)*cwavef(2,ipw) + z_tmp(2)*cwavef(1,ipw)
394 : end do
395 : end do
396 0 : call timab(1305,2,tsec)
397 : end if
398 : ! Apply projbd to cprj_direc
399 0 : scprod=-scprod
400 0 : call timab(1303,1,tsec)
401 0 : call pawcprj_projbd(scprod,cprj_cwavef_bands,cprj_direc)
402 0 : call timab(1303,2,tsec)
403 0 : if (iline==1) then
404 : ! Same correction than for WFs
405 0 : z_tmp = -scprod_csc(2*iband-1:2*iband)*(one-xnorm)/xnorm
406 0 : z_tmp2 = (/one,zero/)
407 : ! cprj = z_tmp*cprjx + z_tmp2*cprj
408 0 : call timab(1302,1,tsec)
409 0 : call pawcprj_zaxpby(z_tmp,z_tmp2,cprj_cwavef,cprj_direc)
410 0 : call timab(1302,2,tsec)
411 : end if
412 :
413 : ! ======================================================================
414 : ! ================= COMPUTE THE CONJUGATE-GRADIENT =====================
415 : ! ======================================================================
416 :
417 0 : call timab(1305,1,tsec)
418 0 : call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,direc_tmp,me_g0,mpi_enreg%comm_spinorfft)
419 0 : call timab(1305,2,tsec)
420 :
421 0 : dotgg=dotgg/xnormd**2
422 :
423 : ! MJV: added 5 Feb 2012 - causes divide by 0 on next iteration of iline
424 0 : if (abs(dotgg) < TINY(0.0_dp)*1.e50_dp) dotgg = TINY(0.0_dp)*1.e50_dp
425 :
426 : ! At first iteration, gamma is set to zero
427 0 : if (iline==1) then
428 0 : gamma=zero
429 0 : dotgp=dotgg
430 0 : call timab(1305,1,tsec)
431 0 : call cg_zcopy(npw*nspinor,direc,conjgr)
432 0 : call timab(1305,2,tsec)
433 0 : call timab(1302,1,tsec)
434 0 : call pawcprj_copy(cprj_direc,cprj_conjgr)
435 0 : call timab(1302,2,tsec)
436 0 : if (prtvol==-level)then
437 0 : write(message,'(a,es27.14e3)')' cgwf: dotgg = ',dotgg
438 0 : call wrtout(std_out,message,'PERS')
439 : end if
440 : else
441 0 : gamma=dotgg/dotgp
442 0 : dotgp=dotgg
443 :
444 0 : if (prtvol==-level)then
445 0 : write(message,'(a,2es27.14e3)')' cgwf: dotgg,gamma = ',dotgg,gamma
446 0 : call wrtout(std_out,message,'PERS')
447 : end if
448 :
449 0 : gamma=gamma*xnormd/xnormd_previous
450 :
451 0 : call timab(1305,1,tsec)
452 : !$OMP PARALLEL DO
453 0 : do ipw=1,npw*nspinor
454 0 : conjgr(1,ipw)=direc(1,ipw)+gamma*conjgr(1,ipw)
455 0 : conjgr(2,ipw)=direc(2,ipw)+gamma*conjgr(2,ipw)
456 : end do
457 0 : call timab(1305,2,tsec)
458 0 : call timab(1302,1,tsec)
459 0 : call pawcprj_axpby(one,gamma,cprj_direc,cprj_conjgr)
460 0 : call timab(1302,2,tsec)
461 : end if
462 0 : call timab(1305,1,tsec)
463 :
464 0 : xnormd_previous=xnormd
465 :
466 : call getcsc(dot,cpopt,conjgr,conjgr,cprj_conjgr,cprj_conjgr,&
467 0 : & gs_hamk,mpi_enreg,1,tim_getcsc)
468 :
469 : ! ======================================================================
470 : ! ============ PROJECTION OF THE CONJUGATED GRADIENT ===================
471 : ! ======================================================================
472 :
473 : call getcsc(dot,cpopt,conjgr,cwavef,cprj_conjgr,cprj_cwavef,&
474 0 : & gs_hamk,mpi_enreg,1,tim_getcsc)
475 0 : dotr=dot(1)
476 0 : doti=dot(2)
477 :
478 : ! Project the conjugated gradient onto the current band
479 0 : call timab(1305,1,tsec)
480 0 : if(istwf_k==1)then
481 :
482 : !$OMP PARALLEL DO
483 0 : do ipw=1,npw*nspinor
484 0 : direc(1,ipw)=conjgr(1,ipw)-(dotr*cwavef(1,ipw)-doti*cwavef(2,ipw))
485 0 : direc(2,ipw)=conjgr(2,ipw)-(dotr*cwavef(2,ipw)+doti*cwavef(1,ipw))
486 : end do
487 : else
488 : !$OMP PARALLEL DO
489 0 : do ipw=1,npw*nspinor
490 0 : direc(1,ipw)=conjgr(1,ipw)-dotr*cwavef(1,ipw)
491 0 : direc(2,ipw)=conjgr(2,ipw)-dotr*cwavef(2,ipw)
492 : end do
493 : end if
494 0 : call timab(1305,2,tsec)
495 0 : call timab(1302,1,tsec)
496 0 : call pawcprj_copy(cprj_conjgr,cprj_direc)
497 0 : z_tmp = (/-dotr,-doti/)
498 0 : z_tmp2 = (/one,zero/)
499 0 : call pawcprj_zaxpby(z_tmp,z_tmp2,cprj_cwavef,cprj_direc)
500 0 : call timab(1302,2,tsec)
501 :
502 : ! ======================================================================
503 : ! ===== COMPUTE CONTRIBUTIONS TO 1ST AND 2ND DERIVATIVES OF ENERGY =====
504 : ! ======================================================================
505 :
506 : ! Compute norm of direc
507 : call getcsc(dot,cpopt,direc,direc,cprj_direc,cprj_direc,&
508 0 : & gs_hamk,mpi_enreg,1,tim_getcsc)
509 0 : xnormd=one/sqrt(abs(dot(1)))
510 :
511 0 : sij_opt=0
512 : ! Compute direc in real space
513 0 : call get_cwavefr(direc,direc_r,gs_hamk,mpi_enreg)
514 : ! Compute dhc = Re{<D|H|C>}
515 : call getchc(z_tmp,cpopt,cwavef,direc,cprj_cwavef,cprj_direc,cwavef_r,direc_r,&
516 0 : & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc)
517 0 : dhc=z_tmp(1)
518 0 : dhc=dhc*xnormd
519 :
520 : ! Compute <D|H|D> or <D|(H-zshift)^2|D>
521 : call getchc(z_tmp,cpopt,direc,direc,cprj_direc,cprj_direc,direc_r,direc_r,&
522 0 : & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc)
523 0 : dhd=z_tmp(1)
524 0 : dhd=dhd*xnormd**2
525 :
526 0 : if (prtvol==-level) then
527 0 : write(message,'(a,3f26.14)') 'cgwf: chc,dhc,dhd=',chc,dhc,dhd
528 0 : call wrtout(std_out,message,'PERS')
529 : end if
530 :
531 : ! ======================================================================
532 : ! ======= COMPUTE MIXING FACTORS - CHECK FOR CONVERGENCE ===============
533 : ! ======================================================================
534 :
535 : ! Compute tan(2 theta),sin(theta) and cos(theta)
536 0 : tan2th=2.0_dp*dhc/(chc-dhd)
537 :
538 0 : if (abs(tan2th)<1.d-05) then
539 0 : costh=1.0_dp-0.125_dp*tan2th**2
540 0 : sinth=0.5_dp*tan2th*(1.0_dp-0.375_dp*tan2th**2)
541 :
542 : ! Check that result is above machine precision
543 0 : if (abs(sinth)<epsilon(0._dp)) then
544 0 : write(message, '(a,es16.4)' ) ' cgwf: converged with tan2th=',tan2th
545 0 : call wrtout(std_out,message,'PERS')
546 : ! Number of one-way 3D ffts skipped
547 0 : nskip=nskip+2*(nline-iline)
548 0 : exit ! Exit from the loop on iline
549 : end if
550 :
551 : else
552 0 : root=sqrt(1.0_dp+tan2th**2)
553 0 : costh=sqrt(0.5_dp+0.5_dp/root)
554 0 : sinth=sign(sqrt(0.5_dp-0.5_dp/root),tan2th)
555 : end if
556 :
557 : ! Check for lower of two possible roots (same sign as curvature at theta where slope is zero)
558 0 : diff=(chc-dhd)
559 : ! Swap c and d if value of diff is positive
560 0 : if (diff>zero) then
561 0 : swap=costh
562 0 : costh=-sinth
563 0 : sinth=swap
564 0 : if(prtvol<0 .or. prtvol>=10)then
565 0 : write(message,*)' Note: swap roots, iline,diff=',iline,diff
566 0 : call wrtout(std_out,message,'PERS')
567 : end if
568 : end if
569 :
570 : ! ============================================================
571 : ! =========== GENERATE NEW wf(G), wf(r) and cprj =============
572 : ! ============================================================
573 :
574 0 : sintn=sinth*xnormd
575 :
576 0 : call timab(1305,1,tsec)
577 : !$OMP PARALLEL DO
578 0 : do ipw=1,npw*nspinor
579 0 : cwavef(1,ipw)=cwavef(1,ipw)*costh+direc(1,ipw)*sintn
580 0 : cwavef(2,ipw)=cwavef(2,ipw)*costh+direc(2,ipw)*sintn
581 : end do
582 0 : do ispinor=1,nspinor
583 0 : do i3=1,gs_hamk%n6
584 0 : do i2=1,gs_hamk%n5
585 0 : do i1=1,gs_hamk%n4
586 0 : cwavef_r(1,i1,i2,i3,ispinor)=cwavef_r(1,i1,i2,i3,ispinor)*costh+direc_r(1,i1,i2,i3,ispinor)*sintn
587 0 : cwavef_r(2,i1,i2,i3,ispinor)=cwavef_r(2,i1,i2,i3,ispinor)*costh+direc_r(2,i1,i2,i3,ispinor)*sintn
588 : end do
589 : end do
590 : end do
591 : end do
592 0 : call timab(1305,2,tsec)
593 0 : if (cprj_update_lvl<0) then
594 0 : call cprj_update_oneband(cwavef,cprj_cwavef,gs_hamk,mpi_enreg,tim_getcprj)
595 : else
596 0 : call timab(1302,1,tsec)
597 0 : call pawcprj_axpby(sintn,costh,cprj_direc,cprj_cwavef)
598 0 : call timab(1302,2,tsec)
599 : end if
600 :
601 : ! ======================================================================
602 : ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
603 : ! ======================================================================
604 :
605 : ! Compute delta(E)
606 0 : deltae=chc*(costh**2-1._dp)+dhd*sinth**2+2._dp*costh*sinth*dhc
607 :
608 : ! Check convergence and eventually exit
609 0 : if (iline==1) then
610 0 : deold=deltae
611 0 : else if (abs(deltae)<tolrde*abs(deold) .and. iline/=nline .and. wfopta10<2)then
612 0 : if(prtvol>=10)then
613 : write(message, '(a,i4,1x,a,1p,e12.4,a,e12.4,a)' ) &
614 0 : & ' cgwf: line',iline,&
615 0 : & ' deltae=',deltae,' < tolrde*',deold,' =>skip lines'
616 0 : call wrtout(std_out,message,'PERS')
617 : end if
618 : ! Update chc before exit
619 : call getchc(z_tmp,cpopt,cwavef,cwavef,cprj_cwavef,cprj_cwavef,cwavef_r,cwavef_r,&
620 0 : & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc)
621 0 : eig(iband)=z_tmp(1)
622 0 : nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
623 0 : exit ! Exit from the loop on iline
624 : end if
625 :
626 : ! Update chc only if last iteration, otherwise it will be done at the beginning of the next one
627 0 : if (iline==nline) then
628 : call getchc(z_tmp,cpopt,cwavef,cwavef,cprj_cwavef,cprj_cwavef,cwavef_r,cwavef_r,&
629 0 : & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc)
630 0 : eig(iband)=z_tmp(1)
631 : end if
632 :
633 : end do ! END LOOP FOR A GIVEN BAND Note that there are three "exit" instructions inside
634 :
635 : else ! nline==0 , needs to provide a residual
636 0 : resid(iband)=-one
637 : end if ! End nline==0 case
638 :
639 : ! ======================================================================
640 : ! =============== END OF CURRENT BAND: CLEANING ========================
641 : ! ======================================================================
642 :
643 : ! At the end of the treatment of a set of bands, write the number of one-way 3D ffts skipped
644 : if (xmpi_paral==0 .and. mpi_enreg%paral_kgb==0 .and. iband==nband .and. prtvol/=0) then
645 : write(message,'(a,i0)')' cgwf: number of one-way 3D ffts skipped in cgwf until now =',nskip
646 : call wrtout(std_out,message,'PERS')
647 : end if
648 :
649 : ! ======================================================================
650 : ! ============= COMPUTE HAMILTONIAN IN WFs SUBSPACE ====================
651 : ! ======================================================================
652 :
653 0 : if (cprj_update_lvl<=2) call cprj_update_oneband(cwavef,cprj_cwavef,gs_hamk,mpi_enreg,tim_getcprj)
654 :
655 : ! Compute local+kinetic part
656 0 : sij_opt=0
657 : call getghc(cpopt,cwavef,cprj_cwavef,direc,scwavef_dum,gs_hamk,gvnlxc,&
658 0 : & eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,3,cwavef_r=cwavef_r)
659 0 : isubh=isubh0
660 0 : call timab(1304,1,tsec)
661 0 : do jband=1,iband
662 0 : cwavef_left => cwavef_bands(:,1+(jband-1)*npw*nspinor:jband*npw*nspinor)
663 0 : call dotprod_g(subham(isubh),subham(isubh+1),istwf_k,npw*nspinor,2,cwavef_left,direc,me_g0,mpi_enreg%comm_spinorfft)
664 0 : isubh=isubh+2
665 : end do
666 0 : call timab(1304,2,tsec)
667 0 : cprj_cwavef_left => cprj_cwavef_bands(:,1:nspinor*iband)
668 : ! Add the nonlocal part
669 : call getchc(subham(isubh0:isubh0+2*iband-1),cpopt,cwavef,cwavef,&
670 : & cprj_cwavef,cprj_cwavef_left,cwavef_r,cwavef_r,&
671 0 : & gs_hamk,zero,mpi_enreg,iband,sij_opt,4)
672 0 : isubh0=isubh
673 :
674 : end do ! End big iband loop.
675 :
676 : ! ======================================================================
677 : ! ============= COMPUTE HAMILTONIAN IN WFs SUBSPACE ====================
678 : ! ======================================================================
679 :
680 : ! if (cprj_update_lvl<=2) call cprj_update(cg,cprj_cwavef_bands,gs_hamk,icg,nband,mpi_enreg,tim_getcprj)
681 : !
682 : ! sij_opt=0
683 : !
684 : ! isubh=1
685 : ! isubh0=1
686 : ! do iband=1,nband
687 : ! cwavef => cwavef_bands(:,1+(iband-1)*npw*nspinor:iband*npw*nspinor)
688 : ! cprj_cwavef => cprj_cwavef_bands(:,nspinor*(iband-1)+1:nspinor*iband)
689 : ! ! Compute local+kinetic part
690 : ! call getghc(cpopt,cwavef,cprj_cwavef,direc,scwavef_dum,gs_hamk,gvnlxc,&
691 : ! & eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,3)
692 : ! isubh=isubh0
693 : ! call timab(1304,1,tsec)
694 : ! do jband=1,iband
695 : ! cwavef_left => cwavef_bands(:,1+(jband-1)*npw*nspinor:jband*npw*nspinor)
696 : ! call dotprod_g(subham(isubh),subham(isubh+1),istwf_k,npw*nspinor,2,cwavef_left,direc,me_g0,mpi_enreg%comm_spinorfft)
697 : ! isubh=isubh+2
698 : ! end do
699 : ! call timab(1304,2,tsec)
700 : ! cprj_cwavef_left => cprj_cwavef_bands(:,1:nspinor*iband)
701 : ! ! Add the nonlocal part
702 : ! call getchc(subham(isubh0:isubh0+2*iband-1),cpopt,cwavef,cwavef,&
703 : ! & cprj_cwavef,cprj_cwavef_left,cwavef_r,cwavef_r,&
704 : ! & gs_hamk,zero,mpi_enreg,iband,sij_opt,4)
705 : ! isubh0=isubh
706 : ! end do
707 :
708 : ! Debugging ouputs
709 0 : if(prtvol==-level)then
710 0 : isubh=1
711 0 : do iband=1,nband
712 0 : do jband=1,iband
713 0 : if (jband<=10) then
714 0 : write(message,'(i7,2f26.14)')isubh,subham(isubh:isubh+1)
715 0 : call wrtout(std_out,message,'PERS')
716 : end if
717 0 : isubh=isubh+2
718 : end do
719 : end do
720 : end if
721 :
722 : ! ===================
723 : ! FINAL DEALLOCATIONS
724 : ! ===================
725 :
726 0 : nullify(cwavef_left)
727 0 : nullify(cwavef_bands)
728 0 : nullify(cwavef)
729 0 : nullify(cprj_cwavef)
730 0 : nullify(cprj_cwavef_left)
731 0 : call pawcprj_free(cprj_direc)
732 0 : call pawcprj_free(cprj_conjgr)
733 0 : ABI_FREE(cprj_direc)
734 0 : ABI_FREE(cprj_conjgr)
735 0 : ABI_FREE(conjgr)
736 0 : ABI_FREE(scwavef_dum)
737 0 : ABI_FREE(gvnlxc)
738 0 : ABI_FREE(direc)
739 0 : ABI_FREE(pcon)
740 0 : ABI_FREE(scprod)
741 0 : ABI_FREE(scprod_csc)
742 :
743 0 : ABI_FREE(direc_tmp)
744 0 : ABI_FREE(cwavef_r)
745 0 : ABI_FREE(direc_r)
746 :
747 : ! Do not delete this line, needed to run with open MP
748 0 : write(unit=message,fmt=*) resid(1)
749 :
750 0 : call timab(1300,2,tsec)
751 :
752 : DBG_EXIT("COLL")
753 :
754 0 : end subroutine cgwf_cprj
755 : !!***
756 :
757 : !!****f* m_cgwf/mksubovl
758 : !! NAME
759 : !! mksubovl
760 : !!
761 : !! FUNCTION
762 : !! Compute the triangular matrix <c_m|S|c_n> for m<=n
763 : !!
764 : !! INPUTS
765 : !! cg(2,mcg)=wavefunction <G|C band,k> coefficients for ALL bands
766 : !! cprj_cwavef_bands=<p_i|c_n> coefficients for all bands n
767 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
768 : !! icg=shift to be applied on the location of data in the array cg
769 : !! nband=number of bands.
770 : !! mpi_enreg=information about MPI parallelization
771 : !!
772 : !! OUTPUT
773 : !! subovl(nband*(nband+1)) = <c_m|S|c_n> for m<=n
774 : !!
775 : !! SIDE EFFECTS
776 : !!
777 : !! NOTES
778 : !!
779 : !! SOURCE
780 : !!
781 0 : subroutine mksubovl(cg,cprj_cwavef_bands,gs_hamk,icg,nband,subovl,mpi_enreg)
782 : !Arguments ------------------------------------
783 : integer,intent(in) :: icg,nband
784 : !arrays
785 : real(dp),intent(out) :: subovl(nband*(nband+1))
786 : type(pawcprj_type),intent(in),target :: cprj_cwavef_bands(:,:)
787 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
788 : type(MPI_type),intent(in) :: mpi_enreg
789 : real(dp),intent(in),target :: cg(:,:)
790 :
791 : !Local variables-------------------------------
792 : integer,parameter :: tim_getcsc=4
793 : integer :: cpopt,iband,isubh,nspinor,wfsize
794 0 : real(dp),pointer :: cwavef(:,:),cwavef_left(:,:)
795 0 : real(dp),pointer :: cwavef_bands(:,:)
796 0 : type(pawcprj_type),pointer :: cprj_cwavef(:,:),cprj_cwavef_left(:,:)
797 :
798 : ! *********************************************************************
799 :
800 0 : nspinor = gs_hamk%nspinor
801 0 : wfsize=gs_hamk%npw_k*nspinor
802 0 : cpopt=2
803 :
804 0 : cwavef_bands => cg(:,1+icg:nband*wfsize+icg)
805 :
806 0 : isubh=1
807 0 : do iband=1,nband
808 0 : cwavef => cwavef_bands(:,1+(iband-1)*wfsize:iband*wfsize)
809 0 : cprj_cwavef => cprj_cwavef_bands(:,nspinor*(iband-1)+1:nspinor*iband)
810 0 : cwavef_left => cwavef_bands(:,1:iband*wfsize)
811 0 : cprj_cwavef_left => cprj_cwavef_bands(:,1:nspinor*iband)
812 : ! Compute csc matrix
813 : call getcsc(subovl(isubh:isubh+2*iband-1),cpopt,cwavef,cwavef_left,&
814 : & cprj_cwavef,cprj_cwavef_left,&
815 0 : & gs_hamk,mpi_enreg,iband,tim_getcsc)
816 0 : isubh=isubh+2*iband
817 : end do
818 :
819 0 : end subroutine mksubovl
820 : !!***
821 :
822 : !!****f* m_cgwf/cprj_update
823 : !! NAME
824 : !! cprj_update
825 : !!
826 : !! FUNCTION
827 : !! Compute the cprj = <p_i|c_n> from input cg (at fixed k point).
828 : !!
829 : !! INPUTS
830 : !! cg(2,mcg)=wavefunction <G|C band,k> coefficients for ALL bands
831 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
832 : !! icg=shift to be applied on the location of data in the array cg
833 : !! nband=number of bands.
834 : !! mpi_enreg=information about MPI parallelization
835 : !!
836 : !! OUTPUT
837 : !! cprj_cwavef_bands=<p_i|c_n> coefficients for all bands n
838 : !!
839 : !! SIDE EFFECTS
840 : !!
841 : !! NOTES
842 : !!
843 : !! SOURCE
844 : !!
845 0 : subroutine cprj_update(cg,cprj_cwavef_bands,gs_hamk,icg,nband,mpi_enreg,tim_getcprj)
846 :
847 : !Arguments ------------------------------------
848 : integer,intent(in) :: icg,nband,tim_getcprj
849 : !arrays
850 : type(pawcprj_type),intent(in),target :: cprj_cwavef_bands(:,:)
851 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
852 : type(MPI_type),intent(in) :: mpi_enreg
853 : real(dp),intent(inout),target :: cg(:,:)
854 :
855 : !Local variables-------------------------------
856 : integer :: choice,iband,wfsize
857 : real(dp) :: tsec(2)
858 0 : real(dp),pointer :: cwavef(:,:),cwavef_bands(:,:)
859 0 : type(pawcprj_type),pointer :: cprj_cwavef(:,:)
860 :
861 0 : wfsize=gs_hamk%npw_k*gs_hamk%nspinor
862 0 : cwavef_bands => cg(:,1+icg:nband*wfsize+icg)
863 :
864 0 : choice = 1
865 0 : if (cprj_cwavef_bands(1,1)%ncpgr==3) then
866 0 : choice = 2
867 : end if
868 :
869 0 : do iband=1,nband
870 0 : cwavef => cwavef_bands(:,1+(iband-1)*wfsize:iband*wfsize)
871 0 : cprj_cwavef => cprj_cwavef_bands(:,gs_hamk%nspinor*(iband-1)+1:gs_hamk%nspinor*iband)
872 :
873 0 : call timab(1290+tim_getcprj,1,tsec)
874 : call getcprj(choice,0,cwavef,cprj_cwavef,&
875 : & gs_hamk%ffnl_k,0,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,gs_hamk%kpg_k,gs_hamk%kpt_k,&
876 : & gs_hamk%lmnmax,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%natom,gs_hamk%nattyp,&
877 : & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,&
878 0 : & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
879 0 : call timab(1290+tim_getcprj,2,tsec)
880 : end do
881 :
882 0 : end subroutine cprj_update
883 : !!***
884 :
885 : !!****f* m_cgwf/cprj_update_oneband
886 : !! NAME
887 : !! cprj_update_oneband
888 : !!
889 : !! FUNCTION
890 : !! Compute the cprj = <p_i|c_n> from input wavefunction (one band only).
891 : !!
892 : !! INPUTS
893 : !! cwavef(2,npw)=a wavefunction <G|C band,k>
894 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
895 : !! icg=shift to be applied on the location of data in the array cg
896 : !! mpi_enreg=information about MPI parallelization
897 : !!
898 : !! OUTPUT
899 : !! cprj_cwavef=<p_i|c_n> coefficients of the input WF
900 : !!
901 : !! SIDE EFFECTS
902 : !!
903 : !! NOTES
904 : !!
905 : !! SOURCE
906 : !!
907 0 : subroutine cprj_update_oneband(cwavef,cprj_cwavef,gs_hamk,mpi_enreg,tim_getcprj)
908 :
909 : !Arguments ------------------------------------
910 : !arrays
911 : integer,intent(in) :: tim_getcprj
912 : real(dp),intent(inout) :: cwavef(:,:)
913 : type(pawcprj_type),intent(inout) :: cprj_cwavef(:,:)
914 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
915 : type(MPI_type),intent(in) :: mpi_enreg
916 :
917 : !Local variables-------------------------------
918 : integer :: choice,wfsize
919 : real(dp) :: tsec(2)
920 :
921 0 : wfsize=gs_hamk%npw_k*gs_hamk%nspinor
922 :
923 0 : choice = 1
924 0 : if (cprj_cwavef(1,1)%ncpgr==3) then
925 0 : choice = 2
926 : end if
927 :
928 0 : call timab(1290+tim_getcprj,1,tsec)
929 : call getcprj(choice,0,cwavef,cprj_cwavef,&
930 : & gs_hamk%ffnl_k,0,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,gs_hamk%kpg_k,gs_hamk%kpt_k,&
931 : & gs_hamk%lmnmax,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%natom,gs_hamk%nattyp,&
932 : & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,&
933 0 : & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
934 0 : call timab(1290+tim_getcprj,2,tsec)
935 :
936 0 : end subroutine cprj_update_oneband
937 : !!***
938 :
939 : !!****f* m_cgwf/cprj_check
940 : !! NAME
941 : !! cprj_check
942 : !!
943 : !! FUNCTION
944 : !! Check if the cprj are consistent with the bands contained in cg.
945 : !! Useful for debugging only.
946 : !!
947 : !! INPUTS
948 : !! cwavef(2,npw)=a wavefunction <G|C band,k>
949 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
950 : !! icg=shift to be applied on the location of data in the array cg
951 : !! message=string to specify the context of the test
952 : !! mpi_enreg=information about MPI parallelization
953 : !! nband=number of bands.
954 : !!
955 : !! OUTPUT
956 : !!
957 : !! SIDE EFFECTS
958 : !!
959 : !! NOTES
960 : !!
961 : !! SOURCE
962 : !!
963 0 : subroutine cprj_check(cg,cprj_cwavef_bands,gs_hamk,icg,nband,message,mpi_enreg)
964 :
965 : !Arguments ------------------------------------
966 : integer,intent(in) :: icg,nband
967 : !arrays
968 : type(pawcprj_type),intent(in),target :: cprj_cwavef_bands(:,:)
969 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
970 : type(MPI_type),intent(in) :: mpi_enreg
971 : real(dp),intent(inout),target :: cg(:,:)
972 : character(len=*),intent(in) :: message
973 :
974 : !Local variables-------------------------------
975 : integer :: choice,iband,ispinor,ncpgr,wfsize
976 0 : real(dp),pointer :: cwavef(:,:),cwavef_bands(:,:)
977 : integer :: iatom
978 : real(dp) :: re,ratio
979 0 : type(pawcprj_type),allocatable :: cprj_tmp(:,:)
980 :
981 0 : write(std_out,'(a)') ''
982 0 : write(std_out,'(2a)') 'cprj_check : ',message
983 :
984 0 : choice = 1
985 0 : ncpgr = 0
986 0 : if (cprj_cwavef_bands(1,1)%ncpgr==3) then
987 0 : choice = 2
988 0 : ncpgr = 3
989 : end if
990 0 : write(std_out,'(a,i3)') 'ncpgr : ',ncpgr
991 :
992 0 : ABI_MALLOC(cprj_tmp,(gs_hamk%natom,gs_hamk%nspinor))
993 0 : call pawcprj_alloc(cprj_tmp,ncpgr,gs_hamk%dimcprj)
994 :
995 0 : wfsize=gs_hamk%npw_k*gs_hamk%nspinor
996 0 : cwavef_bands => cg(:,1+icg:nband*wfsize+icg)
997 :
998 0 : do iband=1,nband
999 0 : cwavef => cwavef_bands(:,1+(iband-1)*wfsize:iband*wfsize)
1000 : call getcprj(choice,0,cwavef,cprj_tmp,&
1001 : & gs_hamk%ffnl_k,0,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,gs_hamk%kpg_k,gs_hamk%kpt_k,&
1002 : & gs_hamk%lmnmax,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%natom,gs_hamk%nattyp,&
1003 : & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,&
1004 0 : & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
1005 0 : do ispinor=1,gs_hamk%nspinor
1006 0 : do iatom=1,gs_hamk%natom
1007 0 : re = sum(abs(cprj_tmp(iatom,ispinor)%cp-cprj_cwavef_bands(iatom,gs_hamk%nspinor*(iband-1)+ispinor)%cp))
1008 0 : if (re>tol6) then
1009 0 : ratio = sum(abs(cprj_tmp(iatom,ispinor)%cp))/sum(abs(cprj_cwavef_bands(iatom,gs_hamk%nspinor*(iband-1)+ispinor)%cp))
1010 0 : write(std_out,'(a)') 'cprj_check:'
1011 0 : write(std_out,'(a,2i5,2es11.3e3)') 'iband,iatom:',iband,iatom,re,ratio
1012 0 : write(std_out,'(a)') ''
1013 0 : flush(std_out)
1014 0 : ABI_ERROR('dif too large')
1015 : end if
1016 0 : if (ncpgr>0) then
1017 0 : re = sum(abs(cprj_tmp(iatom,ispinor)%dcp-cprj_cwavef_bands(iatom,gs_hamk%nspinor*(iband-1)+ispinor)%dcp))
1018 0 : if (re>tol6) then
1019 0 : write(std_out,'(a)') 'cprj_check:'
1020 0 : write(std_out,'(a,2i5,es11.3e3)') 'iband,iatom:',iband,iatom,re
1021 0 : write(std_out,'(a)') ''
1022 0 : flush(std_out)
1023 0 : ABI_ERROR('dif too large (dcp)')
1024 : end if
1025 : end if
1026 : end do
1027 : end do
1028 : end do
1029 :
1030 0 : write(std_out,'(a)') 'cprj_check : ok'
1031 0 : flush(std_out)
1032 :
1033 0 : call pawcprj_free(cprj_tmp)
1034 0 : ABI_FREE(cprj_tmp)
1035 :
1036 0 : end subroutine cprj_check
1037 : !!***
1038 :
1039 : !!****f* m_cgwf/cprj_check_oneband
1040 : !! NAME
1041 : !! cprj_check_oneband
1042 : !!
1043 : !! FUNCTION
1044 : !! Check if the input cprj is consistent with the input WF.
1045 : !! Useful for debugging only.
1046 : !!
1047 : !! INPUTS
1048 : !! cwavef(2,npw)=a wavefunction <G|C band,k>
1049 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
1050 : !! icg=shift to be applied on the location of data in the array cg
1051 : !! message=string to specify the context of the test
1052 : !! mpi_enreg=information about MPI parallelization
1053 : !!
1054 : !! OUTPUT
1055 : !!
1056 : !! SIDE EFFECTS
1057 : !!
1058 : !! NOTES
1059 : !!
1060 : !! SOURCE
1061 : !!
1062 : subroutine cprj_check_oneband(cwavef,cprj_cwavef,gs_hamk,message,mpi_enreg)
1063 :
1064 : !Arguments ------------------------------------
1065 : !arrays
1066 : type(pawcprj_type),intent(in),target :: cprj_cwavef(:,:)
1067 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1068 : type(MPI_type),intent(in) :: mpi_enreg
1069 : real(dp),intent(inout) :: cwavef(:,:)
1070 : character(len=*),intent(in) :: message
1071 :
1072 : !Local variables-------------------------------
1073 : integer :: choice,ncpgr,wfsize
1074 : integer :: iatom
1075 : real(dp) :: re
1076 : type(pawcprj_type),allocatable :: cprj_tmp(:,:)
1077 :
1078 : write(std_out,'(a)') ''
1079 : write(std_out,'(2a)') 'cprj_check (oneband): ',message
1080 :
1081 : choice = 1
1082 : ncpgr = 0
1083 : if (cprj_cwavef(1,1)%ncpgr==3) then
1084 : choice = 2
1085 : ncpgr = 3
1086 : end if
1087 : write(std_out,'(a,i3)') 'ncpgr : ',ncpgr
1088 :
1089 : ABI_MALLOC(cprj_tmp,(gs_hamk%natom,1))
1090 : call pawcprj_alloc(cprj_tmp,ncpgr,gs_hamk%dimcprj)
1091 :
1092 : wfsize=gs_hamk%npw_k*gs_hamk%nspinor
1093 :
1094 : call getcprj(choice,0,cwavef,cprj_tmp,&
1095 : & gs_hamk%ffnl_k,0,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,gs_hamk%kpg_k,gs_hamk%kpt_k,&
1096 : & gs_hamk%lmnmax,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%natom,gs_hamk%nattyp,&
1097 : & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,&
1098 : & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
1099 : do iatom=1,gs_hamk%natom
1100 : re = sum(abs(cprj_tmp(iatom,1)%cp-cprj_cwavef(iatom,1)%cp))
1101 : if (re>tol6) then
1102 : write(std_out,'(a)') 'cprj_check (oneband):'
1103 : write(std_out,'(a,i5,es11.3e3)') 'iatom:',iatom,re
1104 : write(std_out,'(a)') ''
1105 : flush(std_out)
1106 : ABI_ERROR('dif too large')
1107 : end if
1108 : if (ncpgr>0) then
1109 : re = sum(abs(cprj_tmp(iatom,1)%dcp-cprj_cwavef(iatom,1)%dcp))
1110 : if (re>tol6) then
1111 : write(std_out,'(a)') 'cprj_check (oneband):'
1112 : write(std_out,'(a,i5,es11.3e3)') 'iatom:',iatom,re
1113 : write(std_out,'(a)') ''
1114 : flush(std_out)
1115 : ABI_ERROR('dif too large (dcp)')
1116 : end if
1117 : end if
1118 : end do
1119 :
1120 : write(std_out,'(a)') 'cprj_check : ok'
1121 : flush(std_out)
1122 :
1123 : call pawcprj_free(cprj_tmp)
1124 : ABI_FREE(cprj_tmp)
1125 :
1126 : end subroutine cprj_check_oneband
1127 : !!***
1128 :
1129 : !!****f* m_cgwf/get_cprj_id
1130 : !! NAME
1131 : !! get_cprj_id
1132 : !!
1133 : !! FUNCTION
1134 : !! Get an id from cprj array.
1135 : !! Useful for debugging only.
1136 : !!
1137 : !! INPUTS
1138 : !!
1139 : !! OUTPUT
1140 : !!
1141 : !! SIDE EFFECTS
1142 : !!
1143 : !! NOTES
1144 : !!
1145 : !! SOURCE
1146 : !!
1147 0 : real(dp) function get_cprj_id(cprj)
1148 :
1149 : !Arguments ------------------------------------
1150 : !arrays
1151 : type(pawcprj_type),intent(in) :: cprj(:,:)
1152 :
1153 : !Local variables-------------------------------
1154 : integer :: i1,i2
1155 : real(dp) :: id_tmp
1156 :
1157 0 : id_tmp=zero
1158 0 : do i2=1,size(cprj,2)
1159 0 : do i1=1,size(cprj,1)
1160 0 : id_tmp = id_tmp + sum(abs(cprj(i1,i2)%cp))
1161 : end do
1162 : end do
1163 0 : get_cprj_id=id_tmp
1164 :
1165 0 : end function get_cprj_id
1166 : !!***
1167 :
1168 : !!****f* m_cgwf/get_cwavefr
1169 : !! NAME
1170 : !! get_cwavefr
1171 : !!
1172 : !! FUNCTION
1173 : !! Compute the real space wavefunction by FFT of the input WF in reciprocal space.
1174 : !!
1175 : !! INPUTS
1176 : !! cwavef(2,npw)=a wavefunction <G|C band,k>
1177 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
1178 : !! mpi_enreg=information about MPI parallelization
1179 : !!
1180 : !! OUTPUT
1181 : !! cwavef_r(2,n4,n5,n6,nspinor) = real space coefficients of the WF
1182 : !!
1183 : !! SIDE EFFECTS
1184 : !!
1185 : !! NOTES
1186 : !!
1187 : !! SOURCE
1188 : !!
1189 0 : subroutine get_cwavefr(cwavef,cwavef_r,gs_hamk,mpi_enreg)
1190 :
1191 : !Arguments ------------------------------------
1192 : !arrays
1193 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1194 : type(MPI_type),intent(in) :: mpi_enreg
1195 : real(dp),intent(in),target :: cwavef(:,:)
1196 : real(dp),intent(out) :: cwavef_r(:,:,:,:,:)
1197 :
1198 : !Local variables-------------------------------
1199 : integer,parameter :: tim_fourwf=40
1200 : integer :: n4,n5,n6,npw
1201 : real(dp),parameter :: weight_fft = one
1202 0 : real(dp), pointer :: cwavef_fft1(:,:),cwavef_fft2(:,:)
1203 0 : real(dp),allocatable :: denpot_dum(:,:,:),fofgout_dum(:,:)
1204 :
1205 0 : n4=gs_hamk%ngfft(4);n5=gs_hamk%ngfft(5);n6=gs_hamk%ngfft(6)
1206 0 : npw=gs_hamk%npw_k
1207 :
1208 0 : ABI_MALLOC(denpot_dum, (0,0,0))
1209 0 : ABI_MALLOC(fofgout_dum, (0,0))
1210 :
1211 0 : if (gs_hamk%nspinor==1) then
1212 0 : cwavef_fft1 => cwavef
1213 : else
1214 0 : cwavef_fft1 => cwavef(:,1:npw)
1215 0 : cwavef_fft2 => cwavef(:,1+npw:2*npw)
1216 : end if
1217 : call fourwf(0,denpot_dum,cwavef_fft1,fofgout_dum,cwavef_r(:,:,:,:,1),gs_hamk%gbound_k,gs_hamk%gbound_k,gs_hamk%istwf_k,&
1218 : & gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,gs_hamk%npw_fft_k,gs_hamk%npw_fft_k,&
1219 0 : & n4,n5,n6,0,tim_fourwf,weight_fft,weight_fft)
1220 0 : if (gs_hamk%nspinor==2) then
1221 : call fourwf(0,denpot_dum,cwavef_fft2,fofgout_dum,cwavef_r(:,:,:,:,2),gs_hamk%gbound_k,gs_hamk%gbound_k,gs_hamk%istwf_k,&
1222 : & gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,gs_hamk%npw_fft_k,gs_hamk%npw_fft_k,&
1223 0 : & n4,n5,n6,0,tim_fourwf,weight_fft,weight_fft)
1224 : end if
1225 :
1226 0 : ABI_FREE(denpot_dum)
1227 0 : ABI_FREE(fofgout_dum)
1228 :
1229 0 : end subroutine get_cwavefr
1230 : !!***
1231 :
1232 : end module m_cgwf_cprj
1233 : !!***
|
