Branch data 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 : : !!***
|
