TABLE OF CONTENTS
ABINIT/scfopt [ Functions ]
NAME
scfopt
FUNCTION
Compute the next vtrial of the SCF cycle. Possible algorithms are : simple mixing, Anderson (order 1 or 2), Pulay
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG,GMR) 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
cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX iscf= 2 => simple mixing = 3,4 => Anderson mixing = 7 => Pulay mixing istep= number of the step in the SCF cycle mpicomm=the mpi communicator used for the summation comm_atom=the mpi communicator over atoms ; PAW only (optional argument) mpi_summarize=set it to .true. if parallelisation is done over FFT nfft=(effective) number of FFT grid points (for this processor) npawmix=-PAW only- number of spherical part elements to be mixed nspden=number of spin-density components n_fftgr=third dimension of the array f_fftgr n_index=dimension for indices of potential/density (see ivrespc, i_vtrial...) opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential 1 vtrial (and also f_fftgr) actually contains the trial density pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part usepaw= 0 for non paw calculation; =1 for paw calculation
OUTPUT
(see side effects)
SIDE EFFECTS
vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave the input preconditioned residual potential at output, it is the new trial potential . f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid : The input vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(1)). The old vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(2)). The input preconditioned residual potential is in f_fftgr(:,:,i_vrespc(1)) Two input old preconditioned residual potentials in f_fftgr(:,:,i_vrespc(2)) and f_fftgr(:,:,i_vrespc(3)) Before output a permutation of i_vrespc(1), i_vrespc(2) and i_vrespc(3) occurs, without actually copying all the data (change of pointer). i_vrespc(n_index)=index of the preconditioned residual potentials (present and past) in the array f_fftgr i_vtrial(n_index) =indices of the potential (present and past) in the array f_fftgr ==== if usepaw==1 f_paw(npawmix,n_fftgr*mffmem*usepaw)=different functions used for PAW (same as f_fftgr but for spherical part) vpaw(npawmix*usepaw)=at input, the aug. occupancies (rhoij) that gave the input preconditioned residual potential at output, it is the new aug. occupancies.
PARENTS
m_ab7_mixing
CHILDREN
dgetrf,dgetri,dotprodm_v,sqnormm_v,wrtout,xmpi_sum
SOURCE
68 #if defined HAVE_CONFIG_H 69 #include "config.h" 70 #endif 71 72 #include "abi_common.h" 73 74 subroutine scfopt(cplex,f_fftgr,f_paw,iscf,istep,i_vrespc,i_vtrial,& 75 & mpicomm,mpi_summarize,nfft,npawmix,nspden,n_fftgr,& 76 & n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial,errid,errmess, & 77 & comm_atom) ! optional 78 79 use defs_basis 80 use m_linalg_interfaces 81 use m_profiling_abi 82 use m_xmpi 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 'scfopt' 88 use interfaces_14_hidewrite 89 use interfaces_56_mixing, except_this_one => scfopt 90 !End of the abilint section 91 92 implicit none 93 94 !Arguments ------------------------------------ 95 !scalars 96 integer,intent(in) :: cplex,iscf,istep,n_fftgr,n_index,nfft 97 integer,intent(in) :: npawmix,nspden,opt_denpot,pawoptmix,usepaw,mpicomm 98 integer, intent(in),optional :: comm_atom 99 integer,intent(out) :: errid 100 character(len = 500), intent(out) :: errmess 101 logical, intent(in) :: mpi_summarize 102 real(dp), intent(out) :: vresid 103 !arrays 104 integer,intent(inout) :: i_vrespc(n_index),i_vtrial(n_index) 105 real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr) 106 real(dp),intent(inout) :: f_paw(npawmix,n_fftgr*usepaw),vpaw(npawmix*usepaw) 107 real(dp),intent(inout) :: vtrial(cplex*nfft,nspden) 108 !Local variables------------------------------- 109 !scalars 110 integer,parameter :: npulaymax=50 111 integer :: i_vstore,ierr,ifft,ii,index,isp,jj,comm_atom_,niter,npulay,tmp 112 real(dp),save :: prod_resid_old,resid_old,resid_old2 113 real(dp) :: aa1,aa2,bb,cc1,cc2,current,det,lambda,lambda2,resid_best 114 character(len=500) :: message 115 !arrays 116 integer,allocatable :: ipiv(:) 117 real(dp),save :: amat(npulaymax+1,npulaymax+1) 118 real(dp) :: mpibuff(2),prod_resid(1),prod_resid2(1),resid_new(1) 119 real(dp),allocatable :: alpha(:),amatinv(:,:),amat_paw(:),rwork(:) 120 121 ! ************************************************************************* 122 123 !DEBUG 124 !write(std_out,*)' scfopt : enter ; istep,iscf ',istep,iscf 125 !ENDDEBUG 126 127 errid = AB7_NO_ERROR 128 129 comm_atom_=xmpi_comm_self; if(present(comm_atom)) comm_atom_=comm_atom 130 131 i_vstore=i_vtrial(1) 132 if (iscf==4) i_vstore=i_vtrial(2) 133 if (iscf==7) then 134 if (modulo(n_fftgr, 2) == 0 ) then 135 npulay=(n_fftgr-2)/2 136 else 137 npulay=(n_fftgr-1)/2 138 end if 139 i_vstore=i_vtrial(npulay) 140 else 141 npulay=0 142 end if 143 144 !Compute the new residual resid_new, from f_fftgr/f_paw(:,:,i_vrespc(1)) 145 call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr) 146 if (usepaw==1.and.pawoptmix==1) then 147 do index=1,npawmix 148 resid_new(1)=resid_new(1)+f_paw(index,i_vrespc(1))**2 149 end do 150 call xmpi_sum(resid_new(1),comm_atom_,ierr) 151 end if 152 vresid = resid_new(1) 153 154 !_______________________________________________________________ 155 !Here use only the preconditioning, or initialize the other algorithms 156 157 if(istep==1 .or. iscf==2)then 158 159 write(message,'(2a)') ch10,'Simple mixing update:' 160 call wrtout(std_out,message,'COLL') 161 162 write(message,*)' residual square of the potential :',resid_new(1) 163 call wrtout(std_out,message,'COLL') 164 165 ! Store information for later use 166 if (iscf==3.or.iscf==4) resid_old=resid_new(1) 167 if (iscf==7) then 168 amat(:,:)=zero 169 amat(1,1)=resid_new(1) 170 end if 171 172 ! Compute new vtrial (and new rhoij if PAW) 173 if (iscf/=2) f_fftgr(:,:,i_vstore)=vtrial(:,:) 174 vtrial(:,:)=vtrial(:,:)+f_fftgr(:,:,i_vrespc(1)) 175 if (usepaw==1) then 176 if (iscf/=2) f_paw(:,i_vstore)=vpaw(:) 177 vpaw(:)=vpaw(:)+f_paw(:,i_vrespc(1)) 178 end if 179 180 ! _______________________________________________________________ 181 ! Here Anderson algorithm using one previous iteration 182 else if((istep==2 .or. iscf==3).and.iscf/=7)then 183 184 write(message,'(2a)') ch10,'Anderson update:' 185 call wrtout(std_out,message,'COLL') 186 187 write(message,*)' residual square of the potential: ',resid_new(1) 188 call wrtout(std_out,message,'COLL') 189 190 ! Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2)) 191 call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,& 192 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 193 if (usepaw==1.and.pawoptmix==1) then 194 do index=1,npawmix 195 prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2)) 196 end do 197 call xmpi_sum(prod_resid(1),comm_atom_,ierr) 198 end if 199 200 ! Compute mixing factor 201 lambda=(resid_new(1)-prod_resid(1))/(resid_new(1)+resid_old-2*prod_resid(1)) 202 write(message,*)' mixing of old trial potential :',lambda 203 call wrtout(std_out,message,'COLL') 204 205 ! Evaluate best residual square on the line 206 resid_best=(1.0_dp-lambda)*(1.0_dp-lambda)*resid_new(1)& 207 & +(1.0_dp-lambda)*lambda *2*prod_resid(1)& 208 & +lambda *lambda *resid_old 209 write(message,*)' predicted best residual square on the line: ',resid_best 210 call wrtout(std_out,message,'COLL') 211 212 ! Store information for later use 213 if (iscf==4) then 214 prod_resid_old=prod_resid(1) 215 resid_old2=resid_old 216 end if 217 resid_old=resid_new(1) 218 219 ! Save latest trial potential and compute new trial potential 220 do isp=1,nspden 221 do ifft=1,cplex*nfft 222 current=vtrial(ifft,isp) 223 vtrial(ifft,isp)=(one-lambda)*(current +f_fftgr(ifft,isp,i_vrespc(1)))& 224 & +lambda *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2))) 225 f_fftgr(ifft,isp,i_vstore)=current 226 end do 227 end do 228 229 ! PAW: save latest rhoij and compute new rhoij 230 do index=1,npawmix 231 current=vpaw(index) 232 vpaw(index)=(one-lambda)*(current +f_paw(index,i_vrespc(1)))& 233 & +lambda *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2))) 234 f_paw(index,i_vstore)=current 235 end do 236 237 ! _______________________________________________________________ 238 ! Here Anderson algorithm using two previous iterations 239 else if(iscf==4.and.iscf/=7)then 240 241 write(message,'(2a)') ch10,'Anderson (order 2) update:' 242 call wrtout(std_out,message,'COLL') 243 244 write(message,*)' residual square of the potential :',resid_new(1) 245 call wrtout(std_out,message,'COLL') 246 247 ! Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2)) 248 call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,& 249 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 250 if (usepaw==1.and.pawoptmix==1) then 251 do index=1,npawmix 252 prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2)) 253 end do 254 end if 255 256 ! Compute prod_resid2 from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(3)) 257 call dotprodm_v(cplex,1,prod_resid2,i_vrespc(1),i_vrespc(3),mpicomm,mpi_summarize,1,1,& 258 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 259 if (usepaw==1.and.pawoptmix==1) then 260 do index=1,npawmix 261 prod_resid2(1)=prod_resid2(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(3)) 262 end do 263 ! MPI reduction 264 mpibuff(1)=prod_resid(1);mpibuff(2)=prod_resid2(1) 265 call xmpi_sum(mpibuff,comm_atom_,ierr) 266 prod_resid(1)=mpibuff(1);prod_resid2(1)=mpibuff(2) 267 end if 268 269 ! Compute mixing factors 270 aa1=resid_new(1)+resid_old -two*prod_resid (1) 271 aa2=resid_new(1)+resid_old2-two*prod_resid2(1) 272 bb =resid_new(1)+prod_resid_old-prod_resid(1)-prod_resid2(1) 273 cc1=resid_new(1)-prod_resid (1) 274 cc2=resid_new(1)-prod_resid2(1) 275 det=aa1*aa2-bb*bb 276 lambda =(aa2*cc1-bb*cc2)/det 277 lambda2=(aa1*cc2-bb*cc1)/det 278 write(message,*)' mixing of old trial potentials :',lambda,lambda2 279 call wrtout(std_out,message,'COLL') 280 281 ! Store information for later use 282 prod_resid_old=prod_resid(1) 283 resid_old2=resid_old 284 resid_old=resid_new(1) 285 286 ! Save latest trial potential and compute new trial potential 287 do isp=1,nspden 288 do ifft=1,cplex*nfft 289 current=vtrial(ifft,isp) 290 vtrial(ifft,isp)=& 291 & (one-lambda-lambda2)*(current +f_fftgr(ifft,isp,i_vrespc(1)))& 292 & +lambda *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))& 293 & +lambda2 *(f_fftgr(ifft,isp,i_vtrial(2))+f_fftgr(ifft,isp,i_vrespc(3))) 294 f_fftgr(ifft,isp,i_vstore)=current 295 end do 296 end do 297 298 ! PAW: save latest rhoij and compute new rhoij 299 do index=1,npawmix 300 current=vpaw(index) 301 vpaw(index)=& 302 & (one-lambda-lambda2)*(current +f_paw(index,i_vrespc(1)))& 303 & +lambda *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))& 304 & +lambda2 *(f_paw(index,i_vtrial(2))+f_paw(index,i_vrespc(3))) 305 f_paw(index,i_vstore)=current 306 end do 307 308 ! _______________________________________________________________ 309 ! Here Pulay algorithm 310 else if(iscf==7)then 311 312 niter=min(istep,npulay+1) 313 314 write(message,'(2a,i2,a)') ch10,' Pulay update with ',niter-1,' previous iterations:' 315 call wrtout(std_out,message,'COLL') 316 317 if (npulay>npulaymax) then 318 errid = AB7_ERROR_MIXING_CONVERGENCE 319 write(errmess, '(4a)' ) ch10,& 320 & ' scfopt : ERROR - ',ch10,& 321 & ' Too much iterations required for Pulay algorithm (<50) !' 322 return 323 end if 324 325 ! Compute "A" matrix 326 if (istep>npulay+1) then 327 do jj=1,niter-1 328 do ii=1,niter-1 329 amat(ii,jj)=amat(ii+1,jj+1) 330 end do 331 end do 332 end if 333 if (usepaw==1.and.pawoptmix==1) then 334 ABI_ALLOCATE(amat_paw,(niter)) 335 amat_paw(:)=zero 336 do ii=1,niter 337 do index=1,npawmix 338 amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(1+niter-ii)) 339 end do 340 end do 341 call xmpi_sum(amat_paw,comm_atom_,ierr) 342 end if 343 do ii=1,niter 344 call dotprodm_v(cplex,1,amat(ii,niter),i_vrespc(1),i_vrespc(1+niter-ii),mpicomm,mpi_summarize,1,1,& 345 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 346 if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii) 347 if (ii<niter) amat(niter,ii)=amat(ii,niter) 348 end do 349 if (usepaw==1.and.pawoptmix==1)then 350 ABI_DEALLOCATE(amat_paw) 351 end if 352 353 ! Invert "A" matrix 354 ABI_ALLOCATE(amatinv,(niter,niter)) 355 amatinv(1:niter,1:niter)=amat(1:niter,1:niter) 356 ABI_ALLOCATE(ipiv,(niter)) 357 ABI_ALLOCATE(rwork,(niter)) 358 call dgetrf(niter,niter,amatinv,niter,ipiv,ierr) 359 call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr) 360 ABI_DEALLOCATE(ipiv) 361 ABI_DEALLOCATE(rwork) 362 363 ! Compute "alpha" factors 364 ABI_ALLOCATE(alpha,(niter)) 365 alpha=zero 366 det=zero 367 do ii=1,niter 368 do jj=1,niter 369 alpha(ii)=alpha(ii)+amatinv(jj,ii) 370 det=det+amatinv(jj,ii) 371 end do 372 end do 373 alpha(:)=alpha(:)/det 374 ABI_DEALLOCATE(amatinv) 375 write(message,'(a,5(1x,g10.3))')' mixing of old trial potential : alpha(m:m-4)=',(alpha(ii),ii=niter,max(1,niter-4),-1) 376 call wrtout(std_out,message,'COLL') 377 378 ! Save latest trial potential and compute new trial potential 379 do isp=1,nspden 380 do ifft=1,cplex*nfft 381 current=vtrial(ifft,isp) 382 vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,i_vrespc(1))) 383 do ii=niter-1,1,-1 384 vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) & 385 & *(f_fftgr(ifft,isp,i_vtrial(niter-ii))+f_fftgr(ifft,isp,i_vrespc(1+niter-ii))) 386 end do 387 f_fftgr(ifft,isp,i_vstore)=current 388 end do 389 end do 390 391 ! PAW: save latest rhoij and compute new rhoij 392 do index=1,npawmix 393 current=vpaw(index) 394 vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1))) 395 do ii=niter-1,1,-1 396 vpaw(index)=vpaw(index)+alpha(ii) & 397 & *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii))) 398 end do 399 f_paw(index,i_vstore)=current 400 end do 401 402 ABI_DEALLOCATE(alpha) 403 ! _______________________________________________________________ 404 ! End of choice of optimization method 405 end if 406 407 !Permute potential indices 408 if (iscf==3) then 409 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 410 else if (iscf==4) then 411 tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 412 tmp=i_vtrial(2) ; i_vtrial(2)=i_vtrial(1) ; i_vtrial(1)=tmp 413 else if (iscf==7) then 414 tmp=i_vtrial( npulay) 415 do ii= npulay,2,-1 416 i_vtrial(ii)=i_vtrial(ii-1) 417 end do 418 i_vtrial(1)=tmp 419 tmp=i_vrespc(1+npulay) 420 do ii=1+npulay,2,-1 421 i_vrespc(ii)=i_vrespc(ii-1) 422 end do 423 i_vrespc(1)=tmp 424 end if 425 426 end subroutine scfopt