TABLE OF CONTENTS
ABINIT/m_pred_bfgs [ Modules ]
NAME
m_pred_bfgs
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, JCC, SE, FB) 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_pred_bfgs 22 23 use defs_basis 24 use m_abicore 25 use m_abimover 26 use m_abihist 27 use m_xfpack 28 use m_lbfgs 29 use m_errors 30 31 use m_geometry, only : mkrdim, fcart2gred, metric 32 use m_bfgs, only : hessinit, hessupdt, brdene 33 34 implicit none 35 36 private
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) [[cite:Schlegel1982]]. Might be better than ionmov=2 for few degrees of freedom (less than 3 or 4)
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
SOURCE
89 subroutine pred_bfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit) 90 91 !Arguments ------------------------------------ 92 !scalars 93 type(abimover),intent(in) :: ab_mover 94 type(ab_xfh_type),intent(inout) :: ab_xfh 95 type(abihist),intent(inout) :: hist 96 type(abiforstr),intent(in) :: forstr 97 integer,intent(in) :: itime 98 integer,intent(in) :: ionmov 99 integer,intent(in) :: iexit 100 logical,intent(in) :: zDEBUG 101 102 !Local variables------------------------------- 103 !scalars 104 integer :: ihist_prev,ndim,cycl_main 105 integer, parameter :: npul=0 106 integer :: ierr,ii,jj,kk,nitpul 107 real(dp),save :: ucvol0 108 real(dp) :: ucvol,det 109 real(dp) :: etotal,etotal_prev 110 real(dp) :: favg,alpha0 111 112 !arrays 113 integer,allocatable :: ipiv(:) 114 real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:) 115 real(dp),allocatable,save :: vout(:),vout_prev(:) 116 real(dp),allocatable,save ::vinres(:,:),vin1(:,:) 117 real(dp),allocatable :: amat(:,:),amatinv(:,:),alpha(:,:) 118 real(dp),allocatable :: rwork(:) 119 real(dp),save :: acell0(3) ! Initial acell 120 real(dp),save :: rprimd0(3,3) ! Initial rprimd 121 real(dp) :: acell(3) 122 real(dp) :: rprimd(3,3),rprim(3,3) 123 real(dp) :: gprimd(3,3) 124 real(dp) :: gmet(3,3) 125 real(dp) :: rmet(3,3) 126 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom) 127 real(dp) :: xred(3,ab_mover%natom),strten(6) 128 129 !*************************************************************************** 130 !Beginning of executable session 131 !*************************************************************************** 132 133 if(iexit/=0)then 134 ABI_SFREE(vin) 135 ABI_SFREE(vout) 136 ABI_SFREE(vin_prev) 137 ABI_SFREE(vout_prev) 138 ABI_SFREE(vinres) 139 ABI_SFREE(vin1) 140 ABI_SFREE(hessin) 141 return 142 end if 143 144 !write(std_out,*) 'bfgs 01' 145 !########################################################## 146 !### 01. Debugging and Verbose 147 148 if(zDEBUG)then 149 write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),'Debugging and Verbose for pred_bfgs',('-',kk=1,42) 150 write(std_out,*) 'ionmov: ',ionmov 151 write(std_out,*) 'itime: ',itime 152 end if 153 154 !write(std_out,*) 'bfgs 02' 155 !########################################################## 156 !### 02. Compute the dimension of vectors (ndim) 157 158 ndim=3*ab_mover%natom 159 if(ab_mover%optcell==1) ndim=ndim+1 160 if(ab_mover%optcell==2 .or.& 161 & ab_mover%optcell==3) ndim=ndim+6 162 if(ab_mover%optcell>=4) ndim=ndim+3 163 164 if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim 165 166 !write(std_out,*) 'bfgs 03' 167 !########################################################## 168 !### 03. Allocate the vectors vin, vout and hessian matrix 169 170 !Notice that vin, vout, etc could be allocated 171 !From a previous dataset with a different ndim 172 if(itime==1)then 173 ABI_SFREE(vin) 174 ABI_SFREE(vout) 175 ABI_SFREE(vin_prev) 176 ABI_SFREE(vout_prev) 177 ABI_SFREE(vinres) 178 ABI_SFREE(vin1) 179 ABI_SFREE(hessin) 180 if(npul>1) then 181 ABI_MALLOC(vinres,(npul+1,ndim)) 182 ABI_MALLOC(vin1,(npul+1,ndim)) 183 end if 184 ABI_MALLOC(vin,(ndim)) 185 ABI_MALLOC(vout,(ndim)) 186 ABI_MALLOC(vin_prev,(ndim)) 187 ABI_MALLOC(vout_prev,(ndim)) 188 ABI_MALLOC(hessin,(ndim,ndim)) 189 end if 190 191 !write(std_out,*) 'bfgs 04' 192 !########################################################## 193 !### 04. Obtain the present values from the history 194 195 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 196 do ii=1,3 197 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 198 end do 199 200 strten(:)=hist%strten(:,hist%ihist) 201 etotal =hist%etot(hist%ihist) 202 203 !Fill the residual with forces (No preconditioning) 204 !Or the preconditioned forces 205 if (ab_mover%goprecon==0)then 206 call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom) 207 else 208 residual(:,:)=forstr%gred(:,:) 209 end if 210 211 if(zDEBUG)then 212 write (std_out,*) 'residual:' 213 do kk=1,ab_mover%natom 214 write (std_out,*) residual(:,kk) 215 end do 216 write (std_out,*) 'strten:' 217 write (std_out,*) strten(1:3),ch10,strten(4:6) 218 write (std_out,*) 'etotal:' 219 write (std_out,*) etotal 220 end if 221 222 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 223 224 !Save initial values 225 if (itime==1)then 226 acell0(:)=acell(:) 227 rprimd0(:,:)=rprimd(:,:) 228 ucvol0=ucvol 229 end if 230 231 !zDEBUG (UCVOL) 232 if(zDEBUG)then 233 write(std_out,*) 'Volume of cell (ucvol):',ucvol 234 end if 235 236 !Get rid of mean force on whole unit cell, but only if no 237 !generalized constraints are in effect 238 residual_corrected(:,:)=residual(:,:) 239 if(ab_mover%nconeq==0)then 240 do ii=1,3 241 if (ii/=3.or.ab_mover%jellslab==0) then 242 favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom) 243 residual_corrected(ii,:)=residual_corrected(ii,:)-favg 244 end if 245 end do 246 end if 247 248 !write(std_out,*) 'bfgs 05' 249 !########################################################## 250 !### 05. Fill the vectors vin and vout 251 252 !Initialize input vectors : first vin, then vout 253 !The values of vin from the previous iteration 254 !should be the same 255 !if (itime==1)then 256 call xfpack_x2vin(acell, ab_mover%natom, ndim,& 257 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 258 & ab_mover%symrel, ucvol, ucvol0, vin, xred) 259 !end if 260 261 call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,& 262 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,& 263 & vout) 264 265 !write(std_out,*) 'bfgs 06' 266 !########################################################## 267 !### 06. Initialize or update the hessian matrix 268 269 !Initialise the Hessian matrix using gmet 270 if (itime==1)then 271 272 call hessinit(ab_mover, hessin, gmet, ndim, ucvol) 273 274 ! ! Initialize inverse hessian with identity matrix 275 ! ! in cartesian coordinates, which makes use of metric tensor gmet 276 ! ! in reduced coordinates. 277 ! hessin(:,:)=zero 278 ! do ii=1,ab_mover%natom 279 ! do kk=1,3 280 ! do jj=1,3 281 ! ! Warning : implemented in reduced coordinates 282 ! if (ab_mover%iatfix(kk,ii)==0 .and.& 283 ! & ab_mover%iatfix(jj,ii)==0 )then 284 ! hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj) 285 ! end if 286 ! end do 287 ! end do 288 ! end do 289 ! if(ab_mover%optcell/=0)then 290 ! ! These values might lead to too large changes in some cases 291 ! diag=ab_mover%strprecon*30.0_dp/ucvol 292 ! if(ab_mover%optcell==1) diag=diag/three 293 ! do ii=3*ab_mover%natom+1,ndim 294 ! hessin(ii,ii)=diag 295 ! end do 296 ! end if 297 298 if (ab_mover%restartxf/=0) then 299 300 call xfh_recover_new(ab_xfh,ab_mover,acell,cycl_main,residual,& 301 & hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,& 302 & vin_prev,vout,vout_prev,xred) 303 304 end if 305 306 end if 307 308 if(itime>1)then 309 ! Update the hessian matrix, by taking into account the 310 ! current pair (x,f) and the previous one. 311 call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,& 312 & vin_prev,vout,vout_prev) 313 314 end if 315 316 !zDEBUG (vin,vout and hessin before prediction) 317 if(zDEBUG)then 318 write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]' 319 write(std_out,*) 'vin:' 320 do ii=1,ndim,3 321 if (ii+2<=ndim)then 322 write(std_out,*) ii,vin(ii:ii+2) 323 else 324 write(std_out,*) ii,vin(ii:ndim) 325 end if 326 end do 327 write(std_out,*) 'vout:' 328 do ii=1,ndim,3 329 if (ii+2<=ndim)then 330 write(std_out,*) ii,vout(ii:ii+2) 331 else 332 write(std_out,*) ii,vout(ii:ndim) 333 end if 334 end do 335 write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim 336 do kk=1,ndim 337 do jj=1,ndim,3 338 if (jj+2<=ndim)then 339 write(std_out,*) jj,hessin(jj:jj+2,kk) 340 else 341 write(std_out,*) jj,hessin(jj:ndim,kk) 342 end if 343 end do 344 end do 345 end if 346 347 !write(std_out,*) 'bfgs 07' 348 !########################################################## 349 !### 07. Compute the next values 350 351 if(ionmov==2 .or. itime==1)then 352 353 ! Previous cartesian coordinates 354 vin_prev(:)=vin(:) 355 356 ! New atomic cartesian coordinates are obtained from vin, hessin 357 ! and vout 358 359 do ii=1,ndim 360 vin(:)=vin(:)-hessin(:,ii)*vout(ii) 361 end do 362 363 !Pulay mixing for vin 364 nitpul=0 365 if (npul>1) then 366 alpha0=1.0_dp 367 nitpul=min(itime, npul) 368 if (itime>npul) then 369 do jj=1,npul-1 370 vinres(jj,:)=vinres(jj+1,:) 371 vin1(jj,:)=vin1(jj+1,:) 372 end do 373 end if 374 vinres(nitpul,:)=vin(:)-vin_prev(:) 375 vin1(nitpul,:)=vin_prev(:) 376 end if 377 378 if (nitpul>1) then 379 ABI_MALLOC(alpha,(nitpul,ndim)) 380 alpha=zero 381 do kk=1,ndim 382 ABI_MALLOC(amat,(nitpul,nitpul)) 383 ABI_MALLOC(amatinv,(nitpul,nitpul)) 384 amat=zero;amatinv=zero 385 do ii=1,nitpul 386 do jj=ii,nitpul 387 amat(ii,jj)=vinres(jj,kk)*vinres(ii,kk) 388 amat(jj,ii)=amat(ii,jj) 389 end do 390 end do 391 amatinv=amat 392 if (abs(vin(kk)-vin_prev(kk))<tol10) then 393 alpha(:,kk)=zero 394 else 395 ABI_MALLOC(ipiv,(nitpul)) 396 ABI_MALLOC(rwork,(nitpul)) 397 ! amatinv=1.d5*amatinv 398 call dgetrf(nitpul,nitpul,amatinv,nitpul,ipiv,ierr) 399 call dgetri(nitpul,amatinv,nitpul,ipiv,rwork,nitpul,ierr) 400 ! amatinv=1.d5*amatinv 401 ABI_FREE(ipiv) 402 ABI_FREE(rwork) 403 det=zero 404 do ii=1,nitpul 405 do jj=1,nitpul 406 alpha(ii,kk)=alpha(ii,kk)+amatinv(jj,ii) 407 det=det+amatinv(jj,ii) 408 end do 409 end do 410 alpha(:,kk)=alpha(:,kk)/det 411 end if 412 end do 413 ABI_FREE(amat) 414 ABI_FREE(amatinv) 415 vin(:)=vin1(nitpul,:)+alpha0*(vin1(nitpul+1,:)-vin1(nitpul,:)) 416 vin=zero 417 do ii=1,nitpul 418 vin(:)=vin(:)+ alpha(ii,:)*(vin1(ii,:)) 419 end do 420 ABI_FREE(alpha) 421 end if 422 423 424 ! Previous atomic forces 425 vout_prev(:)=vout(:) 426 427 else 428 if(ionmov==3)then 429 ihist_prev = abihist_findIndex(hist,-1) 430 etotal_prev=hist%etot(ihist_prev) 431 ! Here the BFGS algorithm, modified to take into account the energy 432 call brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev) 433 end if 434 435 ! zDEBUG (vin,vout and hessin after prediction) 436 if(zDEBUG)then 437 write(std_out,*) 'Vectors vin and vout [after prediction]' 438 write(std_out,*) 'vin:' 439 do ii=1,ndim,3 440 if (ii+2<=ndim)then 441 write(std_out,*) ii,vin(ii:ii+2) 442 else 443 write(std_out,*) ii,vin(ii:ndim) 444 end if 445 end do 446 write(std_out,*) 'vout:' 447 do ii=1,ndim,3 448 if (ii+2<=ndim)then 449 write(std_out,*) ii,vout(ii:ii+2) 450 else 451 write(std_out,*) ii,vout(ii:ndim) 452 end if 453 end do 454 end if 455 456 457 ! FIXME: this should be outside the if clause on ionmov! 458 ! Implement fixing of atoms : put back old values for fixed 459 ! components 460 do kk=1,ab_mover%natom 461 do jj=1,3 462 ! Warning : implemented in reduced coordinates 463 if ( ab_mover%iatfix(jj,kk)==1) then 464 vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3) 465 end if 466 end do 467 end do 468 end if 469 470 !write(std_out,*) 'bfgs 08' 471 !########################################################## 472 !### 08. Update the history with the prediction 473 474 !Increase indexes 475 hist%ihist = abihist_findIndex(hist,+1) 476 477 !Transfer vin to xred, acell and rprim 478 call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,& 479 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 480 & ab_mover%symrel, ucvol, ucvol0,& 481 & vin, xred) 482 483 if(ab_mover%optcell/=0)then 484 call mkrdim(acell,rprim,rprimd) 485 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 486 end if 487 488 !Fill the history with the variables 489 !xred, acell, rprimd, vel 490 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 491 ihist_prev = abihist_findIndex(hist,-1) 492 hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev) 493 494 if(zDEBUG)then 495 write (std_out,*) 'residual:' 496 do kk=1,ab_mover%natom 497 write (std_out,*) residual(:,kk) 498 end do 499 write (std_out,*) 'strten:' 500 write (std_out,*) strten(1:3),ch10,strten(4:6) 501 write (std_out,*) 'etotal:' 502 write (std_out,*) etotal 503 end if 504 505 end subroutine pred_bfgs
ABINIT/pred_lbfgs [ Functions ]
NAME
pred_lbfgs
FUNCTION
Ionmov predictors (22) Limited-memory Broyden-Fletcher-Goldfarb-Shanno IONMOV 22: Given a starting point xred that is a vector of length 3*natom (reduced nuclei coordinates), and unit cell parameters (acell and rprim) the L-Broyden-Fletcher-Goldfarb-Shanno minimization is performed on the total energy function, using its gradient (atomic forces and stress : gred or fcart and stress) 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 MinPack on netlib.org
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 : (22) Specific kind of BFGS zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
543 subroutine pred_lbfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit) 544 545 !Arguments ------------------------------------ 546 !scalars 547 type(abimover),intent(in) :: ab_mover 548 type(ab_xfh_type),intent(inout) :: ab_xfh 549 type(abihist),intent(inout) :: hist 550 type(abiforstr),intent(in) :: forstr 551 integer,intent(in) :: itime 552 integer,intent(in) :: ionmov 553 integer,intent(in) :: iexit 554 logical,intent(in) :: zDEBUG 555 character(len=500) :: ionmov22_errmsg 556 557 !Local variables------------------------------- 558 !scalars 559 integer :: info,ihist_prev 560 integer :: ndim,cycl_main 561 integer, parameter :: npul=0 562 integer :: ii,jj,kk 563 real(dp),save :: ucvol0 564 real(dp) :: ucvol 565 real(dp) :: etotal 566 real(dp) :: favg 567 568 !arrays 569 real(dp),allocatable :: diag(:) 570 real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:) 571 real(dp),allocatable,save :: vout(:),vout_prev(:) 572 real(dp),allocatable,save ::vinres(:,:),vin1(:,:) 573 real(dp),save :: acell0(3) ! Initial acell 574 real(dp),save :: rprimd0(3,3) ! Initial rprimd 575 real(dp) :: acell(3) 576 real(dp) :: rprimd(3,3),rprim(3,3) 577 real(dp) :: gprimd(3,3) 578 real(dp) :: gmet(3,3) 579 real(dp) :: rmet(3,3) 580 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom) 581 real(dp) :: xred(3,ab_mover%natom) 582 real(dp) :: strten(6) 583 584 !*************************************************************************** 585 !Beginning of executable session 586 !*************************************************************************** 587 588 if(iexit/=0)then 589 call lbfgs_destroy() 590 ABI_SFREE(vin) 591 ABI_SFREE(vout) 592 ABI_SFREE(vin_prev) 593 ABI_SFREE(vout_prev) 594 ABI_SFREE(vinres) 595 ABI_SFREE(vin1) 596 ABI_SFREE(hessin) 597 return 598 end if 599 600 !write(std_out,*) 'bfgs 01' 601 !########################################################## 602 !### 01. Debugging and Verbose 603 604 if(zDEBUG)then 605 write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),'Debugging and Verbose for pred_bfgs',('-',kk=1,42) 606 write(std_out,*) 'ionmov: ',ionmov 607 write(std_out,*) 'itime: ',itime 608 end if 609 610 !write(std_out,*) 'bfgs 02' 611 !########################################################## 612 !### 02. Compute the dimension of vectors (ndim) 613 614 ndim=3*ab_mover%natom 615 if(ab_mover%optcell==1) ndim=ndim+1 616 if(ab_mover%optcell==2 .or.& 617 & ab_mover%optcell==3) ndim=ndim+6 618 if(ab_mover%optcell>=4) ndim=ndim+3 619 620 if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim 621 622 !write(std_out,*) 'bfgs 03' 623 !########################################################## 624 !### 03. Allocate the vectors vin, vout and hessian matrix 625 626 !Notice that vin, vout, etc could be allocated 627 !From a previous dataset with a different ndim 628 if(itime==1)then 629 ABI_SFREE(vin) 630 ABI_SFREE(vout) 631 ABI_SFREE(vin_prev) 632 ABI_SFREE(vout_prev) 633 ABI_SFREE(vinres) 634 ABI_SFREE(vin1) 635 ABI_SFREE(hessin) 636 if(npul>1) then 637 ABI_MALLOC(vinres,(npul+1,ndim)) 638 ABI_MALLOC(vin1,(npul+1,ndim)) 639 end if 640 ABI_MALLOC(vin,(ndim)) 641 ABI_MALLOC(vout,(ndim)) 642 ABI_MALLOC(vin_prev,(ndim)) 643 ABI_MALLOC(vout_prev,(ndim)) 644 ABI_MALLOC(hessin,(ndim,ndim)) 645 end if 646 647 !write(std_out,*) 'bfgs 04' 648 !########################################################## 649 !### 04. Obtain the present values from the history 650 651 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 652 do ii=1,3 653 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 654 end do 655 656 strten(:)=hist%strten(:,hist%ihist) 657 etotal =hist%etot(hist%ihist) 658 659 !Fill the residual with forces (No preconditioning) 660 !Or the preconditioned forces 661 if (ab_mover%goprecon==0)then 662 call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom) 663 else 664 residual(:,:)= forstr%gred(:,:) 665 end if 666 667 if(zDEBUG)then 668 write (std_out,*) 'residual:' 669 do kk=1,ab_mover%natom 670 write (std_out,*) residual(:,kk) 671 end do 672 write (std_out,*) 'strten:' 673 write (std_out,*) strten(1:3),ch10,strten(4:6) 674 write (std_out,*) 'etotal:' 675 write (std_out,*) etotal 676 end if 677 678 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 679 680 !Save initial values 681 if (itime==1)then 682 acell0(:)=acell(:) 683 rprimd0(:,:)=rprimd(:,:) 684 ucvol0=ucvol 685 end if 686 687 !zDEBUG (UCVOL) 688 if(zDEBUG)then 689 write(std_out,*) 'Volume of cell (ucvol):',ucvol 690 end if 691 692 !Get rid of mean force on whole unit cell, but only if no 693 !generalized constraints are in effect 694 residual_corrected(:,:)=residual(:,:) 695 if(ab_mover%nconeq==0)then 696 do ii=1,3 697 if (ii/=3.or.ab_mover%jellslab==0) then 698 favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom) 699 residual_corrected(ii,:)=residual_corrected(ii,:)-favg 700 end if 701 end do 702 end if 703 704 !write(std_out,*) 'bfgs 05' 705 !########################################################## 706 !### 05. Fill the vectors vin and vout 707 708 !Initialize input vectors : first vin, then vout 709 !The values of vin from the previous iteration 710 !should be the same 711 !if (itime==1)then 712 call xfpack_x2vin(acell, ab_mover%natom, ndim,& 713 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 714 & ab_mover%symrel, ucvol, ucvol0, vin, xred) 715 !end if 716 717 call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,& 718 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,& 719 & vout) 720 721 !write(std_out,*) 'bfgs 06' 722 !########################################################## 723 !### 06. Initialize or update the hessian matrix 724 725 !Initialise the Hessian matrix using gmet 726 if (itime==1)then 727 728 ABI_MALLOC(diag,(ndim)) 729 do ii=1,3*ab_mover%natom 730 !diag(ii) = 1.00_dp / rprimd(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1)**2 731 diag(ii) = gmet(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1) 732 end do 733 if(ab_mover%optcell/=0)then 734 ! These values might lead to too large changes in some cases ... 735 do ii=3*ab_mover%natom+1,ndim 736 diag(ii) = ab_mover%strprecon*30.0_dp/ucvol 737 if(ab_mover%optcell==1) diag(ii) = diag(ii) / three 738 end do 739 end if 740 741 !call lbfgs_destroy() 742 call lbfgs_init(ndim,5,diag) 743 ABI_FREE(diag) 744 745 if (ab_mover%restartxf/=0) then 746 747 call xfh_recover_new(ab_xfh,ab_mover,acell,cycl_main,residual,& 748 hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,& 749 vin_prev,vout,vout_prev,xred) 750 751 end if 752 753 end if 754 755 !zDEBUG (vin,vout and hessin before prediction) 756 if(zDEBUG)then 757 write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]' 758 write(std_out,*) 'vin:' 759 do ii=1,ndim,3 760 if (ii+2<=ndim)then 761 write(std_out,*) ii,vin(ii:ii+2) 762 else 763 write(std_out,*) ii,vin(ii:ndim) 764 end if 765 end do 766 write(std_out,*) 'vout:' 767 do ii=1,ndim,3 768 if (ii+2<=ndim)then 769 write(std_out,*) ii,vout(ii:ii+2) 770 else 771 write(std_out,*) ii,vout(ii:ndim) 772 end if 773 end do 774 write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim 775 do kk=1,ndim 776 do jj=1,ndim,3 777 if (jj+2<=ndim)then 778 write(std_out,*) jj,hessin(jj:jj+2,kk) 779 else 780 write(std_out,*) jj,hessin(jj:ndim,kk) 781 end if 782 end do 783 end do 784 end if 785 786 !write(std_out,*) 'bfgs 07' 787 !########################################################## 788 !### 07. Compute the next values 789 790 vin_prev(:) = vin 791 vout_prev(:) = vout 792 info = lbfgs_execute(vin,etotal,vout) 793 794 if (info /= -1) then 795 write (ionmov22_errmsg, '(a,i0,3a)') & 796 'Lbfgs routine failed. Returned value: ', info,ch10, & 797 'Restart your calculation from last step or try a different ionmov' 798 ABI_ERROR_CLASS(ionmov22_errmsg, "Ionmov22Error") 799 end if 800 801 !zDEBUG (vin,vout after prediction) 802 if(zDEBUG)then 803 write(std_out,*) 'Vectors vin and vout [after prediction]' 804 write(std_out,*) 'vin_prev:' 805 do ii=1,ndim,3 806 if (ii+2<=ndim)then 807 write(std_out,*) ii,vin_prev(ii:ii+2) 808 else 809 write(std_out,*) ii,vin_prev(ii:ndim) 810 end if 811 end do 812 write(std_out,*) 'vin:' 813 do ii=1,ndim,3 814 if (ii+2<=ndim)then 815 write(std_out,*) ii,vin(ii:ii+2) 816 else 817 write(std_out,*) ii,vin(ii:ndim) 818 end if 819 end do 820 write(std_out,*) 'vout:' 821 do ii=1,ndim,3 822 if (ii+2<=ndim)then 823 write(std_out,*) ii,vout(ii:ii+2) 824 else 825 write(std_out,*) ii,vout(ii:ndim) 826 end if 827 end do 828 end if 829 830 831 !Implement fixing of atoms : put back old values for fixed 832 !components 833 do kk=1,ab_mover%natom 834 do jj=1,3 835 ! Warning : implemented in reduced coordinates 836 if ( ab_mover%iatfix(jj,kk)==1) then 837 vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3) 838 end if 839 end do 840 end do 841 842 843 !write(std_out,*) 'bfgs 08' 844 !########################################################## 845 !### 08. Update the history with the prediction 846 847 !Increase indexes 848 hist%ihist = abihist_findIndex(hist,+1) 849 850 !Transfer vin to xred, acell and rprim 851 call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,& 852 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 853 & ab_mover%symrel, ucvol, ucvol0,& 854 & vin, xred) 855 856 if(ab_mover%optcell/=0)then 857 call mkrdim(acell,rprim,rprimd) 858 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 859 end if 860 861 !Fill the history with the variables 862 !xcart, xred, acell, rprimd 863 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 864 ihist_prev = abihist_findIndex(hist,-1) 865 hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev) 866 867 if(zDEBUG)then 868 write (std_out,*) 'residual:' 869 do kk=1,ab_mover%natom 870 write (std_out,*) residual(:,kk) 871 end do 872 write (std_out,*) 'strten:' 873 write (std_out,*) strten(1:3),ch10,strten(4:6) 874 write (std_out,*) 'etotal:' 875 write (std_out,*) etotal 876 end if 877 878 end subroutine pred_lbfgs