TABLE OF CONTENTS
ABINIT/extrapwf [ Functions ]
NAME
extrapwf
FUNCTION
Extrapolate wavefunctions for new ionic positions from values of wavefunctions of previous SCF cycle. Use algorithm proposed by T. A. Arias et al. in PRB 45, 1538 (1992)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT,FJ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
atindx(natom)=index table for atoms atindx1(natom)=index table for atoms, inverse of atindx dtset <type(dataset_type)>=all input variables in this dataset istep=number of call the routine kg(3,mpw*mkmem)=reduced planewave coordinates. mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization nattyp(ntypat)=number of atoms of each type in cell. ngfft(18)=contain all needed information about 3D FFT npwarr(nkpt)=number of planewaves in basis at this k point ntypat=number of types of atoms in cell pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data psps<type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translation vectors (bohr) usepaw= 0 for non paw calculation; =1 for paw calculation xred_old(3,natom)=old reduced coordinates for atoms in unit cell ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
SIDE EFFECTS
cg(2,mcg)= plane wave wavefunction coefficient Value from previous SCF cycle is input Extrapolated value is output scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
NOTES
THIS ROUTINE IS NOT USEABLE AT PRESENT. SHOULD BE CAREFULY TESTED AND DEBUGGED (ESPECIALLY WITHIN PAW).
PARENTS
extraprho
CHILDREN
ctocprj,dotprod_g,getph,hermit,metric,pawcprj_alloc,pawcprj_copy pawcprj_free,pawcprj_get,pawcprj_getdim,pawcprj_lincom,pawcprj_put pawcprj_zaxpby,xmpi_allgather,xmpi_alltoallv,zgemm,zhpev
SOURCE
59 #if defined HAVE_CONFIG_H 60 #include "config.h" 61 #endif 62 63 #include "abi_common.h" 64 65 subroutine extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,& 66 & nattyp,ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm) 67 68 use defs_basis 69 use defs_abitypes 70 use m_scf_history 71 use m_xmpi 72 use m_profiling_abi 73 use m_errors 74 use m_cgtools 75 76 use m_numeric_tools, only : hermit 77 use m_geometry, only : metric 78 use m_kg, only : getph 79 use defs_datatypes, only : pseudopotential_type 80 use m_pawtab, only : pawtab_type 81 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, & 82 & pawcprj_free, pawcprj_zaxpby, pawcprj_put, pawcprj_getdim 83 84 !This section has been created automatically by the script Abilint (TD). 85 !Do not modify the following lines by hand. 86 #undef ABI_FUNC 87 #define ABI_FUNC 'extrapwf' 88 use interfaces_32_util 89 use interfaces_66_nonlocal 90 !End of the abilint section 91 92 implicit none 93 94 !Arguments ------------------------------------ 95 !scalars 96 integer,intent(in) :: istep,mcg,mgfft,ntypat,usepaw 97 type(MPI_type),intent(in) :: mpi_enreg 98 type(dataset_type),intent(in) :: dtset 99 type(scf_history_type),intent(inout) :: scf_history 100 type(pseudopotential_type),intent(in) :: psps 101 !arrays 102 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfft(18) 103 integer,intent(in) :: npwarr(dtset%nkpt) 104 real(dp),intent(in) :: rprimd(3,3) 105 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 106 real(dp), intent(inout) :: cg(2,mcg) 107 real(dp),intent(in) :: xred_old(3,dtset%natom) 108 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 109 110 !Local variables------------------------------- 111 !scalars 112 integer :: ia,iat,iatom,iband_max,iband_max1,iband_min,iband_min1,ibd,ibg,iblockbd,iblockbd1,icg,icgb,icgb1,icgb2 113 integer :: ierr,ig,ii,ikpt,ilmn1,ilmn2,inc,ind1,ind2,iorder_cprj 114 integer :: isize,isppol,istep1,istwf_k,itypat,klmn,me_distrb,my_nspinor 115 integer :: nband_k,nblockbd,nprocband,npw_k,npw_nk,spaceComm_band 116 real(dp) :: dotr,dotr1,doti,doti1,eigval 117 !character(len=500) :: message 118 !arrays 119 real(dp) :: alpha(2),beta(2),gmet(3,3),gprimd(3,3),rmet(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom),ucvol 120 integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:),dimcprj(:),npw_block(:),npw_disp(:) 121 real(dp),allocatable :: al(:,:),anm(:),cwavef(:,:),cwavef1(:,:),cwavef_tmp(:,:),deltawf1(:,:),deltawf2(:,:) 122 real(dp),allocatable :: eig(:),evec(:,:) 123 real(dp),allocatable :: unm(:,:,:) 124 real(dp),allocatable :: work(:,:),work1(:,:),wf1(:,:),ylmgr_k(:,:,:),zhpev1(:,:),zhpev2(:) 125 complex(dpc),allocatable :: unm_tmp(:,:),anm_tmp(:,:) 126 type(pawcprj_type),allocatable :: cprj(:,:),cprj_k(:,:),cprj_k1(:,:),cprj_k2(:,:),cprj_k3(:,:),cprj_k4(:,:) 127 !complex(dpc) :: aa 128 129 ! ************************************************************************* 130 131 if (istep==0) return 132 133 !Useful array 134 if (usepaw==1) then 135 ABI_ALLOCATE(dimcprj,(dtset%natom)) 136 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O') 137 end if 138 139 !Metric 140 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 141 142 !History indexes 143 ind1=1;ind2=2 144 145 !First step 146 if (istep==1) then 147 scf_history%cg(:,:,ind1)=cg(:,:) 148 ! scf_history%cg(:,:,ind2)=zero 149 scf_history%cg(:,:,ind2)= cg(:,:) 150 if(usepaw==1) then 151 ! WARNING: THIS SECTION IS USELESS; NOW crpj CAN BE READ FROM SCFCV 152 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old) 153 iatom=0 ; iorder_cprj=0 154 call pawcprj_alloc(scf_history%cprj(:,:,ind1),0,dimcprj) 155 call pawcprj_alloc(scf_history%cprj(:,:,ind2),0,dimcprj) 156 ABI_ALLOCATE(ylmgr_k,(dtset%mpw,3,0)) 157 call ctocprj(atindx,cg,1,scf_history%cprj(:,:,ind1),gmet,gprimd,& 158 & iatom,0,iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,& 159 & dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,& 160 & dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,& 161 & dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,dtset%ntypat,& 162 & dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,0,& 163 & xred_old,ylm,ylmgr_k) 164 ABI_DEALLOCATE(ylmgr_k) 165 ! call pawcprj_set_zero(scf_history%cprj(:,:,ind2)) 166 call pawcprj_copy(scf_history%cprj(:,:,ind1),scf_history%cprj(:,:,ind2)) 167 end if 168 else 169 170 !From 2nd step 171 172 ! Init parallelism 173 me_distrb=mpi_enreg%me_kpt 174 if (mpi_enreg%paral_kgb==1.or.mpi_enreg%paralbd==1) then 175 spaceComm_band=mpi_enreg%comm_band 176 nprocband=mpi_enreg%nproc_band 177 else 178 spaceComm_band=xmpi_comm_self 179 nprocband=1 180 end if 181 182 ! For the moment sequential part only 183 nprocband=1 184 185 ! Additional statements if band-fft parallelism 186 if (nprocband>1) then 187 ABI_ALLOCATE(npw_block,(nprocband)) 188 ABI_ALLOCATE(npw_disp,(nprocband)) 189 ABI_ALLOCATE(bufsize,(nprocband)) 190 ABI_ALLOCATE(bufdisp,(nprocband)) 191 ABI_ALLOCATE(bufsize_wf,(nprocband)) 192 ABI_ALLOCATE(bufdisp_wf,(nprocband)) 193 end if 194 195 icg=0 196 ibg=0 197 198 if(usepaw==1) then 199 ! WARNING: THIS SECTION IS USELESS; NOW cprj CAN BE READ FROM SCFCV 200 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old) 201 ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,scf_history%mcprj)) 202 call pawcprj_alloc(cprj,0,dimcprj) 203 iatom=0 ; iorder_cprj=0 204 ABI_ALLOCATE(ylmgr_k,(dtset%mpw,3,0)) 205 call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,iorder_cprj,& 206 & dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,dtset%mgfft,& 207 & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,& 208 & nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,& 209 & npwarr,dtset%nspinor,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,& 210 & ph1d,psps,rmet,dtset%typat,ucvol,0,xred_old,& 211 & ylm,ylmgr_k) 212 ABI_DEALLOCATE(ylmgr_k) 213 end if ! end usepaw=1 214 215 ! LOOP OVER SPINS 216 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 217 do isppol=1,dtset%nsppol 218 219 ! BIG FAT k POINT LOOP 220 do ikpt=1,dtset%nkpt 221 222 ! Select k point to be treated by this proc 223 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 224 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 225 226 istwf_k=dtset%istwfk(ikpt) 227 228 ! Retrieve number of plane waves 229 npw_k=npwarr(ikpt) 230 if (nprocband>1) then 231 ! Special treatment for band-fft // 232 call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr) 233 npw_nk=sum(npw_block);npw_disp(1)=0 234 do ii=2,nprocband 235 npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1) 236 end do 237 else 238 npw_nk=npw_k 239 end if 240 241 ! Allocate arrays for a wave-function (or a block of WFs) 242 ABI_ALLOCATE(cwavef,(2,npw_nk*my_nspinor)) 243 ABI_ALLOCATE(cwavef1,(2,npw_nk*my_nspinor)) 244 if (nprocband>1) then 245 isize=2*my_nspinor;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:) 246 isize=2*my_nspinor*npw_k;bufsize_wf(:)=isize 247 do ii=1,nprocband 248 bufdisp_wf(ii)=(ii-1)*isize 249 end do 250 end if 251 252 ! Subspace alignment 253 254 ! Loop over bands or blocks of bands 255 nblockbd=nband_k/nprocband 256 icgb=icg 257 258 if(usepaw==1) then 259 ABI_DATATYPE_ALLOCATE( cprj_k,(dtset%natom,my_nspinor*nblockbd)) 260 call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimcprj) 261 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,& 262 & dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,& 263 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 264 ABI_DATATYPE_ALLOCATE( cprj_k1,(dtset%natom,my_nspinor*nblockbd)) 265 call pawcprj_alloc(cprj_k1,scf_history%cprj(1,1,ind1)%ncpgr,dimcprj) 266 call pawcprj_get(atindx1,cprj_k1,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,& 267 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,& 268 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 269 ABI_DATATYPE_ALLOCATE( cprj_k2,(dtset%natom,my_nspinor*nblockbd)) 270 call pawcprj_alloc(cprj_k2,scf_history%cprj(1,1,ind2)%ncpgr,dimcprj) 271 call pawcprj_get(atindx1,cprj_k2,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,& 272 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,& 273 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 274 end if !end usepaw=1 275 276 ABI_ALLOCATE(unm,(2,nblockbd,nblockbd)) 277 unm=zero 278 icgb2=0 279 280 do iblockbd=1,nblockbd 281 iband_min=1+(iblockbd-1)*nprocband 282 iband_max=iblockbd*nprocband 283 284 if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then 285 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) cycle 286 end if 287 288 ! Extract wavefunction information 289 if (nprocband>1) then 290 ! Special treatment for band-fft // 291 ABI_ALLOCATE(cwavef_tmp,(2,npw_k*my_nspinor*nprocband)) 292 do ig=1,npw_k*my_nspinor*nprocband 293 cwavef_tmp(1,ig)=cg(1,ig+icgb) 294 cwavef_tmp(2,ig)=cg(2,ig+icgb) 295 end do 296 call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr) 297 ABI_DEALLOCATE(cwavef_tmp) 298 else 299 do ig=1,npw_k*my_nspinor 300 cwavef(1,ig)=cg(1,ig+icgb) 301 cwavef(2,ig)=cg(2,ig+icgb) 302 end do 303 end if 304 305 icgb1=icg 306 307 do iblockbd1=1,nblockbd 308 iband_min1=1+(iblockbd1-1)*nprocband 309 iband_max1=iblockbd1*nprocband 310 311 if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then 312 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min1,iband_max1,isppol,me_distrb)) cycle 313 end if 314 315 ! Extract wavefunction information 316 317 if (nprocband>1) then 318 ! Special treatment for band-fft // 319 ABI_ALLOCATE(cwavef_tmp,(2,npw_k*my_nspinor*nprocband)) 320 do ig=1,npw_k*my_nspinor*nprocband 321 cwavef_tmp(1,ig)=scf_history%cg(1,ig+icgb1,ind1) 322 cwavef_tmp(2,ig)=scf_history%cg(2,ig+icgb1,ind1) 323 end do 324 call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef1,bufsize,bufdisp,spaceComm_band,ierr) 325 ABI_DEALLOCATE(cwavef_tmp) 326 else 327 do ig=1,npw_k*my_nspinor 328 cwavef1(1,ig)=scf_history%cg(1,ig+icgb1,ind1) 329 cwavef1(2,ig)=scf_history%cg(2,ig+icgb1,ind1) 330 end do 331 end if 332 333 ! Calculate Unm=<psi_nk(t)|S|psi_mk(t-dt)> 334 call dotprod_g(dotr,doti,istwf_k,npw_k*my_nspinor,2,cwavef,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 335 if(usepaw==1) then 336 ia =0 337 do itypat=1,ntypat 338 do iat=1+ia,nattyp(itypat)+ia 339 do ilmn1=1,pawtab(itypat)%lmn_size 340 do ilmn2=1,ilmn1 341 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 342 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+& 343 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)) 344 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-& 345 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)) 346 end do 347 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 348 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 349 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+& 350 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)) 351 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-& 352 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)) 353 end do 354 end do 355 end do 356 ia=ia+nattyp(itypat) 357 end do 358 end if 359 ! unm(1,iblockbd,iblockbd1)=dotr 360 ! unm(2,iblockbd,iblockbd1)=doti 361 unm(1,iblockbd1,iblockbd)=dotr 362 unm(2,iblockbd1,iblockbd)=doti 363 ! End loop over bands iblockbd1 364 icgb1=icgb1+npw_k*my_nspinor*nprocband 365 366 end do 367 368 ! End loop over bands iblockbd 369 icgb2=icgb2+npw_k*my_nspinor*nprocband 370 icgb=icgb+npw_k*my_nspinor*nprocband 371 end do 372 373 ! write(std_out,*) 'UNM' 374 ! do iblockbd=1,nblockbd 375 ! write(std_out,11) (unm(1,iblockbd,iblockbd1),unm(2,iblockbd,iblockbd1),iblockbd1=1,nblockbd) 376 ! end do 377 ! 11 format(12(1x,f9.5),a) 378 ! Compute A=tU^*U 379 ABI_ALLOCATE(unm_tmp,(nblockbd,nblockbd)) 380 ABI_ALLOCATE(anm_tmp,(nblockbd,nblockbd)) 381 ABI_ALLOCATE(anm,(nblockbd*(nblockbd+1))) 382 unm_tmp(:,:)=cmplx(unm(1,:,:),unm(2,:,:),kind=dp) 383 call zgemm('C','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp), unm_tmp,nblockbd, & 384 & unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd) 385 do iblockbd=1,nblockbd 386 do iblockbd1=iblockbd,nblockbd 387 ii=iblockbd1*(iblockbd1-1)+2*(iblockbd-1)+1 388 anm(ii)=real(anm_tmp(iblockbd,iblockbd1)) 389 anm(ii+1)=aimag(anm_tmp(iblockbd,iblockbd1)) 390 end do 391 end do 392 call hermit(anm,anm,ierr,nblockbd) 393 ! aa=dcmplx(0._dp) 394 ! do iblockbd=1,nblockbd 395 ! aa=aa+conjg(unm_tmp(iblockbd,1))*unm_tmp(iblockbd,1) 396 ! end do 397 ! write(std_out,*) 'tU*U', aa 398 ! write(std_out,*) 'ANM_tmp' 399 ! do iblockbd=1,nblockbd 400 ! write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd) 401 ! end do 402 ! write(std_out,*) 'ANM' 403 ! do iblockbd=1,nblockbd*(nblockbd+1) 404 ! write(std_out,11) anm(iblockbd) 405 ! end do 406 407 ! Diagonalize A 408 ABI_ALLOCATE(eig,(nblockbd)) 409 ABI_ALLOCATE(evec,(2*nblockbd,nblockbd)) 410 ABI_ALLOCATE(zhpev1,(2,2*nblockbd-1)) 411 ABI_ALLOCATE(zhpev2,(3*nblockbd-2)) 412 call zhpev('V','U',nblockbd,anm,eig,evec,nblockbd,zhpev1,& 413 & zhpev2,ierr) 414 ABI_DEALLOCATE(anm) 415 ABI_DEALLOCATE(zhpev1) 416 ABI_DEALLOCATE(zhpev2) 417 ! aa=dcmplx(0._dp) 418 ! do iblockbd=1,nblockbd 419 ! aa=aa+anm_tmp(1,iblockbd)*cmplx(evec((2*iblockbd-1),1),evec(2*iblockbd,1),kind=dp) 420 ! end do 421 ! write(std_out,*) 'EIG', aa, eig(1)*evec(1,1),eig(1)*evec(2,1) 422 423 ! Compute A'=evec*tU^/sqrt(eig) 424 call zgemm('C','C',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, & 425 & unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd) 426 do iblockbd=1,nblockbd 427 eigval=dsqrt(eig(iblockbd)) 428 do iblockbd1=1,nblockbd 429 anm_tmp(iblockbd,iblockbd1)=anm_tmp(iblockbd,iblockbd1)/eigval 430 end do 431 end do 432 433 ! Compute tA^A'to come back to the initial subspace for the cg's 434 435 call zgemm('N','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, & 436 & anm_tmp,nblockbd,dcmplx(0._dp),unm_tmp,nblockbd) 437 anm_tmp=unm_tmp 438 ! write(std_out,*) 'ANM_tmp' 439 ! do iblockbd=1,nblockbd 440 ! write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd) 441 ! end do 442 443 ! Wavefunction alignment (istwfk=1 ?) 444 ABI_ALLOCATE(work,(2,npw_nk*my_nspinor*nblockbd)) 445 ABI_ALLOCATE(work1,(2,my_nspinor*nblockbd*npw_nk)) 446 work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind1) 447 call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), & 448 & work1,npw_nk*my_nspinor, & 449 & anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor) 450 scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind1)=work(:,:) 451 452 work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind2) 453 call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), & 454 & work1,npw_nk*my_nspinor, & 455 & anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor) 456 scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind2)=work(:,:) 457 ABI_DEALLOCATE(work1) 458 ! If paw, must also align cprj: 459 if (usepaw==1) then 460 ! New version (MT): 461 ABI_DATATYPE_ALLOCATE(cprj_k3,(dtset%natom,my_nspinor)) 462 ABI_DATATYPE_ALLOCATE(cprj_k4,(dtset%natom,my_nspinor)) 463 call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj) 464 call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj) 465 ABI_ALLOCATE(al,(2,nblockbd)) 466 do iblockbd=1,nblockbd 467 ii=(iblockbd-1)*my_nspinor 468 do iblockbd1=1,nblockbd 469 al(1,iblockbd1)=real (anm_tmp(iblockbd,iblockbd1)) 470 al(2,iblockbd1)=aimag(anm_tmp(iblockbd,iblockbd1)) 471 end do 472 call pawcprj_lincom(al,cprj_k1,cprj_k3,nblockbd) 473 call pawcprj_lincom(al,cprj_k2,cprj_k4,nblockbd) 474 call pawcprj_copy(cprj_k3,cprj_k1(:,ii+1:ii+my_nspinor)) 475 call pawcprj_copy(cprj_k4,cprj_k2(:,ii+1:ii+my_nspinor)) 476 end do 477 ABI_DEALLOCATE(al) 478 ! Old version (FJ): 479 ! allocate( cprj_k3(dtset%natom,my_nspinor*nblockbd)) 480 ! call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj) 481 ! allocate( cprj_k4(dtset%natom,my_nspinor*nblockbd)) 482 ! call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj) 483 ! beta(1)=one;beta(2)=zero 484 ! do iblockbd=1,nblockbd*my_nspinor 485 ! do iblockbd1=1,nblockbd*my_nspinor 486 ! alpha(1)=real(anm_tmp(iblockbd,iblockbd1));alpha(2)=aimag(anm_tmp(iblockbd,iblockbd1)) 487 ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd1:iblockbd1),cprj_k3(:,iblockbd:iblockbd)) 488 ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd1:iblockbd1),cprj_k4(:,iblockbd:iblockbd)) 489 ! end do 490 ! end do 491 ! call pawcprj_copy(cprj_k3,cprj_k1) 492 ! call pawcprj_copy(cprj_k4,cprj_k2) 493 494 call pawcprj_free(cprj_k3) 495 call pawcprj_free(cprj_k4) 496 ABI_DATATYPE_DEALLOCATE(cprj_k3) 497 ABI_DATATYPE_DEALLOCATE(cprj_k4) 498 end if 499 ABI_DEALLOCATE(anm_tmp) 500 ABI_DEALLOCATE(unm_tmp) 501 ABI_DEALLOCATE(work) 502 503 ! Wavefunction extrapolation 504 ibd=0 505 inc=npw_nk*my_nspinor 506 ABI_ALLOCATE(deltawf2,(2,npw_nk*my_nspinor)) 507 ABI_ALLOCATE(wf1,(2,npw_nk*my_nspinor)) 508 ABI_ALLOCATE(deltawf1,(2,npw_nk*my_nspinor)) 509 do iblockbd=1,nblockbd 510 deltawf2(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2) 511 wf1(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1) 512 ! wf1(2,1)=zero;deltawf2(2,1)=zero 513 514 call dotprod_g(dotr,doti,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),cg(:,icg+1+ibd:ibd+icg+inc),& 515 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 516 call dotprod_g(dotr1,doti1,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),wf1,& 517 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 518 if(usepaw==1) then 519 ia =0 520 do itypat=1,ntypat 521 do iat=1+ia,nattyp(itypat)+ia 522 do ilmn1=1,pawtab(itypat)%lmn_size 523 do ilmn2=1,ilmn1 524 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 525 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+& 526 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)) 527 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-& 528 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)) 529 dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+& 530 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)) 531 doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-& 532 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)) 533 end do 534 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 535 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 536 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+& 537 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)) 538 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-& 539 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)) 540 dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+& 541 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)) 542 doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-& 543 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)) 544 end do 545 end do 546 end do 547 ia=ia+nattyp(itypat) 548 end do 549 end if 550 dotr=sqrt(dotr**2+doti**2) 551 dotr1=sqrt(dotr1**2+doti1**2) 552 write(std_out,*)'DOTR, DOTR1',dotr,dotr1 553 dotr=dotr1/dotr 554 write(std_out,*)'DOTR',dotr 555 deltawf1=zero 556 if(dotr>=0.9d0) then 557 deltawf1(:,:)=cg(:,icg+1+ibd:ibd+icg+inc)-wf1(:,:) 558 if(usepaw==1) then 559 alpha(1)=one;alpha(2)=zero 560 beta(1)=-one;beta(2)=zero 561 ia =0 562 call pawcprj_zaxpby(alpha,beta,cprj_k(:,iblockbd:iblockbd),cprj_k1(:,iblockbd:iblockbd)) 563 end if 564 istep1=istep 565 else 566 istep1=1 567 end if 568 scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)=cg(:,icg+1+ibd:ibd+icg+inc) 569 scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)=deltawf1(:,:) 570 if(usepaw==1) then 571 call pawcprj_put(atindx1,cprj_k,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,& 572 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 573 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 574 call pawcprj_put(atindx1,cprj_k1,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,& 575 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 576 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 577 end if 578 579 ! if(istep1>=3) then 580 cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:) & 581 & +scf_history%beta *deltawf2(:,:) 582 583 ! to be used later 584 ! if(usepaw==1) then 585 ! alpha(2)=zero 586 ! beta(1)=one;beta(2)=zero 587 ! alpha(1)=scf_history%alpha 588 ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 589 ! alpha(1)=scf_history%beta 590 ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 591 ! call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,& 592 ! & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 593 ! & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 594 ! end if 595 ! else if (istep1==2) then 596 ! cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:)+scf_history%beta*wf1(:,:) 597 ! ! cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+deltawf1(:,:) 598 ! if(usepaw==1) then 599 ! alpha(2)=zero 600 ! beta(1)=one;beta(2)=zero 601 ! alpha(1)=scf_history%alpha 602 ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 603 ! alpha(1)=scf_history%beta 604 ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 605 ! call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,& 606 ! & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 607 ! & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 608 ! end if 609 ! end if 610 ibd=ibd+inc 611 end do ! end loop on iblockbd 612 613 ABI_DEALLOCATE(deltawf1) 614 ABI_DEALLOCATE(deltawf2) 615 ABI_DEALLOCATE(wf1) 616 ABI_DEALLOCATE(cwavef) 617 ABI_DEALLOCATE(cwavef1) 618 ABI_DEALLOCATE(eig) 619 ABI_DEALLOCATE(evec) 620 ABI_DEALLOCATE(unm) 621 if(usepaw==1) then 622 call pawcprj_free(cprj_k) 623 ABI_DATATYPE_DEALLOCATE(cprj_k) 624 call pawcprj_free(cprj_k1) 625 ABI_DATATYPE_DEALLOCATE(cprj_k1) 626 call pawcprj_free(cprj_k2) 627 ABI_DATATYPE_DEALLOCATE(cprj_k2) 628 end if 629 630 ibg=ibg+my_nspinor*nband_k 631 icg=icg+my_nspinor*nband_k*npw_k 632 633 ! End big k point loop 634 end do 635 ! End loop over spins 636 end do 637 638 if(usepaw==1) then 639 call pawcprj_free(cprj) 640 ABI_DATATYPE_DEALLOCATE(cprj) 641 end if 642 if (nprocband>1) then 643 ABI_DEALLOCATE(npw_block) 644 ABI_DEALLOCATE(npw_disp) 645 ABI_DEALLOCATE(bufsize) 646 ABI_DEALLOCATE(bufdisp) 647 ABI_DEALLOCATE(bufsize_wf) 648 ABI_DEALLOCATE(bufdisp_wf) 649 end if 650 651 end if ! istep>=2 652 653 if (usepaw==1) then 654 ABI_DEALLOCATE(dimcprj) 655 end if 656 657 end subroutine extrapwf