TABLE OF CONTENTS
ABINIT/pred_bfgs [ Functions ]
NAME
pred_bfgs
FUNCTION
Ionmov predictors (2 & 3) Broyden-Fletcher-Goldfarb-Shanno IONMOV 2: Given a starting point xred that is a vector of length 3*natom (reduced nuclei coordinates), and unit cell parameters (acell and rprimd) the Broyden-Fletcher-Goldfarb-Shanno minimization is performed on the total energy function, using its gradient (atomic forces and stresse) as calculated by the routine scfcv. Some atoms can be kept fixed, while the optimization of unit cell parameters is only performed if optcell/=0. The convergence requirement on the atomic forces, dtset%tolmxf, allows an early exit. Otherwise no more than dtset%ntime steps are performed. Returned quantities are xred, and eventually acell and rprim (new ones!). Could see Numerical Recipes (Fortran), 1986, page 307. IONMOV 3: Conduct structural optimization using the Broyden-Fletcher- Goldfarb-Shanno minimization (BFGS), modified to take into account the total energy as well as the gradients (as in usual BFGS). See the paper by Schlegel, J. Comp. Chem. 3, 214 (1982). Might be better than ionmov=2 for few degrees of freedom (less than 3 or 4)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, SE) 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
ab_mover <type(abimover)> : Datatype with all the information needed by the preditor itime : Index of the present iteration ntime : Maximal number of iterations ionmov : (2 or 3) Specific kind of BFGS zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
PARENTS
mover
CHILDREN
brdene,dgetrf,dgetri,fcart2fred,hessinit,hessupdt,hist2var,metric mkrdim,var2hist,xfh_recover_new,xfpack_f2vout,xfpack_vin2x,xfpack_x2vin
SOURCE
63 #if defined HAVE_CONFIG_H 64 #include "config.h" 65 #endif 66 67 #include "abi_common.h" 68 69 subroutine pred_bfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit) 70 71 use defs_basis 72 use m_profiling_abi 73 use m_abimover 74 use m_abihist 75 76 use m_geometry, only : mkrdim, fcart2fred, metric 77 use m_bfgs, only : hessinit, hessupdt, brdene 78 79 !This section has been created automatically by the script Abilint (TD). 80 !Do not modify the following lines by hand. 81 #undef ABI_FUNC 82 #define ABI_FUNC 'pred_bfgs' 83 use interfaces_45_geomoptim, except_this_one => pred_bfgs 84 !End of the abilint section 85 86 implicit none 87 88 !Arguments ------------------------------------ 89 !scalars 90 type(abimover),intent(in) :: ab_mover 91 type(ab_xfh_type),intent(inout) :: ab_xfh 92 type(abihist),intent(inout) :: hist 93 type(abiforstr),intent(in) :: forstr 94 integer,intent(in) :: itime 95 integer,intent(in) :: ionmov 96 integer,intent(in) :: iexit 97 logical,intent(in) :: zDEBUG 98 99 !Local variables------------------------------- 100 !scalars 101 integer :: ihist_prev,ndim,cycl_main 102 integer, parameter :: npul=0 103 integer :: ierr,ii,jj,kk,nitpul 104 real(dp),save :: ucvol0 105 real(dp) :: ucvol,det 106 real(dp) :: etotal,etotal_prev 107 real(dp) :: favg,alpha0 108 109 !arrays 110 integer,allocatable :: ipiv(:) 111 real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:) 112 real(dp),allocatable,save :: vout(:),vout_prev(:) 113 real(dp),allocatable,save ::vinres(:,:),vin1(:,:) 114 real(dp),allocatable :: amat(:,:),amatinv(:,:),alpha(:,:) 115 real(dp),allocatable :: rwork(:) 116 real(dp),save :: acell0(3) ! Initial acell 117 real(dp),save :: rprimd0(3,3) ! Initial rprimd 118 real(dp) :: acell(3) 119 real(dp) :: rprimd(3,3),rprim(3,3) 120 real(dp) :: gprimd(3,3) 121 real(dp) :: gmet(3,3) 122 real(dp) :: rmet(3,3) 123 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom) 124 real(dp) :: xred(3,ab_mover%natom),strten(6) 125 126 !*************************************************************************** 127 !Beginning of executable session 128 !*************************************************************************** 129 130 if(iexit/=0)then 131 if (allocated(vin)) then 132 ABI_DEALLOCATE(vin) 133 end if 134 if (allocated(vout)) then 135 ABI_DEALLOCATE(vout) 136 end if 137 if (allocated(vin_prev)) then 138 ABI_DEALLOCATE(vin_prev) 139 end if 140 if (allocated(vout_prev)) then 141 ABI_DEALLOCATE(vout_prev) 142 end if 143 if (allocated(vinres)) then 144 ABI_DEALLOCATE(vinres) 145 end if 146 if (allocated(vin1)) then 147 ABI_DEALLOCATE(vin1) 148 end if 149 if (allocated(hessin)) then 150 ABI_DEALLOCATE(hessin) 151 end if 152 return 153 end if 154 155 !write(std_out,*) 'bfgs 01' 156 !########################################################## 157 !### 01. Debugging and Verbose 158 159 if(zDEBUG)then 160 write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),& 161 & 'Debugging and Verbose for pred_bfgs',('-',kk=1,42) 162 write(std_out,*) 'ionmov: ',ionmov 163 write(std_out,*) 'itime: ',itime 164 end if 165 166 !write(std_out,*) 'bfgs 02' 167 !########################################################## 168 !### 02. Compute the dimension of vectors (ndim) 169 170 ndim=3*ab_mover%natom 171 if(ab_mover%optcell==1 .or.& 172 & ab_mover%optcell==4 .or.& 173 & ab_mover%optcell==5 .or.& 174 & ab_mover%optcell==6) ndim=ndim+1 175 if(ab_mover%optcell==2 .or.& 176 & ab_mover%optcell==3) ndim=ndim+6 177 if(ab_mover%optcell==7 .or.& 178 & ab_mover%optcell==8 .or.& 179 & ab_mover%optcell==9) ndim=ndim+3 180 181 if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim 182 183 !write(std_out,*) 'bfgs 03' 184 !########################################################## 185 !### 03. Allocate the vectors vin, vout and hessian matrix 186 187 !Notice that vin, vout, etc could be allocated 188 !From a previous dataset with a different ndim 189 if(itime==1)then 190 if (allocated(vin)) then 191 ABI_DEALLOCATE(vin) 192 end if 193 if (allocated(vout)) then 194 ABI_DEALLOCATE(vout) 195 end if 196 if (allocated(vin_prev)) then 197 ABI_DEALLOCATE(vin_prev) 198 end if 199 if (allocated(vout_prev)) then 200 ABI_DEALLOCATE(vout_prev) 201 end if 202 if (allocated(vinres)) then 203 ABI_DEALLOCATE(vinres) 204 end if 205 if (allocated(vin1)) then 206 ABI_DEALLOCATE(vin1) 207 end if 208 if (allocated(hessin)) then 209 ABI_DEALLOCATE(hessin) 210 end if 211 if(npul>1) then 212 ABI_ALLOCATE(vinres,(npul+1,ndim)) 213 ABI_ALLOCATE(vin1,(npul+1,ndim)) 214 end if 215 ABI_ALLOCATE(vin,(ndim)) 216 ABI_ALLOCATE(vout,(ndim)) 217 ABI_ALLOCATE(vin_prev,(ndim)) 218 ABI_ALLOCATE(vout_prev,(ndim)) 219 ABI_ALLOCATE(hessin,(ndim,ndim)) 220 end if 221 222 !write(std_out,*) 'bfgs 04' 223 !########################################################## 224 !### 04. Obtain the present values from the history 225 226 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 227 do ii=1,3 228 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 229 end do 230 231 strten(:)=hist%strten(:,hist%ihist) 232 etotal =hist%etot(hist%ihist) 233 234 !Fill the residual with forces (No preconditioning) 235 !Or the preconditioned forces 236 if (ab_mover%goprecon==0)then 237 call fcart2fred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom) 238 else 239 residual(:,:)=forstr%fred(:,:) 240 end if 241 242 if(zDEBUG)then 243 write (std_out,*) 'residual:' 244 do kk=1,ab_mover%natom 245 write (std_out,*) residual(:,kk) 246 end do 247 write (std_out,*) 'strten:' 248 write (std_out,*) strten(1:3),ch10,strten(4:6) 249 write (std_out,*) 'etotal:' 250 write (std_out,*) etotal 251 end if 252 253 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 254 255 !Save initial values 256 if (itime==1)then 257 acell0(:)=acell(:) 258 rprimd0(:,:)=rprimd(:,:) 259 ucvol0=ucvol 260 end if 261 262 !zDEBUG (UCVOL) 263 if(zDEBUG)then 264 write(std_out,*) 'Volume of cell (ucvol):',ucvol 265 end if 266 267 !Get rid of mean force on whole unit cell, but only if no 268 !generalized constraints are in effect 269 residual_corrected(:,:)=residual(:,:) 270 if(ab_mover%nconeq==0)then 271 do ii=1,3 272 if (ii/=3.or.ab_mover%jellslab==0) then 273 favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom) 274 residual_corrected(ii,:)=residual_corrected(ii,:)-favg 275 end if 276 end do 277 end if 278 279 !write(std_out,*) 'bfgs 05' 280 !########################################################## 281 !### 05. Fill the vectors vin and vout 282 283 !Initialize input vectors : first vin, then vout 284 !The values of vin from the previous iteration 285 !should be the same 286 !if (itime==1)then 287 call xfpack_x2vin(acell, acell0, ab_mover%natom, ndim,& 288 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 289 & ab_mover%symrel, ucvol, ucvol0, vin, xred) 290 !end if 291 292 call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,& 293 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,& 294 & vout) 295 296 !write(std_out,*) 'bfgs 06' 297 !########################################################## 298 !### 06. Initialize or update the hessian matrix 299 300 !Initialise the Hessian matrix using gmet 301 if (itime==1)then 302 303 call hessinit(ab_mover, hessin, gmet, ndim, ucvol) 304 305 ! ! Initialize inverse hessian with identity matrix 306 ! ! in cartesian coordinates, which makes use of metric tensor gmet 307 ! ! in reduced coordinates. 308 ! hessin(:,:)=zero 309 ! do ii=1,ab_mover%natom 310 ! do kk=1,3 311 ! do jj=1,3 312 ! ! Warning : implemented in reduced coordinates 313 ! if (ab_mover%iatfix(kk,ii)==0 .and.& 314 ! & ab_mover%iatfix(jj,ii)==0 )then 315 ! hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj) 316 ! end if 317 ! end do 318 ! end do 319 ! end do 320 ! if(ab_mover%optcell/=0)then 321 ! ! These values might lead to too large changes in some cases 322 ! diag=ab_mover%strprecon*30.0_dp/ucvol 323 ! if(ab_mover%optcell==1) diag=diag/three 324 ! do ii=3*ab_mover%natom+1,ndim 325 ! hessin(ii,ii)=diag 326 ! end do 327 ! end if 328 329 if (ab_mover%restartxf/=0) then 330 331 call xfh_recover_new(ab_xfh,ab_mover,acell,acell0,cycl_main,residual,& 332 & hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,& 333 & vin_prev,vout,vout_prev,xred) 334 335 end if 336 337 end if 338 339 if(itime>1)then 340 ! Update the hessian matrix, by taking into account the 341 ! current pair (x,f) and the previous one. 342 call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,& 343 & vin_prev,vout,vout_prev) 344 345 end if 346 347 !zDEBUG (vin,vout and hessin before prediction) 348 if(zDEBUG)then 349 write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]' 350 write(std_out,*) 'vin:' 351 do ii=1,ndim,3 352 if (ii+2<=ndim)then 353 write(std_out,*) ii,vin(ii:ii+2) 354 else 355 write(std_out,*) ii,vin(ii:ndim) 356 end if 357 end do 358 write(std_out,*) 'vout:' 359 do ii=1,ndim,3 360 if (ii+2<=ndim)then 361 write(std_out,*) ii,vout(ii:ii+2) 362 else 363 write(std_out,*) ii,vout(ii:ndim) 364 end if 365 end do 366 write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim 367 do kk=1,ndim 368 do jj=1,ndim,3 369 if (jj+2<=ndim)then 370 write(std_out,*) jj,hessin(jj:jj+2,kk) 371 else 372 write(std_out,*) jj,hessin(jj:ndim,kk) 373 end if 374 end do 375 end do 376 end if 377 378 !write(std_out,*) 'bfgs 07' 379 !########################################################## 380 !### 07. Compute the next values 381 382 if(ionmov==2 .or. itime==1)then 383 384 ! Previous cartesian coordinates 385 vin_prev(:)=vin(:) 386 387 ! New atomic cartesian coordinates are obtained from vin, hessin 388 ! and vout 389 390 do ii=1,ndim 391 vin(:)=vin(:)-hessin(:,ii)*vout(ii) 392 end do 393 394 !Pulay mixing for vin 395 nitpul=0 396 if (npul>1) then 397 alpha0=1.0_dp 398 nitpul=min(itime, npul) 399 if (itime>npul) then 400 do jj=1,npul-1 401 vinres(jj,:)=vinres(jj+1,:) 402 vin1(jj,:)=vin1(jj+1,:) 403 end do 404 end if 405 vinres(nitpul,:)=vin(:)-vin_prev(:) 406 vin1(nitpul,:)=vin_prev(:) 407 end if 408 409 if (nitpul>1) then 410 ABI_ALLOCATE(alpha,(nitpul,ndim)) 411 alpha=zero 412 do kk=1,ndim 413 ABI_ALLOCATE(amat,(nitpul,nitpul)) 414 ABI_ALLOCATE(amatinv,(nitpul,nitpul)) 415 amat=zero;amatinv=zero 416 do ii=1,nitpul 417 do jj=ii,nitpul 418 amat(ii,jj)=vinres(jj,kk)*vinres(ii,kk) 419 amat(jj,ii)=amat(ii,jj) 420 end do 421 end do 422 amatinv=amat 423 if (abs(vin(kk)-vin_prev(kk))<tol10) then 424 alpha(:,kk)=zero 425 else 426 ABI_ALLOCATE(ipiv,(nitpul)) 427 ABI_ALLOCATE(rwork,(nitpul)) 428 ! amatinv=1.d5*amatinv 429 call dgetrf(nitpul,nitpul,amatinv,nitpul,ipiv,ierr) 430 call dgetri(nitpul,amatinv,nitpul,ipiv,rwork,nitpul,ierr) 431 ! amatinv=1.d5*amatinv 432 ABI_DEALLOCATE(ipiv) 433 ABI_DEALLOCATE(rwork) 434 det=zero 435 do ii=1,nitpul 436 do jj=1,nitpul 437 alpha(ii,kk)=alpha(ii,kk)+amatinv(jj,ii) 438 det=det+amatinv(jj,ii) 439 end do 440 end do 441 alpha(:,kk)=alpha(:,kk)/det 442 end if 443 end do 444 ABI_DEALLOCATE(amat) 445 ABI_DEALLOCATE(amatinv) 446 vin(:)=vin1(nitpul,:)+alpha0*(vin1(nitpul+1,:)-vin1(nitpul,:)) 447 vin=zero 448 do ii=1,nitpul 449 vin(:)=vin(:)+ alpha(ii,:)*(vin1(ii,:)) 450 end do 451 ABI_DEALLOCATE(alpha) 452 end if 453 454 455 ! Previous atomic forces 456 vout_prev(:)=vout(:) 457 458 else 459 if(ionmov==3)then 460 ihist_prev = abihist_findIndex(hist,-1) 461 etotal_prev=hist%etot(ihist_prev) 462 ! Here the BFGS algorithm, modified to take into account the energy 463 call brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev) 464 end if 465 466 ! zDEBUG (vin,vout and hessin after prediction) 467 if(zDEBUG)then 468 write(std_out,*) 'Vectors vin and vout [after prediction]' 469 write(std_out,*) 'vin:' 470 do ii=1,ndim,3 471 if (ii+2<=ndim)then 472 write(std_out,*) ii,vin(ii:ii+2) 473 else 474 write(std_out,*) ii,vin(ii:ndim) 475 end if 476 end do 477 write(std_out,*) 'vout:' 478 do ii=1,ndim,3 479 if (ii+2<=ndim)then 480 write(std_out,*) ii,vout(ii:ii+2) 481 else 482 write(std_out,*) ii,vout(ii:ndim) 483 end if 484 end do 485 end if 486 487 488 ! FIXME: this should be outside the if clause on ionmov! 489 ! Implement fixing of atoms : put back old values for fixed 490 ! components 491 do kk=1,ab_mover%natom 492 do jj=1,3 493 ! Warning : implemented in reduced coordinates 494 if ( ab_mover%iatfix(jj,kk)==1) then 495 vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3) 496 end if 497 end do 498 end do 499 end if 500 501 !write(std_out,*) 'bfgs 08' 502 !########################################################## 503 !### 08. Update the history with the prediction 504 505 !Increase indexes 506 hist%ihist = abihist_findIndex(hist,+1) 507 508 !Transfer vin to xred, acell and rprim 509 call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,& 510 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 511 & ab_mover%symrel, ucvol, ucvol0,& 512 & vin, xred) 513 514 if(ab_mover%optcell/=0)then 515 call mkrdim(acell,rprim,rprimd) 516 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 517 end if 518 519 !Fill the history with the variables 520 !xred, acell, rprimd, vel 521 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 522 ihist_prev = abihist_findIndex(hist,-1) 523 hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev) 524 525 if(zDEBUG)then 526 write (std_out,*) 'residual:' 527 do kk=1,ab_mover%natom 528 write (std_out,*) residual(:,kk) 529 end do 530 write (std_out,*) 'strten:' 531 write (std_out,*) strten(1:3),ch10,strten(4:6) 532 write (std_out,*) 'etotal:' 533 write (std_out,*) etotal 534 end if 535 536 end subroutine pred_bfgs