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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) 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 atindx1(natom)=index table for atoms, inverse of atindx cg(2,mcg)= plane wave wavefunction coefficient 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 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
PARENTS
scfcv
CHILDREN
atm2fft,extrapwf,jellium,pawrhoij_alloc,wrtout
SOURCE
72 #if defined HAVE_CONFIG_H 73 #include "config.h" 74 #endif 75 76 #include "abi_common.h" 77 78 subroutine extraprho(atindx,atindx1,cg,dtset,gmet,gprimd,gsqcut,istep,& 79 & kg,mcg,mgfft,mpi_enreg,mqgrid,my_natom,nattyp,nfft,ngfft,npwarr,ntypat,pawrhoij,& 80 & pawtab,ph1d,psps,qgrid,rhor,rprimd,scf_history,ucvol,usepaw,& 81 & xred_new,xred_old,ylm,zion,znucl) 82 83 use defs_basis 84 use defs_datatypes 85 use defs_abitypes 86 use m_scf_history 87 use m_errors 88 89 use m_atomdata, only : atom_length 90 use m_pawtab, only : pawtab_type 91 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc 92 93 !This section has been created automatically by the script Abilint (TD). 94 !Do not modify the following lines by hand. 95 #undef ABI_FUNC 96 #define ABI_FUNC 'extraprho' 97 use interfaces_14_hidewrite 98 use interfaces_64_psp 99 use interfaces_67_common, except_this_one => extraprho 100 !End of the abilint section 101 102 implicit none 103 104 !Arguments ------------------------------------ 105 !scalars 106 integer,intent(in) :: istep,mcg,mgfft,my_natom,mqgrid,nfft,ntypat,usepaw 107 real(dp),intent(in) :: gsqcut,ucvol 108 type(MPI_type),intent(in) :: mpi_enreg 109 type(dataset_type),intent(in) :: dtset 110 type(scf_history_type),intent(inout) :: scf_history 111 type(pseudopotential_type),intent(in) :: psps 112 !arrays 113 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem) 114 integer,intent(in) :: nattyp(ntypat),ngfft(18),npwarr(dtset%nkpt) 115 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom) 116 real(dp),intent(in) :: qgrid(mqgrid),rprimd(3,3) 117 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 118 real(dp),intent(in) :: zion(ntypat),znucl(ntypat) 119 real(dp), intent(inout) :: cg(2,mcg) 120 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),xred_new(3,dtset%natom) 121 real(dp),intent(in) :: xred_old(3,dtset%natom) 122 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw) 123 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 124 125 !Local variables------------------------------- 126 !scalars 127 integer :: cplex,dplex,iatom,ii,ind1,ind1new,ind2,ind2new,irhoij,ispden,itypat,jrhoij,klmn 128 integer :: lmn2_size,nselect,nspden_rhoij,optatm,optdyfr,opteltfr,optgr,option,optn,optn2 129 integer :: optstr,optv 130 real(dp) :: a11,a12,a22,a33,alpha,b1,b2,beta,detA,fact,ratio1,ratio2 131 logical :: hasmoved,usegauss 132 character(len=500) :: message 133 !arrays 134 integer :: dummy3(3) 135 real(dp) :: diff_t(3),diff_tmdt(3),diff_tpdt(3),dummy2(2) 136 real(dp) :: dummy_in(0) 137 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0) 138 real(dp) :: strn_dummy6(6),strv_dummy6(6) 139 real(dp),allocatable :: deltarho(:),gauss(:,:),rhoijtmp(:,:),work1(:) 140 real(dp),allocatable :: work2(:,:),work3(:,:),xred_tpdt(:,:) 141 142 ! ************************************************************************* 143 144 !--------------------------------------------------------------- 145 !----------- Inits 146 !--------------------------------------------------------------- 147 148 !History indexes 149 ind1=scf_history%hindex(1) 150 ind2=scf_history%hindex(2) 151 152 !Compatibility tests 153 if (ind1==0.and.ind2>0)then 154 MSG_BUG(' Incompatible history indexes !') 155 end if 156 157 !Rotated values of history indexes 158 if (ind1>0.and.ind2>0) then 159 ind1new=ind2;ind2new=ind1 160 else if (ind1>0.and.ind2==0) then 161 ind1new=3-ind1;ind2new=ind1 162 else if (ind1==0.and.ind2==0) then 163 ind1new=1;ind2new=0 164 end if 165 166 !Compute ionic positions at t+dt in red. coordinates 167 !Has to take the boundary conditions into account 168 ABI_ALLOCATE(xred_tpdt,(3,dtset%natom)) 169 do iatom=1,dtset%natom 170 xred_tpdt(1,iatom)=xred_old(1,iatom)+mod(xred_new(1,iatom)-xred_old(1,iatom)+1.5_dp,one)-half 171 xred_tpdt(2,iatom)=xred_old(2,iatom)+mod(xred_new(2,iatom)-xred_old(2,iatom)+1.5_dp,one)-half 172 xred_tpdt(3,iatom)=xred_old(3,iatom)+mod(xred_new(3,iatom)-xred_old(3,iatom)+1.5_dp,one)-half 173 end do 174 175 !--------------------------------------------------------------- 176 !----------- Compute Alpha and Beta 177 !----------- see (4) in Comp. Phys. Comm. 118 (1999), 31-33 178 !--------------------------------------------------------------- 179 180 !Compute a_ij matrix 181 a11=zero;a12=zero;a22=zero;a33=zero;b1=zero;b2=zero 182 diff_t=zero;diff_tmdt=zero;diff_tpdt=zero 183 do iatom=1,dtset%natom 184 185 diff_tpdt(1:3)=xred_tpdt(1:3,iatom)-xred_old(1:3,iatom) 186 if (ind1>0) then 187 diff_t(1:3)=scf_history%xreddiff(1:3,iatom,ind1) 188 if (ind2>0) diff_tmdt(1:3)=scf_history%xreddiff(1:3,iatom,ind2) 189 end if 190 do ii=1,3 191 a11=a11+diff_t(ii)**2 192 a22=a22+diff_tmdt(ii)**2 193 a33=a33+diff_tpdt(ii)**2 194 a12=a12+diff_t(ii) *diff_tmdt(ii) 195 b1 =b1 +diff_t(ii) *diff_tpdt(ii) 196 b2 =b2 +diff_tmdt(ii)*diff_tpdt(ii) 197 end do 198 199 ! Store reduced coordinates diffs in SCF history 200 scf_history%xreddiff(1:3,iatom,ind1new)=diff_tpdt(1:3) 201 202 end do 203 ABI_DEALLOCATE(xred_tpdt) 204 hasmoved=(a11>=tol10.or.a22>=tol10.or.a33>=tol10) 205 206 !Compute alpha and beta 207 alpha=zero;beta=zero 208 if (hasmoved.and.ind1>0) then 209 ratio1=one;if (abs(a33)>=tol10) ratio1=(a11+a33-two*b1)/a33 210 ratio2=one;if (abs(a33)>=tol10) ratio2=(a11+a33-two*b2)/a33 211 detA=a11*a22-a12**2 212 if (abs(a11)>=tol10.and.(abs(a22)<tol10.or.abs(detA)<tol10)) then 213 alpha=b1/a11 214 else if (abs(a22)>=tol10.and.(abs(a11)<tol10.or.abs(detA)<tol10)) then 215 beta=b2/a22 216 else if (abs(ratio1)+abs(ratio2)<tol6) then 217 if (ind2>0) then 218 alpha=two;beta=-one 219 else 220 alpha=one 221 end if 222 write(message,'(3a,f4.1,a,f4.1)')& 223 & 'Ionic positions lead to a collinear system !',ch10,& 224 & 'Mixing coeffs have been set to: alpha=',alpha,' beta=',beta 225 MSG_WARNING(message) 226 else if (abs(a11)>=tol10.and.abs(a22)>=tol10) then 227 alpha=(b1*a22-b2*a12)/detA 228 beta =(b2*a11-b1*a12)/detA 229 end if 230 end if 231 232 !--------------------------------------------------------------- 233 !----------- Contribution from delta_rho(t), delta_rho(t-dt) 234 !----------- and delta_rho(t-2dt) to predicted rho(t+dt) 235 !--------------------------------------------------------------- 236 237 !deltarho(t+dt) <- deltarho(t) + alpha.[deltarho(t)-deltarho(t-dt)] 238 !+ beta .[deltarho(t-dt)-deltarho(t-2dt)] 239 !Note: scf_history%deltarhor is updated at the same time 240 241 ABI_ALLOCATE(deltarho,(nfft)) 242 do ispden=1,dtset%nspden 243 244 if (ispden==1) then 245 deltarho(:)=rhor(:,ispden)-scf_history%atmrho_last(:) 246 else if (ispden==2.and.dtset%nspden==2) then 247 deltarho(:)=rhor(:,ispden)-half*scf_history%atmrho_last(:) 248 end if 249 250 251 ! rho(t+dt) <- deltarho(t) + alpha.deltarho(t) 252 if (dtset%nspden/=4.or.ispden==1) then 253 rhor(:,ispden)=(one+alpha)*deltarho(:) 254 else 255 rhor(:,ispden)=(one+alpha)*rhor(:,ispden) 256 end if 257 258 if (hasmoved) then 259 260 ! rho(t+dt) <- -alpha.deltarho(t-dt) + beta.deltarho(t-dt) 261 if (abs(beta-alpha)>tol14.and.ind1>0) then 262 rhor(:,ispden)=rhor(:,ispden)+(beta-alpha)*scf_history%deltarhor(:,ispden,ind1) 263 end if 264 265 ! rho(t+dt) <- -beta.deltarho(t-2dt) 266 if (abs(beta)>tol14.and.ind2>0) then 267 rhor(:,ispden)=rhor(:,ispden)-beta*scf_history%deltarhor(:,ispden,ind2) 268 end if 269 270 end if 271 272 ! Store deltarho(t) in history 273 if (dtset%nspden/=4.or.ispden==1) then 274 scf_history%deltarhor(:,ispden,ind1new)=deltarho(:) 275 else 276 scf_history%deltarhor(:,ispden,ind1new)=rhor(:,ispden) 277 end if 278 279 end do 280 281 ABI_DEALLOCATE(deltarho) 282 283 !--------------------------------------------------------------- 284 !----------- Contribution from rho_at(t+dt) to predicted rho(t+dt) 285 !--------------------------------------------------------------- 286 287 !Determine whether a gaussian atomic density has to be used or not 288 !MG: Note that there's a small inconsistency between initro and extraprho because in initrho 289 ! we use `use_gaussian(ntypat)`. 290 usegauss=.true. 291 if (usepaw==0) usegauss = any(.not.psps%nctab(1:ntypat)%has_tvale) 292 if (usepaw==1) usegauss=(minval(pawtab(1:ntypat)%has_tvale)==0) 293 if (usegauss) then 294 optn2=3 295 ABI_ALLOCATE(gauss,(2,ntypat)) 296 do itypat=1,ntypat 297 gauss(1,itypat)=zion(itypat) 298 gauss(2,itypat) = atom_length(dtset%densty(itypat,1),zion(itypat),znucl(itypat)) 299 end do 300 call wrtout(std_out," Extrapolating rho(t+dt) using gaussian functions as atomic densities", "COLL") 301 else 302 optn2=2 303 ABI_ALLOCATE(gauss,(2,0)) 304 call wrtout(std_out," Extrapolating rho(t+dt) using atomic densities taken from pseudos", "COLL") 305 end if 306 307 !Compute rho_at(t+dt) as sum of atomic densities 308 !Note: scf_history%atmrho_last is updated at the same time 309 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=0;optn=1 310 call atm2fft(atindx1,scf_history%atmrho_last,dummy_out1,dummy_out2,dummy_out3,& 311 & dummy_out4,gauss,gmet,gprimd,dummy_out5,dummy_out6,gsqcut,mgfft,mqgrid,dtset%natom,nattyp,& 312 & nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,& 313 & dummy3,dummy_in,strn_dummy6,strv_dummy6,ucvol,usepaw,dummy_in,dummy_in,dummy_in,dummy2,dummy_in,& 314 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 315 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 316 ABI_DEALLOCATE(gauss) 317 318 !Take eventually into account jellium slab 319 if (dtset%jellslab/=0) then 320 option=2 321 ABI_ALLOCATE(work1,(nfft)) 322 ABI_ALLOCATE(work2,(nfft,1)) 323 ABI_ALLOCATE(work3,(2,nfft)) 324 work2(:,1)=scf_history%atmrho_last(:) 325 call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,1,option,dtset%paral_kgb,& 326 & dtset%slabwsrad,work3,work2,rprimd,work1,dtset%slabzbeg,dtset%slabzend) 327 scf_history%atmrho_last(:)=work2(:,1) 328 ABI_DEALLOCATE(work1) 329 ABI_DEALLOCATE(work2) 330 ABI_DEALLOCATE(work3) 331 end if 332 333 !Add rho_at(t+dt) to rho(t+dt) 334 rhor(:,1)=rhor(:,1)+scf_history%atmrho_last(:) 335 if (dtset%nspden==2) rhor(:,2)=rhor(:,2)+half*scf_history%atmrho_last(:) 336 337 !--------------------------------------------------------------- 338 !----------- Extrapolation of PAW rhoij occupancy matrixes 339 !--------------------------------------------------------------- 340 341 if (usepaw==1) then 342 343 if (ind2==0) then 344 nspden_rhoij=dtset%nspden;if (dtset%pawspnorb>0.and.dtset%nspinor==2) nspden_rhoij=4 345 call pawrhoij_alloc(scf_history%pawrhoij(:,ind1new),dtset%pawcpxocc,nspden_rhoij,& 346 & dtset%nspinor,dtset%nsppol,dtset%typat,pawtab=pawtab,& 347 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 348 end if 349 350 do iatom=1,my_natom 351 352 nspden_rhoij=pawrhoij(iatom)%nspden 353 lmn2_size=pawrhoij(iatom)%lmn2_size 354 cplex=pawrhoij(iatom)%cplex;dplex=cplex-1 355 356 if (hasmoved) then 357 ABI_ALLOCATE(rhoijtmp,(cplex*lmn2_size,nspden_rhoij)) 358 rhoijtmp=zero 359 360 do ispden=1,nspden_rhoij 361 362 ! rhoij(t+dt) <- rhoij(t) + alpha.rhoij(t) 363 fact=one+alpha 364 jrhoij=1 365 do irhoij=1,pawrhoij(iatom)%nrhoijsel 366 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 367 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 368 & +fact*pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 369 jrhoij=jrhoij+cplex 370 end do 371 372 ! rhoij(t+dt) <- -alpha.rhoij(t-dt) + beta.rhoij(t-dt) 373 if (abs(beta-alpha)>tol14.and.ind1>0) then 374 fact=beta-alpha 375 jrhoij=1 376 do irhoij=1,scf_history%pawrhoij(iatom,ind1)%nrhoijsel 377 klmn=cplex*scf_history%pawrhoij(iatom,ind1)%rhoijselect(irhoij)-dplex 378 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 379 & +fact*scf_history%pawrhoij(iatom,ind1)%rhoijp(jrhoij:jrhoij+dplex,ispden) 380 jrhoij=jrhoij+cplex 381 end do 382 end if 383 384 ! rho(t+dt) <- -beta.rhoij(t-2dt) 385 if (abs(beta)>tol14.and.ind2>0) then 386 fact=-beta 387 jrhoij=1 388 do irhoij=1,scf_history%pawrhoij(iatom,ind2)%nrhoijsel 389 klmn=cplex*scf_history%pawrhoij(iatom,ind2)%rhoijselect(irhoij)-dplex 390 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 391 & +fact*scf_history%pawrhoij(iatom,ind2)%rhoijp(jrhoij:jrhoij+dplex,ispden) 392 jrhoij=jrhoij+cplex 393 end do 394 end if 395 396 end do !ispden 397 end if !hasmoved 398 399 ! Store rhoij(t) in history 400 ! (cannot use pawrhoij_copy here because update for single atom) 401 nselect=pawrhoij(iatom)%nrhoijsel 402 scf_history%pawrhoij(iatom,ind1new)%nrhoijsel=nselect 403 scf_history%pawrhoij(iatom,ind1new)%rhoijselect(1:nselect)=pawrhoij(iatom)%rhoijselect(1:nselect) 404 scf_history%pawrhoij(iatom,ind1new)%rhoijp(1:cplex*nselect,1:nspden_rhoij)= & 405 & pawrhoij(iatom)%rhoijp(1:cplex*nselect,1:nspden_rhoij) 406 407 ! Select non-zero values of rhoij(t+dt) 408 if (hasmoved) then 409 nselect=0 410 if (cplex==1) then 411 do klmn=1,lmn2_size 412 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 413 nselect=nselect+1 414 pawrhoij(iatom)%rhoijselect(nselect)=klmn 415 do ispden=1,nspden_rhoij 416 pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden) 417 end do 418 end if 419 end do 420 else 421 do klmn=1,lmn2_size 422 if (any(abs(rhoijtmp(2*klmn-1:2*klmn,:))>tol10)) then 423 nselect=nselect+1 424 pawrhoij(iatom)%rhoijselect(nselect)=klmn 425 do ispden=1,nspden_rhoij 426 pawrhoij(iatom)%rhoijp(2*nselect-1:2*nselect,ispden)=rhoijtmp(2*klmn-1:2*klmn,ispden) 427 end do 428 end if 429 end do 430 end if 431 pawrhoij(iatom)%nrhoijsel=nselect 432 ABI_DEALLOCATE(rhoijtmp) 433 end if 434 435 end do !iatom 436 end if !usepaw 437 438 scf_history%alpha=alpha 439 scf_history%beta=beta 440 441 !--------------------------------------------------------------- 442 !----------- End 443 !--------------------------------------------------------------- 444 445 if(scf_history%usecg==1) then 446 if (hasmoved) then 447 call extrapwf(atindx,atindx1,cg,dtset,istep,kg,mcg,mgfft,mpi_enreg,nattyp,& 448 & ngfft,npwarr,ntypat,pawtab,psps,rprimd,scf_history,usepaw,xred_old,ylm) 449 else 450 scf_history%cg(:,:,2)=zero 451 end if 452 453 end if 454 455 !Rotate history indexes 456 scf_history%hindex(1)=ind1new 457 scf_history%hindex(2)=ind2new 458 459 end subroutine extraprho