TABLE OF CONTENTS
ABINIT/extraprho [ Functions ]
NAME
extraprho
FUNCTION
Extrapolate electronic density for new ionic positions from values of density of previous SCF cycle. Use algorithm proposed by D. Alfe in Comp. Phys. Comm. 118 (1999), 31-33 [[cite:Alfe1999]]
INPUTS
atindx atindx1(natom)=index table for atoms, inverse of atindx cg(2,mcg)= plane wave wavefunction coefficient cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn> dtset <type(dataset_type)>=all input variables in this dataset | densty(ntypat,4)=parameters for initialisation of the gaussian density | jellslab,slabzbeg,slabzend,slabwsrad=parameters for jellium slab | natom=number of atoms in cell. | nspden=number of spin-density components gmet(3,3)=reciprocal space metric gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff value on G**2 for sphere inside fft box istep=number of call the routine kg(3,mpw*mkmem)=reduced planewave coordinates. mgfft=maximum size of 1D FFTs mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization mqgrid=number of grid pts in q array for f(q) spline. my_natom=number of atoms treated by current processor nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) 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 ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information psps<type(pseudopotential_type)>=variables related to pseudopotentials qgrid(mqgrid)=q grid for spline from 0 to qmax rprimd(3,3)=dimensional primitive translation vectors (bohr) ucvol=unit cell volume (bohr**3). usepaw= 0 for non paw calculation; =1 for paw calculation xred_new(3,natom)=new reduced coordinates for atoms in unit cell 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 zion(ntypat)=charge on each type of atom znucl(ntypat)=atomic number of each atom type
SIDE EFFECTS
pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data Value from previous SCF cycle is input Extrapolated value is output rhor(nfft,nspden)=the density from previous SCF cycle is input the extrapolated density is output scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
SOURCE
116 subroutine extraprho(atindx,atindx1,cg,cprj,dtset,gmet,gprimd,gsqcut,istep,& 117 & kg,mcg,mcprj,mgfft,mpi_enreg,mqgrid,my_natom,nattyp,nfft,ngfft,npwarr,ntypat,pawrhoij,& 118 & pawtab,ph1d,psps,qgrid,rhor,rprimd,scf_history,ucvol,usepaw,& 119 & xred_new,xred_old,ylm,zion,znucl) 120 121 !Arguments ------------------------------------ 122 !scalars 123 integer,intent(in) :: istep,mcg,mcprj,mgfft,my_natom,mqgrid,nfft,ntypat,usepaw 124 real(dp),intent(in) :: gsqcut,ucvol 125 type(MPI_type),intent(in) :: mpi_enreg 126 type(dataset_type),intent(in) :: dtset 127 type(scf_history_type),intent(inout) :: scf_history 128 type(pseudopotential_type),intent(in) :: psps 129 !arrays 130 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem) 131 integer,intent(in) :: nattyp(ntypat),ngfft(18),npwarr(dtset%nkpt) 132 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom) 133 real(dp),intent(in) :: qgrid(mqgrid) 134 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 135 real(dp),intent(in) :: zion(ntypat),znucl(ntypat) 136 real(dp), intent(inout) :: cg(2,mcg) 137 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),rprimd(3,3),xred_new(3,dtset%natom) 138 real(dp),intent(in) :: xred_old(3,dtset%natom) 139 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw) 140 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 141 type(pawcprj_type),intent(inout) :: cprj(:,:) 142 143 !Local variables------------------------------- 144 !scalars 145 integer :: cplex_rhoij,dplex,iatom,ii,ind1,ind1new,ind2,ind2new,iq,iq0,irhoij,ispden,itypat,jrhoij,klmn 146 integer :: lmn2_size,nselect,nspden_rhoij,optatm,optdyfr,opteltfr,optgr,option,optn,optn2 147 integer :: optstr,optv,qphase_rhoij 148 real(dp) :: a11,a12,a22,a33,alpha,b1,b2,beta,detA,fact,ratio1,ratio2 149 logical :: hasmoved,usegauss 150 character(len=500) :: message 151 !arrays 152 integer :: dummy3(3) 153 real(dp) :: diff_t(3),diff_tmdt(3),diff_tpdt(3),dummy2(2) 154 real(dp) :: dummy_in(0) 155 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0) 156 real(dp) :: strn_dummy6(6),strv_dummy6(6) 157 real(dp),allocatable :: deltarho(:),gauss(:,:),rhoijtmp(:,:),work1(:) 158 real(dp),allocatable :: work2(:,:),work3(:,:),xred_tpdt(:,:) 159 160 ! ************************************************************************* 161 162 !--------------------------------------------------------------- 163 !----------- Inits 164 !--------------------------------------------------------------- 165 166 !History indexes 167 ind1=scf_history%hindex(1) 168 ind2=scf_history%hindex(2) 169 170 !Compatibility tests 171 if (ind1==0.and.ind2>0)then 172 ABI_BUG(' Incompatible history indexes !') 173 end if 174 175 !Rotated values of history indexes 176 if (ind1>0.and.ind2>0) then 177 ind1new=ind2;ind2new=ind1 178 else if (ind1>0.and.ind2==0) then 179 ind1new=3-ind1;ind2new=ind1 180 else if (ind1==0.and.ind2==0) then 181 ind1new=1;ind2new=0 182 end if 183 184 !Compute ionic positions at t+dt in red. coordinates 185 !Has to take the boundary conditions into account 186 ABI_MALLOC(xred_tpdt,(3,dtset%natom)) 187 do iatom=1,dtset%natom 188 xred_tpdt(1,iatom)=xred_old(1,iatom)+mod(xred_new(1,iatom)-xred_old(1,iatom)+1.5_dp,one)-half 189 xred_tpdt(2,iatom)=xred_old(2,iatom)+mod(xred_new(2,iatom)-xred_old(2,iatom)+1.5_dp,one)-half 190 xred_tpdt(3,iatom)=xred_old(3,iatom)+mod(xred_new(3,iatom)-xred_old(3,iatom)+1.5_dp,one)-half 191 end do 192 193 !--------------------------------------------------------------- 194 !----------- Compute Alpha and Beta 195 !----------- see (4) in Comp. Phys. Comm. 118 (1999), 31-33 [[cite:Alfe1999]] 196 !--------------------------------------------------------------- 197 198 !Compute a_ij matrix 199 a11=zero;a12=zero;a22=zero;a33=zero;b1=zero;b2=zero 200 diff_t=zero;diff_tmdt=zero;diff_tpdt=zero 201 do iatom=1,dtset%natom 202 203 diff_tpdt(1:3)=xred_tpdt(1:3,iatom)-xred_old(1:3,iatom) 204 if (ind1>0) then 205 diff_t(1:3)=scf_history%xreddiff(1:3,iatom,ind1) 206 if (ind2>0) diff_tmdt(1:3)=scf_history%xreddiff(1:3,iatom,ind2) 207 end if 208 do ii=1,3 209 a11=a11+diff_t(ii)**2 210 a22=a22+diff_tmdt(ii)**2 211 a33=a33+diff_tpdt(ii)**2 212 a12=a12+diff_t(ii) *diff_tmdt(ii) 213 b1 =b1 +diff_t(ii) *diff_tpdt(ii) 214 b2 =b2 +diff_tmdt(ii)*diff_tpdt(ii) 215 end do 216 217 ! Store reduced coordinates diffs in SCF history 218 scf_history%xreddiff(1:3,iatom,ind1new)=diff_tpdt(1:3) 219 220 end do 221 ABI_FREE(xred_tpdt) 222 hasmoved=(a11>=tol10.or.a22>=tol10.or.a33>=tol10) 223 224 !Compute alpha and beta 225 alpha=zero;beta=zero 226 if (hasmoved.and.ind1>0) then 227 ratio1=one;if (abs(a33)>=tol10) ratio1=(a11+a33-two*b1)/a33 228 ratio2=one;if (abs(a33)>=tol10) ratio2=(a11+a33-two*b2)/a33 229 detA=a11*a22-a12**2 230 if (abs(a11)>=tol10.and.(abs(a22)<tol10.or.abs(detA)<tol10)) then 231 alpha=b1/a11 232 else if (abs(a22)>=tol10.and.(abs(a11)<tol10.or.abs(detA)<tol10)) then 233 beta=b2/a22 234 else if (abs(ratio1)+abs(ratio2)<tol6) then 235 if (ind2>0) then 236 alpha=two;beta=-one 237 else 238 alpha=one 239 end if 240 write(message,'(3a,f4.1,a,f4.1)')& 241 & 'Ionic positions lead to a collinear system !',ch10,& 242 & 'Mixing coeffs have been set to: alpha=',alpha,' beta=',beta 243 ABI_WARNING(message) 244 else if (abs(a11)>=tol10.and.abs(a22)>=tol10) then 245 alpha=(b1*a22-b2*a12)/detA 246 beta =(b2*a11-b1*a12)/detA 247 end if 248 end if 249 250 251 !--------------------------------------------------------------- 252 !----------- Contribution from delta_rho(t), delta_rho(t-dt) 253 !----------- and delta_rho(t-2dt) to predicted rho(t+dt) 254 !--------------------------------------------------------------- 255 256 !deltarho(t+dt) <- deltarho(t) + alpha.[deltarho(t)-deltarho(t-dt)] 257 !+ beta .[deltarho(t-dt)-deltarho(t-2dt)] 258 !Note: scf_history%deltarhor is updated at the same time 259 260 ABI_MALLOC(deltarho,(nfft)) 261 do ispden=1,dtset%nspden 262 263 if (ispden==1) then 264 deltarho(:)=rhor(:,ispden)-scf_history%atmrho_last(:) 265 else if (ispden==2.and.dtset%nspden==2) then 266 deltarho(:)=rhor(:,ispden)-half*scf_history%atmrho_last(:) 267 end if 268 269 270 ! rho(t+dt) <- deltarho(t) + alpha.deltarho(t) 271 if (dtset%nspden/=4.or.ispden==1) then 272 rhor(:,ispden)=(one+alpha)*deltarho(:) 273 else 274 rhor(:,ispden)=(one+alpha)*rhor(:,ispden) 275 end if 276 277 if (hasmoved) then 278 279 ! rho(t+dt) <- -alpha.deltarho(t-dt) + beta.deltarho(t-dt) 280 if (abs(beta-alpha)>tol14.and.ind1>0) then 281 rhor(:,ispden)=rhor(:,ispden)+(beta-alpha)*scf_history%deltarhor(:,ispden,ind1) 282 end if 283 284 ! rho(t+dt) <- -beta.deltarho(t-2dt) 285 if (abs(beta)>tol14.and.ind2>0) then 286 rhor(:,ispden)=rhor(:,ispden)-beta*scf_history%deltarhor(:,ispden,ind2) 287 end if 288 289 end if 290 291 ! Store deltarho(t) in history 292 if (dtset%nspden/=4.or.ispden==1) then 293 scf_history%deltarhor(:,ispden,ind1new)=deltarho(:) 294 else 295 scf_history%deltarhor(:,ispden,ind1new)=rhor(:,ispden) 296 end if 297 298 end do 299 300 ABI_FREE(deltarho) 301 302 !--------------------------------------------------------------- 303 !----------- Contribution from rho_at(t+dt) to predicted rho(t+dt) 304 !--------------------------------------------------------------- 305 306 !Determine whether a gaussian atomic density has to be used or not 307 !MG: Note that there's a small inconsistency between initro and extraprho because in initrho 308 ! we use `use_gaussian(ntypat)`. 309 usegauss=.true. 310 if (usepaw==0) usegauss = any(.not.psps%nctab(1:ntypat)%has_tvale) 311 if (usepaw==1) usegauss=(minval(pawtab(1:ntypat)%has_tvale)==0) 312 if (usegauss) then 313 optn2=3 314 ABI_MALLOC(gauss,(2,ntypat)) 315 do itypat=1,ntypat 316 gauss(1,itypat)=zion(itypat) 317 gauss(2,itypat) = atom_length(dtset%densty(itypat,1),zion(itypat),znucl(itypat)) 318 end do 319 call wrtout(std_out," Extrapolating rho(t+dt) using gaussian functions as atomic densities", "COLL") 320 else 321 optn2=2 322 ABI_MALLOC(gauss,(2,0)) 323 call wrtout(std_out," Extrapolating rho(t+dt) using atomic densities taken from pseudos", "COLL") 324 end if 325 326 !Compute rho_at(t+dt) as sum of atomic densities 327 !Note: scf_history%atmrho_last is updated at the same time 328 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=0;optn=1 329 call atm2fft(atindx1,scf_history%atmrho_last,dummy_out1,dummy_out2,dummy_out3,& 330 & dummy_out4,gauss,gmet,gprimd,dummy_out5,dummy_out6,gsqcut,mgfft,mqgrid,dtset%natom,nattyp,& 331 & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,& 332 & dummy3,dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,usepaw,dummy_in,dummy_in,dummy_in,dummy2,dummy_in,& 333 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 334 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 335 ABI_FREE(gauss) 336 337 !Take eventually into account jellium slab 338 if (dtset%jellslab/=0) then 339 option=2 340 ABI_MALLOC(work1,(nfft)) 341 ABI_MALLOC(work2,(nfft,1)) 342 ABI_MALLOC(work3,(2,nfft)) 343 work2(:,1)=scf_history%atmrho_last(:) 344 call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,1,option,& 345 & dtset%slabwsrad,work3,work2,rprimd,work1,dtset%slabzbeg,dtset%slabzend) 346 scf_history%atmrho_last(:)=work2(:,1) 347 ABI_FREE(work1) 348 ABI_FREE(work2) 349 ABI_FREE(work3) 350 end if 351 352 !Add rho_at(t+dt) to rho(t+dt) 353 rhor(:,1)=rhor(:,1)+scf_history%atmrho_last(:) 354 if (dtset%nspden==2) rhor(:,2)=rhor(:,2)+half*scf_history%atmrho_last(:) 355 356 !--------------------------------------------------------------- 357 !----------- Extrapolation of PAW rhoij occupancy matrixes 358 !--------------------------------------------------------------- 359 360 if (usepaw==1) then 361 362 if (ind2==0) then 363 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 364 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc) 365 call pawrhoij_alloc(scf_history%pawrhoij(:,ind1new),cplex_rhoij,nspden_rhoij,& 366 & dtset%nspinor,dtset%nsppol,dtset%typat,pawtab=pawtab,& 367 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 368 end if 369 370 do iatom=1,my_natom 371 372 nspden_rhoij=pawrhoij(iatom)%nspden 373 lmn2_size=pawrhoij(iatom)%lmn2_size 374 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij;dplex=cplex_rhoij-1 375 qphase_rhoij=pawrhoij(iatom)%qphase 376 377 if (hasmoved) then 378 ABI_MALLOC(rhoijtmp,(cplex_rhoij*qphase_rhoij*lmn2_size,nspden_rhoij)) 379 rhoijtmp=zero 380 381 do ispden=1,nspden_rhoij 382 do iq=1,qphase_rhoij 383 iq0=merge(0,cplex_rhoij*lmn2_size,iq==1) 384 385 ! rhoij(t+dt) <- rhoij(t) + alpha.rhoij(t) 386 fact=one+alpha 387 jrhoij=1+iq0 388 do irhoij=1,pawrhoij(iatom)%nrhoijsel 389 klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex+iq0 390 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 391 & +fact*pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 392 jrhoij=jrhoij+cplex_rhoij 393 end do 394 395 ! rhoij(t+dt) <- -alpha.rhoij(t-dt) + beta.rhoij(t-dt) 396 if (abs(beta-alpha)>tol14.and.ind1>0) then 397 fact=beta-alpha 398 jrhoij=1+iq0 399 do irhoij=1,scf_history%pawrhoij(iatom,ind1)%nrhoijsel 400 klmn=cplex_rhoij*scf_history%pawrhoij(iatom,ind1)%rhoijselect(irhoij)-dplex+iq0 401 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 402 & +fact*scf_history%pawrhoij(iatom,ind1)%rhoijp(jrhoij:jrhoij+dplex,ispden) 403 jrhoij=jrhoij+cplex_rhoij 404 end do 405 end if 406 407 ! rho(t+dt) <- -beta.rhoij(t-2dt) 408 if (abs(beta)>tol14.and.ind2>0) then 409 fact=-beta 410 jrhoij=1+iq0 411 do irhoij=1,scf_history%pawrhoij(iatom,ind2)%nrhoijsel 412 klmn=cplex_rhoij*scf_history%pawrhoij(iatom,ind2)%rhoijselect(irhoij)-dplex+iq0 413 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 414 & +fact*scf_history%pawrhoij(iatom,ind2)%rhoijp(jrhoij:jrhoij+dplex,ispden) 415 jrhoij=jrhoij+cplex_rhoij 416 end do 417 end if 418 419 end do ! iq 420 end do !ispden 421 end if !hasmoved 422 423 ! Store rhoij(t) in history 424 ! (cannot use pawrhoij_copy here because update for single atom) 425 nselect=pawrhoij(iatom)%nrhoijsel 426 scf_history%pawrhoij(iatom,ind1new)%nrhoijsel=nselect 427 scf_history%pawrhoij(iatom,ind1new)%rhoijselect(:)=0 428 scf_history%pawrhoij(iatom,ind1new)%rhoijselect(1:nselect)=pawrhoij(iatom)%rhoijselect(1:nselect) 429 scf_history%pawrhoij(iatom,ind1new)%rhoijp(1:cplex_rhoij*nselect,1:nspden_rhoij)= & 430 & pawrhoij(iatom)%rhoijp(1:cplex_rhoij*nselect,1:nspden_rhoij) 431 432 ! Select non-zero values of rhoij(t+dt) 433 if (hasmoved) then 434 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,pawrhoij(iatom)%nrhoijsel,& 435 & cplex_rhoij,qphase_rhoij,lmn2_size,nspden_rhoij,& 436 & rhoij_input=rhoijtmp) 437 ABI_FREE(rhoijtmp) 438 end if 439 440 end do !iatom 441 end if !usepaw 442 443 444 scf_history%alpha=alpha 445 scf_history%beta=beta 446 447 448 449 !--------------------------------------------------------------- 450 !----------- End 451 !--------------------------------------------------------------- 452 453 if(scf_history%usecg==1) then 454 if (hasmoved) then 455 if (dtset%extrapwf==1) then 456 call extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,nattyp,& 457 & ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm) 458 elseif(dtset%extrapwf==2) then 459 scf_history%hindex(1)=ind1 460 scf_history%hindex(2)=ind2 461 scf_history%hindex(3)=ind1new 462 call extrapwf_biortho(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,& 463 & nattyp,npwarr,pawtab,scf_history) 464 end if 465 else 466 scf_history%cg(:,:,2)=zero 467 end if 468 469 end if 470 !Rotate history indexes 471 scf_history%hindex(1)=ind1new 472 scf_history%hindex(2)=ind2new 473 474 475 end subroutine extraprho
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) [[cite:Arias1992]]
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).
SOURCE
520 subroutine extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,& 521 & nattyp,ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm) 522 523 !Arguments ------------------------------------ 524 !scalars 525 integer,intent(in) :: istep,mcg,mgfft,ntypat,usepaw 526 type(MPI_type),intent(in) :: mpi_enreg 527 type(dataset_type),intent(in) :: dtset 528 type(scf_history_type),intent(inout) :: scf_history 529 type(pseudopotential_type),intent(in) :: psps 530 !arrays 531 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfft(18) 532 integer,intent(in) :: npwarr(dtset%nkpt) 533 real(dp),intent(in) :: rprimd(3,3) 534 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 535 real(dp), intent(inout) :: cg(2,mcg) 536 real(dp),intent(in) :: xred_old(3,dtset%natom) 537 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 538 539 540 !Local variables------------------------------- 541 !scalars 542 integer :: ia,iat,iatom,iband_max,iband_max1,iband_min,iband_min1,ibd,ibg,iblockbd,iblockbd1,icg,icgb,icgb1,icgb2 543 integer :: ierr,ig,ii,ikpt,ilmn1,ilmn2,inc,ind1,ind2,iorder_cprj 544 integer :: isize,isppol,istep1,istwf_k,itypat,klmn,me_distrb,my_nspinor 545 integer :: nband_k,nblockbd,nprocband,npw_k,npw_nk,spaceComm_band 546 real(dp) :: dotr,dotr1,doti,doti1,eigval 547 !character(len=500) :: message 548 !arrays 549 real(dp) :: alpha(2),beta(2),gmet(3,3),gprimd(3,3),rmet(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom),ucvol 550 integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:),dimcprj(:),npw_block(:),npw_disp(:) 551 real(dp),allocatable :: al(:,:),anm(:),cwavef(:,:),cwavef1(:,:),cwavef_tmp(:,:),deltawf1(:,:),deltawf2(:,:) 552 real(dp),allocatable :: eig(:),evec(:,:) 553 real(dp),allocatable :: unm(:,:,:) 554 real(dp),allocatable :: work(:,:),work1(:,:),wf1(:,:),ylmgr_k(:,:,:),zhpev1(:,:),zhpev2(:) 555 complex(dpc),allocatable :: unm_tmp(:,:),anm_tmp(:,:) 556 type(pawcprj_type),allocatable :: cprj(:,:),cprj_k(:,:),cprj_k1(:,:),cprj_k2(:,:),cprj_k3(:,:),cprj_k4(:,:) 557 !complex(dpc) :: aa 558 559 ! ************************************************************************* 560 561 if (istep==0) return 562 563 !Useful array 564 if (usepaw==1) then 565 ABI_MALLOC(dimcprj,(dtset%natom)) 566 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O') 567 end if 568 569 !Metric 570 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 571 572 !History indexes 573 ind1=1;ind2=2 574 575 !First step 576 if (istep==1) then 577 scf_history%cg(:,:,ind1)=cg(:,:) 578 ! scf_history%cg(:,:,ind2)=zero 579 scf_history%cg(:,:,ind2)= cg(:,:) 580 if(usepaw==1) then 581 ! WARNING: THIS SECTION IS USELESS; NOW crpj CAN BE READ FROM SCFCV 582 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old) 583 iatom=0 ; iorder_cprj=0 584 call pawcprj_alloc(scf_history%cprj(:,:,ind1),0,dimcprj) 585 call pawcprj_alloc(scf_history%cprj(:,:,ind2),0,dimcprj) 586 ABI_MALLOC(ylmgr_k,(dtset%mpw,3,0)) 587 call ctocprj(atindx,cg,1,scf_history%cprj(:,:,ind1),gmet,gprimd,& 588 & iatom,0,iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,& 589 & dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,& 590 & dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,& 591 & dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,dtset%ntypat,& 592 & dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,0,& 593 & xred_old,ylm,ylmgr_k) 594 ABI_FREE(ylmgr_k) 595 ! call pawcprj_set_zero(scf_history%cprj(:,:,ind2)) 596 call pawcprj_copy(scf_history%cprj(:,:,ind1),scf_history%cprj(:,:,ind2)) 597 end if 598 else 599 600 !From 2nd step 601 602 ! Init parallelism 603 me_distrb=mpi_enreg%me_kpt 604 if (mpi_enreg%paral_kgb==1.or.mpi_enreg%paralbd==1) then 605 spaceComm_band=mpi_enreg%comm_band 606 nprocband=mpi_enreg%nproc_band 607 else 608 spaceComm_band=xmpi_comm_self 609 nprocband=1 610 end if 611 612 ! For the moment sequential part only 613 nprocband=1 614 615 ! Additional statements if band-fft parallelism 616 if (nprocband>1) then 617 ABI_MALLOC(npw_block,(nprocband)) 618 ABI_MALLOC(npw_disp,(nprocband)) 619 ABI_MALLOC(bufsize,(nprocband)) 620 ABI_MALLOC(bufdisp,(nprocband)) 621 ABI_MALLOC(bufsize_wf,(nprocband)) 622 ABI_MALLOC(bufdisp_wf,(nprocband)) 623 end if 624 625 icg=0 626 ibg=0 627 628 if(usepaw==1) then 629 ! WARNING: THIS SECTION IS USELESS; NOW cprj CAN BE READ FROM SCFCV 630 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred_old) 631 ABI_MALLOC(cprj,(dtset%natom,scf_history%mcprj)) 632 call pawcprj_alloc(cprj,0,dimcprj) 633 iatom=0 ; iorder_cprj=0 634 ABI_MALLOC(ylmgr_k,(dtset%mpw,3,0)) 635 call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,iorder_cprj,& 636 & dtset%istwfk,kg,dtset%kptns,mcg,scf_history%mcprj,dtset%mgfft,& 637 & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,& 638 & nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,& 639 & npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,& 640 & ph1d,psps,rmet,dtset%typat,ucvol,0,xred_old,& 641 & ylm,ylmgr_k) 642 ABI_FREE(ylmgr_k) 643 end if ! end usepaw=1 644 645 ! LOOP OVER SPINS 646 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 647 do isppol=1,dtset%nsppol 648 649 ! BIG FAT k POINT LOOP 650 do ikpt=1,dtset%nkpt 651 652 ! Select k point to be treated by this proc 653 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 654 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 655 656 istwf_k=dtset%istwfk(ikpt) 657 658 ! Retrieve number of plane waves 659 npw_k=npwarr(ikpt) 660 if (nprocband>1) then 661 ! Special treatment for band-fft // 662 call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr) 663 npw_nk=sum(npw_block);npw_disp(1)=0 664 do ii=2,nprocband 665 npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1) 666 end do 667 else 668 npw_nk=npw_k 669 end if 670 671 ! Allocate arrays for a wave-function (or a block of WFs) 672 ABI_MALLOC(cwavef,(2,npw_nk*my_nspinor)) 673 ABI_MALLOC(cwavef1,(2,npw_nk*my_nspinor)) 674 if (nprocband>1) then 675 isize=2*my_nspinor;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:) 676 isize=2*my_nspinor*npw_k;bufsize_wf(:)=isize 677 do ii=1,nprocband 678 bufdisp_wf(ii)=(ii-1)*isize 679 end do 680 end if 681 682 ! Subspace alignment 683 684 ! Loop over bands or blocks of bands 685 nblockbd=nband_k/nprocband 686 icgb=icg 687 688 if(usepaw==1) then 689 ABI_MALLOC( cprj_k,(dtset%natom,my_nspinor*nblockbd)) 690 call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimcprj) 691 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,& 692 & dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,& 693 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 694 ABI_MALLOC( cprj_k1,(dtset%natom,my_nspinor*nblockbd)) 695 call pawcprj_alloc(cprj_k1,scf_history%cprj(1,1,ind1)%ncpgr,dimcprj) 696 call pawcprj_get(atindx1,cprj_k1,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,& 697 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,& 698 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 699 ABI_MALLOC( cprj_k2,(dtset%natom,my_nspinor*nblockbd)) 700 call pawcprj_alloc(cprj_k2,scf_history%cprj(1,1,ind2)%ncpgr,dimcprj) 701 call pawcprj_get(atindx1,cprj_k2,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,& 702 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,my_nspinor,dtset%nsppol,0,& 703 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 704 end if !end usepaw=1 705 706 ABI_MALLOC(unm,(2,nblockbd,nblockbd)) 707 unm=zero 708 icgb2=0 709 710 do iblockbd=1,nblockbd 711 iband_min=1+(iblockbd-1)*nprocband 712 iband_max=iblockbd*nprocband 713 714 if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then 715 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) cycle 716 end if 717 718 ! Extract wavefunction information 719 if (nprocband>1) then 720 ! Special treatment for band-fft // 721 ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*nprocband)) 722 do ig=1,npw_k*my_nspinor*nprocband 723 cwavef_tmp(1,ig)=cg(1,ig+icgb) 724 cwavef_tmp(2,ig)=cg(2,ig+icgb) 725 end do 726 call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr) 727 ABI_FREE(cwavef_tmp) 728 else 729 do ig=1,npw_k*my_nspinor 730 cwavef(1,ig)=cg(1,ig+icgb) 731 cwavef(2,ig)=cg(2,ig+icgb) 732 end do 733 end if 734 735 icgb1=icg 736 737 do iblockbd1=1,nblockbd 738 iband_min1=1+(iblockbd1-1)*nprocband 739 iband_max1=iblockbd1*nprocband 740 741 if(xmpi_paral==1.and.mpi_enreg%paral_kgb/=1) then 742 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min1,iband_max1,isppol,me_distrb)) cycle 743 end if 744 745 ! Extract wavefunction information 746 747 if (nprocband>1) then 748 ! Special treatment for band-fft // 749 ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*nprocband)) 750 do ig=1,npw_k*my_nspinor*nprocband 751 cwavef_tmp(1,ig)=scf_history%cg(1,ig+icgb1,ind1) 752 cwavef_tmp(2,ig)=scf_history%cg(2,ig+icgb1,ind1) 753 end do 754 call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef1,bufsize,bufdisp,spaceComm_band,ierr) 755 ABI_FREE(cwavef_tmp) 756 else 757 do ig=1,npw_k*my_nspinor 758 cwavef1(1,ig)=scf_history%cg(1,ig+icgb1,ind1) 759 cwavef1(2,ig)=scf_history%cg(2,ig+icgb1,ind1) 760 end do 761 end if 762 763 ! Calculate Unm=<psi_nk(t)|S|psi_mk(t-dt)> 764 call dotprod_g(dotr,doti,istwf_k,npw_k*my_nspinor,2,cwavef,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 765 if(usepaw==1) then 766 ia =0 767 do itypat=1,ntypat 768 do iat=1+ia,nattyp(itypat)+ia 769 do ilmn1=1,pawtab(itypat)%lmn_size 770 do ilmn2=1,ilmn1 771 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 772 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+& 773 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)) 774 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-& 775 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)) 776 end do 777 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 778 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 779 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)+& 780 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)) 781 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd1)%cp(2,ilmn2)-& 782 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd1)%cp(1,ilmn2)) 783 end do 784 end do 785 end do 786 ia=ia+nattyp(itypat) 787 end do 788 end if 789 ! unm(1,iblockbd,iblockbd1)=dotr 790 ! unm(2,iblockbd,iblockbd1)=doti 791 unm(1,iblockbd1,iblockbd)=dotr 792 unm(2,iblockbd1,iblockbd)=doti 793 ! End loop over bands iblockbd1 794 icgb1=icgb1+npw_k*my_nspinor*nprocband 795 796 end do 797 798 ! End loop over bands iblockbd 799 icgb2=icgb2+npw_k*my_nspinor*nprocband 800 icgb=icgb+npw_k*my_nspinor*nprocband 801 end do 802 803 ! write(std_out,*) 'UNM' 804 ! do iblockbd=1,nblockbd 805 ! write(std_out,11) (unm(1,iblockbd,iblockbd1),unm(2,iblockbd,iblockbd1),iblockbd1=1,nblockbd) 806 ! end do 807 ! 11 format(12(1x,f9.5),a) 808 ! Compute A=tU^*U 809 ABI_MALLOC(unm_tmp,(nblockbd,nblockbd)) 810 ABI_MALLOC(anm_tmp,(nblockbd,nblockbd)) 811 ABI_MALLOC(anm,(nblockbd*(nblockbd+1))) 812 unm_tmp(:,:)=cmplx(unm(1,:,:),unm(2,:,:),kind=dp) 813 call zgemm('C','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp), unm_tmp,nblockbd, & 814 & unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd) 815 do iblockbd=1,nblockbd 816 do iblockbd1=iblockbd,nblockbd 817 ii=iblockbd1*(iblockbd1-1)+2*(iblockbd-1)+1 818 anm(ii)=real(anm_tmp(iblockbd,iblockbd1)) 819 anm(ii+1)=aimag(anm_tmp(iblockbd,iblockbd1)) 820 end do 821 end do 822 call hermit(anm,anm,ierr,nblockbd) 823 ! aa=dcmplx(0._dp) 824 ! do iblockbd=1,nblockbd 825 ! aa=aa+conjg(unm_tmp(iblockbd,1))*unm_tmp(iblockbd,1) 826 ! end do 827 ! write(std_out,*) 'tU*U', aa 828 ! write(std_out,*) 'ANM_tmp' 829 ! do iblockbd=1,nblockbd 830 ! write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd) 831 ! end do 832 ! write(std_out,*) 'ANM' 833 ! do iblockbd=1,nblockbd*(nblockbd+1) 834 ! write(std_out,11) anm(iblockbd) 835 ! end do 836 837 ! Diagonalize A 838 ABI_MALLOC(eig,(nblockbd)) 839 ABI_MALLOC(evec,(2*nblockbd,nblockbd)) 840 ABI_MALLOC(zhpev1,(2,2*nblockbd-1)) 841 ABI_MALLOC(zhpev2,(3*nblockbd-2)) 842 call zhpev('V','U',nblockbd,anm,eig,evec,nblockbd,zhpev1,& 843 & zhpev2,ierr) 844 ABI_FREE(anm) 845 ABI_FREE(zhpev1) 846 ABI_FREE(zhpev2) 847 ! aa=dcmplx(0._dp) 848 ! do iblockbd=1,nblockbd 849 ! aa=aa+anm_tmp(1,iblockbd)*cmplx(evec((2*iblockbd-1),1),evec(2*iblockbd,1),kind=dp) 850 ! end do 851 ! write(std_out,*) 'EIG', aa, eig(1)*evec(1,1),eig(1)*evec(2,1) 852 853 ! Compute A'=evec*tU^/sqrt(eig) 854 call zgemm('C','C',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, & 855 & unm_tmp,nblockbd,dcmplx(0._dp),anm_tmp,nblockbd) 856 do iblockbd=1,nblockbd 857 eigval=dsqrt(eig(iblockbd)) 858 do iblockbd1=1,nblockbd 859 anm_tmp(iblockbd,iblockbd1)=anm_tmp(iblockbd,iblockbd1)/eigval 860 end do 861 end do 862 863 ! Compute tA^A'to come back to the initial subspace for the cg's 864 865 call zgemm('N','N',nblockbd,nblockbd,nblockbd,dcmplx(1._dp),evec,nblockbd, & 866 & anm_tmp,nblockbd,dcmplx(0._dp),unm_tmp,nblockbd) 867 anm_tmp=unm_tmp 868 ! write(std_out,*) 'ANM_tmp' 869 ! do iblockbd=1,nblockbd 870 ! write(std_out,11) (anm_tmp(iblockbd,iblockbd1),iblockbd1=1,nblockbd) 871 ! end do 872 873 ! Wavefunction alignment (istwfk=1 ?) 874 ABI_MALLOC(work,(2,npw_nk*my_nspinor*nblockbd)) 875 ABI_MALLOC(work1,(2,my_nspinor*nblockbd*npw_nk)) 876 work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind1) 877 call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), & 878 & work1,npw_nk*my_nspinor, & 879 & anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor) 880 scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind1)=work(:,:) 881 882 work1(:,:)=scf_history%cg(:,icg+1:icg+my_nspinor*nblockbd*npw_nk,ind2) 883 call zgemm('N','N',npw_nk*my_nspinor,nblockbd,nblockbd,dcmplx(1._dp), & 884 & work1,npw_nk*my_nspinor, & 885 & anm_tmp,nblockbd,dcmplx(0._dp),work,npw_nk*my_nspinor) 886 scf_history%cg(:,1+icg:npw_nk*my_nspinor*nblockbd+icg,ind2)=work(:,:) 887 ABI_FREE(work1) 888 ! If paw, must also align cprj: 889 if (usepaw==1) then 890 ! New version (MT): 891 ABI_MALLOC(cprj_k3,(dtset%natom,my_nspinor)) 892 ABI_MALLOC(cprj_k4,(dtset%natom,my_nspinor)) 893 call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj) 894 call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj) 895 ABI_MALLOC(al,(2,nblockbd)) 896 do iblockbd=1,nblockbd 897 ii=(iblockbd-1)*my_nspinor 898 do iblockbd1=1,nblockbd 899 al(1,iblockbd1)=real (anm_tmp(iblockbd,iblockbd1)) 900 al(2,iblockbd1)=aimag(anm_tmp(iblockbd,iblockbd1)) 901 end do 902 call pawcprj_lincom(al,cprj_k1,cprj_k3,nblockbd) 903 call pawcprj_lincom(al,cprj_k2,cprj_k4,nblockbd) 904 call pawcprj_copy(cprj_k3,cprj_k1(:,ii+1:ii+my_nspinor)) 905 call pawcprj_copy(cprj_k4,cprj_k2(:,ii+1:ii+my_nspinor)) 906 end do 907 ABI_FREE(al) 908 ! Old version (FJ): 909 ! allocate( cprj_k3(dtset%natom,my_nspinor*nblockbd)) 910 ! call pawcprj_alloc(cprj_k3,cprj_k1(1,1)%ncpgr,dimcprj) 911 ! allocate( cprj_k4(dtset%natom,my_nspinor*nblockbd)) 912 ! call pawcprj_alloc(cprj_k4,cprj_k2(1,1)%ncpgr,dimcprj) 913 ! beta(1)=one;beta(2)=zero 914 ! do iblockbd=1,nblockbd*my_nspinor 915 ! do iblockbd1=1,nblockbd*my_nspinor 916 ! alpha(1)=real(anm_tmp(iblockbd,iblockbd1));alpha(2)=aimag(anm_tmp(iblockbd,iblockbd1)) 917 ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd1:iblockbd1),cprj_k3(:,iblockbd:iblockbd)) 918 ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd1:iblockbd1),cprj_k4(:,iblockbd:iblockbd)) 919 ! end do 920 ! end do 921 ! call pawcprj_copy(cprj_k3,cprj_k1) 922 ! call pawcprj_copy(cprj_k4,cprj_k2) 923 924 call pawcprj_free(cprj_k3) 925 call pawcprj_free(cprj_k4) 926 ABI_FREE(cprj_k3) 927 ABI_FREE(cprj_k4) 928 end if 929 ABI_FREE(anm_tmp) 930 ABI_FREE(unm_tmp) 931 ABI_FREE(work) 932 933 ! Wavefunction extrapolation 934 ibd=0 935 inc=npw_nk*my_nspinor 936 ABI_MALLOC(deltawf2,(2,npw_nk*my_nspinor)) 937 ABI_MALLOC(wf1,(2,npw_nk*my_nspinor)) 938 ABI_MALLOC(deltawf1,(2,npw_nk*my_nspinor)) 939 do iblockbd=1,nblockbd 940 deltawf2(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2) 941 wf1(:,:)=scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1) 942 ! wf1(2,1)=zero;deltawf2(2,1)=zero 943 944 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),& 945 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 946 call dotprod_g(dotr1,doti1,istwf_k,npw_nk*my_nspinor,2,cg(:,icg+1+ibd:ibd+icg+inc),wf1,& 947 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 948 if(usepaw==1) then 949 ia =0 950 do itypat=1,ntypat 951 do iat=1+ia,nattyp(itypat)+ia 952 do ilmn1=1,pawtab(itypat)%lmn_size 953 do ilmn2=1,ilmn1 954 klmn=((ilmn1-1)*ilmn1)/2+ilmn2 955 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+& 956 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)) 957 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-& 958 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)) 959 dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+& 960 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)) 961 doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-& 962 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)) 963 end do 964 do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size 965 klmn=((ilmn2-1)*ilmn2)/2+ilmn1 966 dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)+& 967 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)) 968 doti=doti+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k(iat,iblockbd)%cp(2,ilmn2)-& 969 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k(iat,iblockbd)%cp(1,ilmn2)) 970 dotr1=dotr1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)+& 971 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)) 972 doti1=doti1+pawtab(itypat)%sij(klmn)*(cprj_k(iat,iblockbd)%cp(1,ilmn1)*cprj_k1(iat,iblockbd)%cp(2,ilmn2)-& 973 & cprj_k(iat,iblockbd)%cp(2,ilmn1)*cprj_k1(iat,iblockbd)%cp(1,ilmn2)) 974 end do 975 end do 976 end do 977 ia=ia+nattyp(itypat) 978 end do 979 end if 980 dotr=sqrt(dotr**2+doti**2) 981 dotr1=sqrt(dotr1**2+doti1**2) 982 write(std_out,*)'DOTR, DOTR1',dotr,dotr1 983 dotr=dotr1/dotr 984 write(std_out,*)'DOTR',dotr 985 deltawf1=zero 986 if(dotr>=0.9d0) then 987 deltawf1(:,:)=cg(:,icg+1+ibd:ibd+icg+inc)-wf1(:,:) 988 if(usepaw==1) then 989 alpha(1)=one;alpha(2)=zero 990 beta(1)=-one;beta(2)=zero 991 ia =0 992 call pawcprj_zaxpby(alpha,beta,cprj_k(:,iblockbd:iblockbd),cprj_k1(:,iblockbd:iblockbd)) 993 end if 994 istep1=istep 995 else 996 istep1=1 997 end if 998 scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind1)=cg(:,icg+1+ibd:ibd+icg+inc) 999 scf_history%cg(:,1+icg+ibd:icg+ibd+inc,ind2)=deltawf1(:,:) 1000 if(usepaw==1) then 1001 call pawcprj_put(atindx1,cprj_k,scf_history%cprj(:,:,ind1),dtset%natom,1,ibg,ikpt,1,isppol,& 1002 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 1003 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1004 call pawcprj_put(atindx1,cprj_k1,scf_history%cprj(:,:,ind2),dtset%natom,1,ibg,ikpt,1,isppol,& 1005 & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 1006 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1007 end if 1008 1009 ! if(istep1>=3) then 1010 cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:) & 1011 & +scf_history%beta *deltawf2(:,:) 1012 1013 ! to be used later 1014 ! if(usepaw==1) then 1015 ! alpha(2)=zero 1016 ! beta(1)=one;beta(2)=zero 1017 ! alpha(1)=scf_history%alpha 1018 ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 1019 ! alpha(1)=scf_history%beta 1020 ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 1021 ! call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,& 1022 ! & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 1023 ! & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1024 ! end if 1025 ! else if (istep1==2) then 1026 ! cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+scf_history%alpha*deltawf1(:,:)+scf_history%beta*wf1(:,:) 1027 ! ! cg(:,icg+1+ibd:ibd+icg+inc)=cg(:,icg+1+ibd:ibd+icg+inc)+deltawf1(:,:) 1028 ! if(usepaw==1) then 1029 ! alpha(2)=zero 1030 ! beta(1)=one;beta(2)=zero 1031 ! alpha(1)=scf_history%alpha 1032 ! call pawcprj_zaxpby(alpha,beta,cprj_k1(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 1033 ! alpha(1)=scf_history%beta 1034 ! call pawcprj_zaxpby(alpha,beta,cprj_k2(:,iblockbd:iblockbd),cprj_k(:,iblockbd:iblockbd)) 1035 ! call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,& 1036 ! & dtset%mband,dtset%mkmem,dtset%natom,nblockbd,nblockbd,dimcprj,my_nspinor,dtset%nsppol,0,& 1037 ! & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1038 ! end if 1039 ! end if 1040 ibd=ibd+inc 1041 end do ! end loop on iblockbd 1042 1043 ABI_FREE(deltawf1) 1044 ABI_FREE(deltawf2) 1045 ABI_FREE(wf1) 1046 ABI_FREE(cwavef) 1047 ABI_FREE(cwavef1) 1048 ABI_FREE(eig) 1049 ABI_FREE(evec) 1050 ABI_FREE(unm) 1051 if(usepaw==1) then 1052 call pawcprj_free(cprj_k) 1053 ABI_FREE(cprj_k) 1054 call pawcprj_free(cprj_k1) 1055 ABI_FREE(cprj_k1) 1056 call pawcprj_free(cprj_k2) 1057 ABI_FREE(cprj_k2) 1058 end if 1059 1060 ibg=ibg+my_nspinor*nband_k 1061 icg=icg+my_nspinor*nband_k*npw_k 1062 1063 ! End big k point loop 1064 end do 1065 ! End loop over spins 1066 end do 1067 1068 if(usepaw==1) then 1069 call pawcprj_free(cprj) 1070 ABI_FREE(cprj) 1071 end if 1072 if (nprocband>1) then 1073 ABI_FREE(npw_block) 1074 ABI_FREE(npw_disp) 1075 ABI_FREE(bufsize) 1076 ABI_FREE(bufdisp) 1077 ABI_FREE(bufsize_wf) 1078 ABI_FREE(bufdisp_wf) 1079 end if 1080 1081 end if ! istep>=2 1082 1083 if (usepaw==1) then 1084 ABI_FREE(dimcprj) 1085 end if 1086 1087 end subroutine extrapwf
ABINIT/extrapwf_biortho [ Functions ]
NAME
extrapwf_biortho
FUNCTION
Extrapolate wavefunctions for new ionic positions from values of wavefunctions of previous SCF cycle. Use biorthogonal algorithm proposed XG
INPUTS
atindx1(dtset%natom)=index table for atoms, inverse of atindx dtset <type(dataset_type)>=all input variables in this dataset istep=number of call the routine mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of cprj array mpi_enreg=information about MPI parallelization nattyp(dtset%ntypat)=number of atoms of each type in cell. npwarr(nkpt)=number of planewaves in basis at this k point pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
SIDE EFFECTS
cg(2,mcg)= plane wave wavefunction coefficient Value from previous SCF cycle is input and stored in some form Extrapolated value is output cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors Value from previous SCF cycle is input and stored in some form Extrapolated value is output scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles
SOURCE
1123 subroutine extrapwf_biortho(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,& 1124 & nattyp,npwarr,pawtab,scf_history_wf) 1125 1126 !use m_scf_history 1127 use m_cgcprj, only : dotprod_set_cgcprj,cgcprj_cholesky,lincom_cgcprj 1128 1129 !Arguments ------------------------------------ 1130 !scalars 1131 integer,intent(in) :: istep,mcg,mcprj 1132 type(MPI_type),intent(in) :: mpi_enreg 1133 type(dataset_type),intent(in) :: dtset 1134 type(scf_history_type),intent(inout) :: scf_history_wf 1135 !arrays 1136 integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat) 1137 integer,intent(in) :: npwarr(dtset%nkpt) 1138 real(dp), intent(inout) :: cg(2,mcg) 1139 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj) 1140 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 1141 1142 !Local variables------------------------------- 1143 !scalars 1144 integer :: hermitian 1145 integer :: ibdmix,ibg,ibg_hist,icg,icg_hist !,iband 1146 integer :: ierr,ikpt,indh,ind1,ind2,ind1new,inplace 1147 integer :: isppol,istwf_k,kk,me_distrb,mband,my_nspinor,mcprj_k 1148 integer :: nband_k,nbdmix,nbdmax,npw_k,ntypat 1149 integer :: spaceComm_band,usepaw 1150 real(dp) :: alpha,beta !,dotr,doti 1151 1152 !arrays 1153 integer,allocatable :: ipiv(:),dimcprj(:) 1154 real(dp),allocatable ::psi_ortho(:,:),mmn(:,:,:) 1155 real(dp),allocatable :: smn(:,:,:) 1156 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:) 1157 1158 ! ************************************************************************* 1159 1160 if (istep==0) return 1161 1162 ntypat=dtset%ntypat 1163 usepaw=dtset%usepaw 1164 mband=dtset%mband 1165 nbdmax=dtset%mband 1166 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 1167 me_distrb=mpi_enreg%me_kpt 1168 spaceComm_band=xmpi_comm_self 1169 1170 !scf_history_wf%alpha contains dtset%wfmix 1171 alpha=scf_history_wf%alpha 1172 beta=scf_history_wf%beta 1173 ind1=scf_history_wf%hindex(1) 1174 ind2=scf_history_wf%hindex(2) 1175 ind1new=scf_history_wf%hindex(3) 1176 icg=0 1177 icg_hist=0 1178 ibg=0 1179 ibg_hist=0 1180 1181 !Useful array 1182 ABI_MALLOC(dimcprj,(dtset%natom)) 1183 if (usepaw==1) then 1184 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O') 1185 end if 1186 1187 if(istep==1)then 1188 do indh=1,scf_history_wf%history_size 1189 call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj) 1190 end do 1191 end if 1192 1193 mcprj_k=my_nspinor*nbdmax 1194 ABI_MALLOC(cprj_k,(dtset%natom,mcprj_k)) 1195 ABI_MALLOC(cprj_kh,(dtset%natom,mcprj_k)) 1196 1197 if(usepaw==1) then 1198 call pawcprj_alloc(cprj_k,0,dimcprj) 1199 call pawcprj_alloc(cprj_kh,0,dimcprj) 1200 end if 1201 ABI_MALLOC(smn,(2,nbdmax,nbdmax)) 1202 ABI_MALLOC(mmn,(2,nbdmax,nbdmax)) 1203 1204 !Explanation for the index for the wavefunction stored in scf_history_wf 1205 !The reference is the cg+cprj output after the wf optimization at istep 1. 1206 !For wavefunction mixing for molecular dynamics, we use the same mixing as for the density in extraprho. To keep the same indexes, 1207 ! we choose to take indh=3 for the reference. 1208 1209 !First step 1210 if (istep==1) then 1211 1212 indh=3 ! This input wavefunction is the reference 1213 1214 ! LOOP OVER SPINS 1215 do isppol=1,dtset%nsppol 1216 1217 ! BIG FAT k POINT LOOP 1218 do ikpt=1,dtset%nkpt 1219 1220 ! Select k point to be treated by this proc 1221 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1222 nbdmix=min(nband_k,nbdmax) 1223 1224 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 1225 1226 npw_k=npwarr(ikpt) 1227 1228 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 1229 1230 if(usepaw==1) then 1231 ! scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix) 1232 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,mband,& 1233 & dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,& 1234 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 1235 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,0,isppol,& 1236 & mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 1237 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1238 end if 1239 1240 ! Update the counters 1241 ibg=ibg+my_nspinor*nband_k 1242 ibg_hist=ibg_hist+my_nspinor*nbdmix 1243 icg=icg+my_nspinor*nband_k*npw_k 1244 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 1245 1246 end do 1247 end do 1248 1249 else 1250 ! From istep==2 1251 if (istep==2) ind1=3 1252 if (istep==3) ind2=3 1253 ! biorthogonalization 1254 indh=3 ! This input wavefunction is the reference 1255 1256 ! LOOP OVER SPINS 1257 do isppol=1,dtset%nsppol 1258 1259 ! BIG FAT k POINT LOOP 1260 do ikpt=1,dtset%nkpt 1261 1262 ! Select k point to be treated by this proc 1263 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1264 nbdmix=min(nband_k,nbdmax) 1265 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 1266 istwf_k=dtset%istwfk(ikpt) 1267 npw_k=npwarr(ikpt) 1268 ABI_MALLOC(psi_ortho,(2,npw_k*my_nspinor*nbdmix)) 1269 psi_ortho=zero 1270 ! Biorthogonalization 1271 1272 if(usepaw==1) then 1273 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,mband,& 1274 & dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,& 1275 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 1276 call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,0,isppol,& 1277 & mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,my_nspinor,dtset%nsppol,0,& 1278 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 1279 end if !end usepaw=1 1280 1281 hermitian=0 1282 1283 call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,& 1284 & 0,0,icg_hist,icg,ikpt,isppol,istwf_k,nbdmax,mcg,mcg,mcprj_k,mcprj_k,dtset%mkmem,& 1285 & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn(:,1:nbdmix,1:nbdmix),usepaw) 1286 1287 ! Invert S matrix, that is NOT hermitian. 1288 ! Calculate M=S^-1 1289 mmn=zero 1290 do kk=1,nbdmix 1291 mmn(1,kk,kk)=one 1292 end do 1293 1294 ABI_MALLOC(ipiv,(nbdmix)) 1295 ! The smn is destroyed by the following inverse call 1296 call zgesv(nbdmix,nbdmix,smn,nbdmax,ipiv,mmn,nbdmax,ierr) 1297 ABI_FREE(ipiv) 1298 !DEBUG 1299 if(ierr/=0)then 1300 ABI_ERROR(' The call to cgesv general inversion routine failed') 1301 end if 1302 !ENDDEBUG 1303 1304 ! The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place 1305 inplace=0 1306 call lincom_cgcprj(mmn(:,1:nbdmix,1:nbdmix),cg,cprj_k,dimcprj,& 1307 & icg,inplace,mcg,mcprj_k,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,& 1308 & cgout=psi_ortho,cprjout=cprj_kh,icgout=0) 1309 1310 !!!TEST 1311 ! if (usepaw==0) then 1312 ! do iband=1,nband_k 1313 ! call dotprod_g(dotr,doti,istwf_k,npw_k,2,scf_history_wf%cg(:,icg+1+my_nspinor*npw_k:icg+2*my_nspinor*npw_k,indh),& 1314 !& psi_ortho(:,1+(iband-1)*my_nspinor*npw_k:iband*my_nspinor*npw_k),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1315 ! write(80+mpi_enreg%me,*) dotr,doti 1316 ! flush(80+mpi_enreg%me) 1317 ! end do 1318 ! else 1319 ! hermitian=0 1320 ! call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),psi_ortho,scf_history_wf%cprj(:,:,indh)& 1321 !& ,cprj_kh,dimcprj,hermitian,0,0,icg_hist,icg_hist,ikpt,isppol,istwf_k,nbdmax,mcg,mcg,mcprj_k,mcprj_k,dtset%mkmem,& 1322 !& mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn(:,1:nbdmix,1:nbdmix),usepaw) 1323 ! write(90+mpi_enreg%me,*) smn 1324 ! flush(90+mpi_enreg%me) 1325 ! end if 1326 !!!TEST 1327 ! The biorthogonalised set of wavefunctions is now stored at the proper place 1328 1329 ! cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero 1330 1331 ! psi(t+dt) <- psi(t) + alpha.psi(t) 1332 cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=(one+alpha)*psi_ortho(:,1:my_nspinor*npw_k*nbdmix) 1333 if(usepaw==1) then 1334 do ibdmix=1,nbdmix 1335 call pawcprj_axpby(one+alpha,zero,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix)) 1336 end do ! end loop on ibdmix 1337 end if 1338 ! psi(t+dt) <- -alpha.psi(t-dt) + beta.psi(t-dt) 1339 if (abs(beta-alpha)>tol14.and.ind1>0) then 1340 cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)+& 1341 & (beta-alpha)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind1) 1342 if(usepaw==1) then 1343 do ibdmix=1,nbdmix 1344 call pawcprj_axpby(beta-alpha,one,scf_history_wf%cprj(:,ibg_hist+ibdmix:ibg_hist+ibdmix,ind1),& 1345 & cprj_k(:,ibdmix:ibdmix)) 1346 end do ! end loop on ibdmix 1347 end if 1348 end if 1349 1350 ! psi(t+dt) <- -beta.psi(t-2dt) 1351 if (abs(beta)>tol14.and.ind2>0) then 1352 cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)& 1353 & -beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind2) 1354 if(usepaw==1) then 1355 do ibdmix=1,nbdmix 1356 call pawcprj_axpby(-beta,one,scf_history_wf%cprj(:,ibg_hist+ibdmix:ibg_hist+ibdmix,ind2),cprj_k(:,ibdmix:ibdmix)) 1357 end do ! end loop on ibdmix 1358 end if 1359 end if 1360 1361 ! Store psi(t) in history 1362 scf_history_wf%cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix,ind1new)=psi_ortho(:,1:my_nspinor*npw_k*nbdmix) 1363 if(usepaw==1) then 1364 call pawcprj_put(atindx1,cprj_kh,scf_history_wf%cprj(:,:,ind1new),dtset%natom,1,ibg_hist,ikpt,0,isppol,& 1365 & nbdmix,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 1366 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1367 end if 1368 1369 ! Back to usual orthonormalization for the cg and cprj_k 1370 call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,mcprj_k,dtset%mkmem,& 1371 & mpi_enreg,dtset%natom,nattyp,nbdmax,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw) 1372 1373 ! Need to transfer cprj_k to cprj 1374 if(usepaw==1) then 1375 call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,0,isppol,& 1376 & mband,dtset%mkmem,dtset%natom,nbdmax,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 1377 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 1378 end if 1379 1380 ibg=ibg+my_nspinor*nband_k 1381 ibg_hist=ibg_hist+my_nspinor*nbdmix 1382 icg=icg+my_nspinor*nband_k*npw_k 1383 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 1384 ABI_FREE(psi_ortho) 1385 ! End big k point loop 1386 end do 1387 ! End loop over spins 1388 end do 1389 1390 1391 end if !istep>1 1392 1393 1394 if(usepaw==1) then 1395 call pawcprj_free(cprj_k) 1396 call pawcprj_free(cprj_kh) 1397 end if 1398 ABI_FREE(cprj_k) 1399 ABI_FREE(cprj_kh) 1400 ABI_FREE(dimcprj) 1401 ABI_FREE(mmn) 1402 ABI_FREE(smn) 1403 1404 1405 1406 end subroutine extrapwf_biortho
ABINIT/m_extraprho [ Modules ]
NAME
m_extraprho
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_extraprho 22 23 use defs_basis 24 use m_abicore 25 use m_scf_history 26 use m_errors 27 use m_xmpi 28 use m_cgtools 29 use m_dtset 30 31 use defs_datatypes, only : pseudopotential_type 32 use defs_abitypes, only : MPI_type 33 use m_atomdata, only : atom_length 34 use m_numeric_tools, only : hermit 35 use m_geometry, only : metric 36 use m_kg, only : getph 37 use m_jellium, only : jellium 38 use m_atm2fft, only : atm2fft 39 use m_pawtab, only : pawtab_type 40 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_inquire_dim, pawrhoij_filter 41 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, & 42 pawcprj_free, pawcprj_zaxpby,pawcprj_axpby, pawcprj_put, pawcprj_getdim 43 use m_mpinfo, only : proc_distrb_cycle 44 use m_cgprj, only : ctocprj 45 46 implicit none 47 48 private