TABLE OF CONTENTS
- ABINIT/m_cgwf_cprj
- m_cgwf/cgwf_cprj
- m_cgwf/cprj_check
- m_cgwf/cprj_check_oneband
- m_cgwf/cprj_update
- m_cgwf/cprj_update_oneband
- m_cgwf/get_cprj_id
- m_cgwf/get_cwavefr
- m_cgwf/mksubovl
ABINIT/m_cgwf_cprj [ Modules ]
NAME
m_cgwf_cprj
FUNCTION
Conjugate-gradient eigensolver.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (DCA, XG, GMR, MT, MVeithen, ISouza, JIniguez) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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
m_cgwf/cgwf_cprj [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
cgwf_cprj
FUNCTION
Update all wavefunction |C>, non self-consistently. Uses a conjugate-gradient algorithm. This version is available for PAW only, as it needs the cprj in memory.
INPUTS
gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg mcg=second dimension of the cg array mpi_enreg=information about MPI parallelization nband=number of bands. nline=number of line minimizations per band. ortalg=governs the choice of the algorithm for orthogonalisation. prtvol=control print volume and debugging output tolrde=tolerance on the ratio of differences of energies (for the line minimisation) tolwfr=tolerance on largest wf residual wfoptalg=govern the choice of algorithm for wf optimisation (0, 1, 10 and 11 : in the present routine, usual CG algorithm ;
OUTPUT
resid(nband)=wf residual for new states=|(H-e)|C>|^2 (hartree^2)
SIDE EFFECTS
cg(2,mcg)
at input =wavefunction <G|C band,k> coefficients for ALL bands
at output same as input except that
the current band, with number 'band' has been updated
cprj_cwavef_bands=<p_i|c_n> coefficients for all bands n, they are updated when the WFs change.
quit= if 1, proceeds to smooth ending of the job.
NOTES
SOURCE
99 subroutine cgwf_cprj(cg,cprj_cwavef_bands,cprj_update_lvl,eig,& 100 & 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 real(dp),allocatable :: conjgr(:,:),gvnlxc(:,:) 130 real(dp), pointer :: cwavef(:,:),cwavef_left(:,:),cwavef_bands(:,:) 131 real(dp), allocatable :: cwavef_r(:,:,:,:,:) 132 real(dp),allocatable,target :: direc(:,:) 133 real(dp),allocatable :: direc_tmp(:,:),pcon(:),scprod(:,:),scwavef_dum(:,:) 134 real(dp),allocatable :: direc_r(:,:,:,:,:),scprod_csc(:) 135 real(dp) :: z_tmp(2),z_tmp2(2) 136 type(pawcprj_type),pointer :: cprj_cwavef(:,:),cprj_cwavef_left(:,:) 137 type(pawcprj_type),allocatable :: cprj_direc(:,:),cprj_conjgr(:,:) 138 139 ! ********************************************************************* 140 141 DBG_ENTER("COLL") 142 143 !Starting the routine 144 call timab(1300,1,tsec) 145 146 !Some checks 147 !Only wfoptalg==10 for now 148 if (wfoptalg/=10) then 149 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 me_g0 = mpi_enreg%me_g0 158 159 ! cprj are already in memory 160 cpopt = 2 161 162 !Initializations and allocations 163 istwf_k=gs_hamk%istwf_k 164 wfopta10=mod(wfoptalg,10) 165 166 optekin=0;if (wfoptalg>=10) optekin=1 167 natom=gs_hamk%natom 168 npw=gs_hamk%npw_k 169 nspinor=gs_hamk%nspinor 170 171 ABI_MALLOC(pcon,(npw)) 172 ABI_MALLOC(conjgr,(2,npw*nspinor)) 173 ABI_MALLOC(scwavef_dum,(0,0)) 174 ABI_MALLOC(direc,(2,npw*nspinor)) 175 ABI_MALLOC(scprod,(2,nband)) 176 ABI_MALLOC(scprod_csc,(2*nband)) 177 ABI_MALLOC(direc_tmp,(2,npw*nspinor)) 178 ABI_MALLOC(gvnlxc,(0,0)) 179 180 ABI_MALLOC(cprj_direc ,(natom,nspinor)) 181 ABI_MALLOC(cprj_conjgr ,(natom,nspinor)) 182 183 n4=gs_hamk%ngfft(4);n5=gs_hamk%ngfft(5);n6=gs_hamk%ngfft(6) 184 ABI_MALLOC(cwavef_r,(2,n4,n5,n6,nspinor)) 185 ABI_MALLOC(direc_r, (2,n4,n5,n6,nspinor)) 186 187 ncpgr = 0 ! no need of gradients here... 188 call pawcprj_alloc(cprj_direc,ncpgr,gs_hamk%dimcprj) 189 call pawcprj_alloc(cprj_conjgr,ncpgr,gs_hamk%dimcprj) 190 191 cwavef_bands => cg(:,1+icg:nband*npw*nspinor+icg) 192 193 if (cprj_update_lvl==-2)then 194 call cprj_update(cg,cprj_cwavef_bands,gs_hamk,icg,nband,mpi_enreg,tim_getcprj) 195 endif 196 197 isubh=1 198 isubh0=1 199 200 ! Big iband loop 201 do iband=1,nband 202 203 ! ====================================================================== 204 ! ========== INITIALISATION OF MINIMIZATION ITERATIONS ================= 205 ! ====================================================================== 206 207 if (prtvol>=10) then ! Tell us what is going on: 208 write(message, '(a,i6,2x,a,i3,a)' )' --- cgwf is called for band',iband,'for',nline,' lines' 209 call wrtout(std_out,message,'PERS') 210 end if 211 212 dotgp=one 213 214 ! Extraction of the vector that is iteratively updated 215 cwavef => cwavef_bands(:,1+(iband-1)*npw*nspinor:iband*npw*nspinor) 216 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 & gs_hamk,mpi_enreg,1,tim_getcsc) 221 xnorm=one/sqrt(dot(1)) 222 z_tmp = (/xnorm,zero/) 223 ! cwavef = xnorm * cwavef 224 call timab(1305,1,tsec) 225 call cg_zscal(npw*nspinor,z_tmp,cwavef) 226 call timab(1305,2,tsec) 227 ! cprj = xnorm * cprj 228 call timab(1302,1,tsec) 229 call pawcprj_axpby(zero,xnorm,cprj_cwavef,cprj_cwavef) 230 call timab(1302,2,tsec) 231 232 ! Compute wavefunction in real space 233 call get_cwavefr(cwavef,cwavef_r,gs_hamk,mpi_enreg) 234 235 if (prtvol==-level) then 236 write(message,'(a,f26.14)')' cgwf: xnorm = ',xnorm 237 call wrtout(std_out,message,'PERS') 238 end if 239 240 ! ====================================================================== 241 ! ====== BEGIN LOOP FOR A GIVEN BAND: MINIMIZATION ITERATIONS ========== 242 ! ====================================================================== 243 if(nline/=0)then 244 do iline=1,nline 245 246 ! === COMPUTE THE RESIDUAL === 247 ! Compute lambda = <C|H|C> 248 sij_opt=0 249 call getchc(z_tmp,cpopt,cwavef,cwavef,cprj_cwavef,cprj_cwavef,cwavef_r,cwavef_r,& 250 & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc) 251 chc=z_tmp(1) 252 lam0=chc 253 eval=chc 254 eig(iband)=chc 255 256 ! Check that lam0 is decreasing on succeeding lines: 257 if (iline==1) then 258 lamold=lam0 259 else 260 if (lam0 > lamold+tol12) then 261 write(message, '(a,i8,a,1p,e14.6,a1,3x,a,1p,e14.6,a1)')& 262 & 'New trial energy at line ',iline,' = ',lam0,ch10,& 263 & 'is higher than former =',lamold,ch10 264 ABI_WARNING(message) 265 end if 266 lamold=lam0 267 end if 268 269 ! ====================================================================== 270 ! =========== COMPUTE THE STEEPEST DESCENT DIRECTION =================== 271 ! ====================================================================== 272 273 sij_opt=-1 274 call getghc(cpopt,cwavef,cprj_cwavef,direc,scwavef_dum,gs_hamk,gvnlxc,& 275 & eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,type_calc,cwavef_r=cwavef_r) 276 277 ! Compute residual (squared) norm 278 call timab(1305,1,tsec) 279 call sqnorm_g(resid(iband),istwf_k,npw*nspinor,direc,me_g0,mpi_enreg%comm_fft) 280 call timab(1305,2,tsec) 281 282 if (prtvol==-level) then 283 write(message,'(a,i0,f26.14,es27.14e3)')' cgwf: iline,eval,resid = ',iline,eval,resid(iband) 284 call wrtout(std_out,message,'PERS') 285 end if 286 xnormd=1/sqrt(resid(iband)) 287 z_tmp = (/xnormd,zero/) 288 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 if (resid(iband)<tolwfr) then 296 if (prtvol>=10) then 297 write(message, '(a,i4,a,i2,a,es12.4)' ) & 298 & ' cgwf: band ',iband,' converged after ',iline,' line minimizations: resid =',resid(iband) 299 call wrtout(std_out,message,'PERS') 300 end if 301 nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped 302 exit ! Exit from the loop on iline 303 end if 304 305 ! If user require exiting the job, stop line minimisations 306 if (quit==1) then 307 write(message, '(a,i0)' )' cgwf: user require exiting => skip update of band ',iband 308 call wrtout(std_out,message,'PERS') 309 310 nskip=nskip+(nline-iline+1) ! Number of two-way 3D ffts skipped 311 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 if(ortalg>=0)then 324 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 & gs_hamk,mpi_enreg,nband,tim_getcsc) 327 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 & 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 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 call timab(1305,1,tsec) 342 if (wfoptalg<10.or.iline==1) then 343 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 do ispinor=1,nspinor 346 igs=(ispinor-1)*npw 347 !$OMP PARALLEL DO 348 do ipw=1+igs,npw+igs 349 direc(1,ipw)=direc(1,ipw)*pcon(ipw-igs) 350 direc(2,ipw)=direc(2,ipw)*pcon(ipw-igs) 351 end do 352 end do 353 end if 354 call timab(1305,2,tsec) 355 356 ! ======= PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION ============== 357 ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================ 358 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 & 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 if (iline==1) then 363 ! We compensate the normalization of the current band 364 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 scprod = reshape(scprod_csc,(/2,nband/)) 368 call projbd(cg,direc,-1,icg,icg,istwf_k,mcg,mcg,nband,npw,nspinor,& 369 & direc,scprod,1,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft) 370 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 z_tmp = -scprod_csc(2*iband-1:2*iband)*(one-xnorm)/xnorm 388 call timab(1305,1,tsec) 389 do ispinor=1,nspinor 390 igs=(ispinor-1)*npw 391 do ipw=1+igs,npw+igs 392 direc(1,ipw)=direc(1,ipw) + z_tmp(1)*cwavef(1,ipw) - z_tmp(2)*cwavef(2,ipw) 393 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 call timab(1305,2,tsec) 397 end if 398 ! Apply projbd to cprj_direc 399 scprod=-scprod 400 call timab(1303,1,tsec) 401 call pawcprj_projbd(scprod,cprj_cwavef_bands,cprj_direc) 402 call timab(1303,2,tsec) 403 if (iline==1) then 404 ! Same correction than for WFs 405 z_tmp = -scprod_csc(2*iband-1:2*iband)*(one-xnorm)/xnorm 406 z_tmp2 = (/one,zero/) 407 ! cprj = z_tmp*cprjx + z_tmp2*cprj 408 call timab(1302,1,tsec) 409 call pawcprj_zaxpby(z_tmp,z_tmp2,cprj_cwavef,cprj_direc) 410 call timab(1302,2,tsec) 411 end if 412 413 ! ====================================================================== 414 ! ================= COMPUTE THE CONJUGATE-GRADIENT ===================== 415 ! ====================================================================== 416 417 call timab(1305,1,tsec) 418 call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,direc_tmp,me_g0,mpi_enreg%comm_spinorfft) 419 call timab(1305,2,tsec) 420 421 dotgg=dotgg/xnormd**2 422 423 ! MJV: added 5 Feb 2012 - causes divide by 0 on next iteration of iline 424 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 if (iline==1) then 428 gamma=zero 429 dotgp=dotgg 430 call timab(1305,1,tsec) 431 call cg_zcopy(npw*nspinor,direc,conjgr) 432 call timab(1305,2,tsec) 433 call timab(1302,1,tsec) 434 call pawcprj_copy(cprj_direc,cprj_conjgr) 435 call timab(1302,2,tsec) 436 if (prtvol==-level)then 437 write(message,'(a,es27.14e3)')' cgwf: dotgg = ',dotgg 438 call wrtout(std_out,message,'PERS') 439 end if 440 else 441 gamma=dotgg/dotgp 442 dotgp=dotgg 443 444 if (prtvol==-level)then 445 write(message,'(a,2es27.14e3)')' cgwf: dotgg,gamma = ',dotgg,gamma 446 call wrtout(std_out,message,'PERS') 447 end if 448 449 gamma=gamma*xnormd/xnormd_previous 450 451 call timab(1305,1,tsec) 452 !$OMP PARALLEL DO 453 do ipw=1,npw*nspinor 454 conjgr(1,ipw)=direc(1,ipw)+gamma*conjgr(1,ipw) 455 conjgr(2,ipw)=direc(2,ipw)+gamma*conjgr(2,ipw) 456 end do 457 call timab(1305,2,tsec) 458 call timab(1302,1,tsec) 459 call pawcprj_axpby(one,gamma,cprj_direc,cprj_conjgr) 460 call timab(1302,2,tsec) 461 end if 462 call timab(1305,1,tsec) 463 464 xnormd_previous=xnormd 465 466 call getcsc(dot,cpopt,conjgr,conjgr,cprj_conjgr,cprj_conjgr,& 467 & 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 & gs_hamk,mpi_enreg,1,tim_getcsc) 475 dotr=dot(1) 476 doti=dot(2) 477 478 ! Project the conjugated gradient onto the current band 479 call timab(1305,1,tsec) 480 if(istwf_k==1)then 481 482 !$OMP PARALLEL DO 483 do ipw=1,npw*nspinor 484 direc(1,ipw)=conjgr(1,ipw)-(dotr*cwavef(1,ipw)-doti*cwavef(2,ipw)) 485 direc(2,ipw)=conjgr(2,ipw)-(dotr*cwavef(2,ipw)+doti*cwavef(1,ipw)) 486 end do 487 else 488 !$OMP PARALLEL DO 489 do ipw=1,npw*nspinor 490 direc(1,ipw)=conjgr(1,ipw)-dotr*cwavef(1,ipw) 491 direc(2,ipw)=conjgr(2,ipw)-dotr*cwavef(2,ipw) 492 end do 493 end if 494 call timab(1305,2,tsec) 495 call timab(1302,1,tsec) 496 call pawcprj_copy(cprj_conjgr,cprj_direc) 497 z_tmp = (/-dotr,-doti/) 498 z_tmp2 = (/one,zero/) 499 call pawcprj_zaxpby(z_tmp,z_tmp2,cprj_cwavef,cprj_direc) 500 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 & gs_hamk,mpi_enreg,1,tim_getcsc) 509 xnormd=one/sqrt(abs(dot(1))) 510 511 sij_opt=0 512 ! Compute direc in real space 513 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 & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc) 517 dhc=z_tmp(1) 518 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 & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc) 523 dhd=z_tmp(1) 524 dhd=dhd*xnormd**2 525 526 if (prtvol==-level) then 527 write(message,'(a,3f26.14)') 'cgwf: chc,dhc,dhd=',chc,dhc,dhd 528 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 tan2th=2.0_dp*dhc/(chc-dhd) 537 538 if (abs(tan2th)<1.d-05) then 539 costh=1.0_dp-0.125_dp*tan2th**2 540 sinth=0.5_dp*tan2th*(1.0_dp-0.375_dp*tan2th**2) 541 542 ! Check that result is above machine precision 543 if (abs(sinth)<epsilon(0._dp)) then 544 write(message, '(a,es16.4)' ) ' cgwf: converged with tan2th=',tan2th 545 call wrtout(std_out,message,'PERS') 546 ! Number of one-way 3D ffts skipped 547 nskip=nskip+2*(nline-iline) 548 exit ! Exit from the loop on iline 549 end if 550 551 else 552 root=sqrt(1.0_dp+tan2th**2) 553 costh=sqrt(0.5_dp+0.5_dp/root) 554 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 diff=(chc-dhd) 559 ! Swap c and d if value of diff is positive 560 if (diff>zero) then 561 swap=costh 562 costh=-sinth 563 sinth=swap 564 if(prtvol<0 .or. prtvol>=10)then 565 write(message,*)' Note: swap roots, iline,diff=',iline,diff 566 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 sintn=sinth*xnormd 575 576 call timab(1305,1,tsec) 577 !$OMP PARALLEL DO 578 do ipw=1,npw*nspinor 579 cwavef(1,ipw)=cwavef(1,ipw)*costh+direc(1,ipw)*sintn 580 cwavef(2,ipw)=cwavef(2,ipw)*costh+direc(2,ipw)*sintn 581 end do 582 do ispinor=1,nspinor 583 do i3=1,gs_hamk%n6 584 do i2=1,gs_hamk%n5 585 do i1=1,gs_hamk%n4 586 cwavef_r(1,i1,i2,i3,ispinor)=cwavef_r(1,i1,i2,i3,ispinor)*costh+direc_r(1,i1,i2,i3,ispinor)*sintn 587 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 call timab(1305,2,tsec) 593 if (cprj_update_lvl<0) then 594 call cprj_update_oneband(cwavef,cprj_cwavef,gs_hamk,mpi_enreg,tim_getcprj) 595 else 596 call timab(1302,1,tsec) 597 call pawcprj_axpby(sintn,costh,cprj_direc,cprj_cwavef) 598 call timab(1302,2,tsec) 599 end if 600 601 ! ====================================================================== 602 ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY =================== 603 ! ====================================================================== 604 605 ! Compute delta(E) 606 deltae=chc*(costh**2-1._dp)+dhd*sinth**2+2._dp*costh*sinth*dhc 607 608 ! Check convergence and eventually exit 609 if (iline==1) then 610 deold=deltae 611 else if (abs(deltae)<tolrde*abs(deold) .and. iline/=nline .and. wfopta10<2)then 612 if(prtvol>=10)then 613 write(message, '(a,i4,1x,a,1p,e12.4,a,e12.4,a)' ) & 614 & ' cgwf: line',iline,& 615 & ' deltae=',deltae,' < tolrde*',deold,' =>skip lines' 616 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 & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc) 621 eig(iband)=z_tmp(1) 622 nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped 623 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 if (iline==nline) then 628 call getchc(z_tmp,cpopt,cwavef,cwavef,cprj_cwavef,cprj_cwavef,cwavef_r,cwavef_r,& 629 & gs_hamk,zero,mpi_enreg,1,sij_opt,type_calc) 630 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 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 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 sij_opt=0 657 call getghc(cpopt,cwavef,cprj_cwavef,direc,scwavef_dum,gs_hamk,gvnlxc,& 658 & eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,3,cwavef_r=cwavef_r) 659 isubh=isubh0 660 call timab(1304,1,tsec) 661 do jband=1,iband 662 cwavef_left => cwavef_bands(:,1+(jband-1)*npw*nspinor:jband*npw*nspinor) 663 call dotprod_g(subham(isubh),subham(isubh+1),istwf_k,npw*nspinor,2,cwavef_left,direc,me_g0,mpi_enreg%comm_spinorfft) 664 isubh=isubh+2 665 end do 666 call timab(1304,2,tsec) 667 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 & gs_hamk,zero,mpi_enreg,iband,sij_opt,4) 672 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 if(prtvol==-level)then 710 isubh=1 711 do iband=1,nband 712 do jband=1,iband 713 if (jband<=10) then 714 write(message,'(i7,2f26.14)')isubh,subham(isubh:isubh+1) 715 call wrtout(std_out,message,'PERS') 716 end if 717 isubh=isubh+2 718 end do 719 end do 720 end if 721 722 ! =================== 723 ! FINAL DEALLOCATIONS 724 ! =================== 725 726 nullify(cwavef_left) 727 nullify(cwavef_bands) 728 nullify(cwavef) 729 nullify(cprj_cwavef) 730 nullify(cprj_cwavef_left) 731 call pawcprj_free(cprj_direc) 732 call pawcprj_free(cprj_conjgr) 733 ABI_FREE(cprj_direc) 734 ABI_FREE(cprj_conjgr) 735 ABI_FREE(conjgr) 736 ABI_FREE(scwavef_dum) 737 ABI_FREE(gvnlxc) 738 ABI_FREE(direc) 739 ABI_FREE(pcon) 740 ABI_FREE(scprod) 741 ABI_FREE(scprod_csc) 742 743 ABI_FREE(direc_tmp) 744 ABI_FREE(cwavef_r) 745 ABI_FREE(direc_r) 746 747 ! Do not delete this line, needed to run with open MP 748 write(unit=message,fmt=*) resid(1) 749 750 call timab(1300,2,tsec) 751 752 DBG_EXIT("COLL") 753 754 end subroutine cgwf_cprj
m_cgwf/cprj_check [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
cprj_check
FUNCTION
Check if the cprj are consistent with the bands contained in cg. Useful for debugging only.
INPUTS
cwavef(2,npw)=a wavefunction <G|C band,k> gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg message=string to specify the context of the test mpi_enreg=information about MPI parallelization nband=number of bands.
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
963 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 real(dp),pointer :: cwavef(:,:),cwavef_bands(:,:) 977 integer :: iatom 978 real(dp) :: re,ratio 979 type(pawcprj_type),allocatable :: cprj_tmp(:,:) 980 981 write(std_out,'(a)') '' 982 write(std_out,'(2a)') 'cprj_check : ',message 983 984 choice = 1 985 ncpgr = 0 986 if (cprj_cwavef_bands(1,1)%ncpgr==3) then 987 choice = 2 988 ncpgr = 3 989 end if 990 write(std_out,'(a,i3)') 'ncpgr : ',ncpgr 991 992 ABI_MALLOC(cprj_tmp,(gs_hamk%natom,gs_hamk%nspinor)) 993 call pawcprj_alloc(cprj_tmp,ncpgr,gs_hamk%dimcprj) 994 995 wfsize=gs_hamk%npw_k*gs_hamk%nspinor 996 cwavef_bands => cg(:,1+icg:nband*wfsize+icg) 997 998 do iband=1,nband 999 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,gs_hamk%natom,gs_hamk%nattyp,& 1003 & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,& 1004 & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm) 1005 do ispinor=1,gs_hamk%nspinor 1006 do iatom=1,gs_hamk%natom 1007 re = sum(abs(cprj_tmp(iatom,ispinor)%cp-cprj_cwavef_bands(iatom,gs_hamk%nspinor*(iband-1)+ispinor)%cp)) 1008 if (re>tol6) then 1009 ratio = sum(abs(cprj_tmp(iatom,ispinor)%cp))/sum(abs(cprj_cwavef_bands(iatom,gs_hamk%nspinor*(iband-1)+ispinor)%cp)) 1010 write(std_out,'(a)') 'cprj_check:' 1011 write(std_out,'(a,2i5,2es11.3e3)') 'iband,iatom:',iband,iatom,re,ratio 1012 write(std_out,'(a)') '' 1013 flush(std_out) 1014 ABI_ERROR('dif too large') 1015 end if 1016 if (ncpgr>0) then 1017 re = sum(abs(cprj_tmp(iatom,ispinor)%dcp-cprj_cwavef_bands(iatom,gs_hamk%nspinor*(iband-1)+ispinor)%dcp)) 1018 if (re>tol6) then 1019 write(std_out,'(a)') 'cprj_check:' 1020 write(std_out,'(a,2i5,es11.3e3)') 'iband,iatom:',iband,iatom,re 1021 write(std_out,'(a)') '' 1022 flush(std_out) 1023 ABI_ERROR('dif too large (dcp)') 1024 end if 1025 end if 1026 end do 1027 end do 1028 end do 1029 1030 write(std_out,'(a)') 'cprj_check : ok' 1031 flush(std_out) 1032 1033 call pawcprj_free(cprj_tmp) 1034 ABI_FREE(cprj_tmp) 1035 1036 end subroutine cprj_check
m_cgwf/cprj_check_oneband [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
cprj_check_oneband
FUNCTION
Check if the input cprj is consistent with the input WF. Useful for debugging only.
INPUTS
cwavef(2,npw)=a wavefunction <G|C band,k> gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg message=string to specify the context of the test mpi_enreg=information about MPI parallelization
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
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,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
m_cgwf/cprj_update [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
cprj_update
FUNCTION
Compute the cprj = <p_i|c_n> from input cg (at fixed k point).
INPUTS
cg(2,mcg)=wavefunction <G|C band,k> coefficients for ALL bands gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg nband=number of bands. mpi_enreg=information about MPI parallelization
OUTPUT
cprj_cwavef_bands=<p_i|c_n> coefficients for all bands n
SIDE EFFECTS
NOTES
SOURCE
845 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 real(dp),pointer :: cwavef(:,:),cwavef_bands(:,:) 859 type(pawcprj_type),pointer :: cprj_cwavef(:,:) 860 861 wfsize=gs_hamk%npw_k*gs_hamk%nspinor 862 cwavef_bands => cg(:,1+icg:nband*wfsize+icg) 863 864 choice = 1 865 if (cprj_cwavef_bands(1,1)%ncpgr==3) then 866 choice = 2 867 end if 868 869 do iband=1,nband 870 cwavef => cwavef_bands(:,1+(iband-1)*wfsize:iband*wfsize) 871 cprj_cwavef => cprj_cwavef_bands(:,gs_hamk%nspinor*(iband-1)+1:gs_hamk%nspinor*iband) 872 873 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,gs_hamk%natom,gs_hamk%nattyp,& 877 & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,& 878 & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm) 879 call timab(1290+tim_getcprj,2,tsec) 880 end do 881 882 end subroutine cprj_update
m_cgwf/cprj_update_oneband [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
cprj_update_oneband
FUNCTION
Compute the cprj = <p_i|c_n> from input wavefunction (one band only).
INPUTS
cwavef(2,npw)=a wavefunction <G|C band,k> gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg mpi_enreg=information about MPI parallelization
OUTPUT
cprj_cwavef=<p_i|c_n> coefficients of the input WF
SIDE EFFECTS
NOTES
SOURCE
907 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 wfsize=gs_hamk%npw_k*gs_hamk%nspinor 922 923 choice = 1 924 if (cprj_cwavef(1,1)%ncpgr==3) then 925 choice = 2 926 end if 927 928 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,gs_hamk%natom,gs_hamk%nattyp,& 932 & gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,& 933 & gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm) 934 call timab(1290+tim_getcprj,2,tsec) 935 936 end subroutine cprj_update_oneband
m_cgwf/get_cprj_id [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
get_cprj_id
FUNCTION
Get an id from cprj array. Useful for debugging only.
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
1147 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 id_tmp=zero 1158 do i2=1,size(cprj,2) 1159 do i1=1,size(cprj,1) 1160 id_tmp = id_tmp + sum(abs(cprj(i1,i2)%cp)) 1161 end do 1162 end do 1163 get_cprj_id=id_tmp 1164 1165 end function get_cprj_id
m_cgwf/get_cwavefr [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
get_cwavefr
FUNCTION
Compute the real space wavefunction by FFT of the input WF in reciprocal space.
INPUTS
cwavef(2,npw)=a wavefunction <G|C band,k> gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k mpi_enreg=information about MPI parallelization
OUTPUT
cwavef_r(2,n4,n5,n6,nspinor) = real space coefficients of the WF
SIDE EFFECTS
NOTES
SOURCE
1189 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 real(dp), pointer :: cwavef_fft1(:,:),cwavef_fft2(:,:) 1203 real(dp),allocatable :: denpot_dum(:,:,:),fofgout_dum(:,:) 1204 1205 n4=gs_hamk%ngfft(4);n5=gs_hamk%ngfft(5);n6=gs_hamk%ngfft(6) 1206 npw=gs_hamk%npw_k 1207 1208 ABI_MALLOC(denpot_dum, (0,0,0)) 1209 ABI_MALLOC(fofgout_dum, (0,0)) 1210 1211 if (gs_hamk%nspinor==1) then 1212 cwavef_fft1 => cwavef 1213 else 1214 cwavef_fft1 => cwavef(:,1:npw) 1215 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 & n4,n5,n6,0,tim_fourwf,weight_fft,weight_fft) 1220 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 & n4,n5,n6,0,tim_fourwf,weight_fft,weight_fft) 1224 end if 1225 1226 ABI_FREE(denpot_dum) 1227 ABI_FREE(fofgout_dum) 1228 1229 end subroutine get_cwavefr
m_cgwf/mksubovl [ Functions ]
[ Top ] [ m_cgwf ] [ Functions ]
NAME
mksubovl
FUNCTION
Compute the triangular matrix <c_m|S|c_n> for m<=n
INPUTS
cg(2,mcg)=wavefunction <G|C band,k> coefficients for ALL bands cprj_cwavef_bands=<p_i|c_n> coefficients for all bands n gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg nband=number of bands. mpi_enreg=information about MPI parallelization
OUTPUT
subovl(nband*(nband+1)) = <c_m|S|c_n> for m<=n
SIDE EFFECTS
NOTES
SOURCE
781 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 real(dp),pointer :: cwavef(:,:),cwavef_left(:,:) 795 real(dp),pointer :: cwavef_bands(:,:) 796 type(pawcprj_type),pointer :: cprj_cwavef(:,:),cprj_cwavef_left(:,:) 797 798 ! ********************************************************************* 799 800 nspinor = gs_hamk%nspinor 801 wfsize=gs_hamk%npw_k*nspinor 802 cpopt=2 803 804 cwavef_bands => cg(:,1+icg:nband*wfsize+icg) 805 806 isubh=1 807 do iband=1,nband 808 cwavef => cwavef_bands(:,1+(iband-1)*wfsize:iband*wfsize) 809 cprj_cwavef => cprj_cwavef_bands(:,nspinor*(iband-1)+1:nspinor*iband) 810 cwavef_left => cwavef_bands(:,1:iband*wfsize) 811 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 & gs_hamk,mpi_enreg,iband,tim_getcsc) 816 isubh=isubh+2*iband 817 end do 818 819 end subroutine mksubovl