TABLE OF CONTENTS
ABINIT/m_pred_diisrelax [ Modules ]
NAME
m_pred_diisrelax
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_pred_diisrelax 23 24 use defs_basis 25 use m_abicore 26 use m_abimover 27 use m_abihist 28 use m_linalg_interfaces 29 30 use m_geometry, only : xcart2xred, xred2xcart 31 use m_bfgs, only : hessinit, hessupdt 32 33 implicit none 34 35 private
ABINIT/pred_diisrelax [ Functions ]
NAME
pred_diisrelax
FUNCTION
Ionmov predictor (20) Direct inversion of the iterative subspace IONMOV 20: Given a starting point xred that is a vector of length 3*natom (reduced nuclei coordinates), and unit cell parameters (rprimd) this routine uses the DIIS (direct inversion of the iterative subspace) to minize the gradient (forces) on atoms. The preconditioning used to compute errors from gradients is using an inversed hessian matrix obtained by a BFGS algorithm. This method is known to converge to the nearest point where gradients vanish. This is efficient to refine positions around a saddle point for instance.
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 zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces,acell, rprimd, stresses
SOURCE
75 subroutine pred_diisrelax(ab_mover,hist,itime,ntime,zDEBUG,iexit) 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: itime 80 integer,intent(in) :: ntime 81 integer,intent(in) :: iexit 82 logical,intent(in) :: zDEBUG 83 type(abimover),intent(in) :: ab_mover 84 type(abihist),intent(inout),target :: hist 85 86 !Local variables------------------------------- 87 !scalars 88 integer :: ihist_prev,ndim,nhist,shift,diisSize 89 integer :: ii,jj,kk,info 90 real(dp) :: etotal 91 real(dp) :: suma 92 !arrays 93 real(dp) :: acell(3) 94 real(dp) :: rprimd(3,3) 95 real(dp) :: xred(3,ab_mover%natom) 96 real(dp) :: xcart(3,ab_mover%natom) 97 real(dp) :: strten(6) 98 real(dp) :: ident(3,3) 99 real(dp),allocatable,save :: hessin(:,:) 100 ! DIISRELAX SPECIFIC 101 ! error: Store the supposed error 102 ! steps. it is required to compute the DIIS matrix. 103 ! diisMatrix: Store the matrix used to compute the coefficients. 104 ! diisCoeff: Store the coefficients computed from diisMatrix. 105 ! workMatrix: Lapack work array. 106 ! workArray: Lapack work array. 107 ! ipiv: Lapack work array. 108 integer, allocatable :: ipiv(:) 109 real(dp) :: fcart_tmp(3*ab_mover%natom) 110 real(dp) :: error_tmp(3*ab_mover%natom) 111 real(dp) :: xcart_tmp(3*ab_mover%natom) 112 real(dp), allocatable,save :: error(:, :, :) 113 real(dp), allocatable :: xcart_hist(:,:,:) 114 real(dp), allocatable :: diisMatrix(:, :) 115 real(dp), allocatable :: diisCoeff(:) 116 real(dp), allocatable :: workArray(:) 117 real(dp), allocatable :: workMatrix(:, :) 118 real(dp),pointer :: fcart_hist(:,:,:) 119 120 !*************************************************************************** 121 !Beginning of executable session 122 !*************************************************************************** 123 124 if(iexit/=0)then 125 ABI_SFREE(ipiv) 126 ABI_SFREE(error) 127 ABI_SFREE(diisMatrix) 128 ABI_SFREE(diisCoeff) 129 ABI_SFREE(workArray) 130 ABI_SFREE(workMatrix) 131 ABI_SFREE(hessin) 132 return 133 end if 134 135 !write(std_out,*) 'diisrelax 01' 136 !########################################################## 137 !### 01. Debugging and Verbose 138 139 if(zDEBUG)then 140 write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),& 141 & 'Debugging and Verbose for pred_diisrelax',('-',kk=1,37) 142 write(std_out,*) 'ionmov: ',20 143 write(std_out,*) 'itime: ',itime 144 end if 145 146 !write(std_out,*) 'diisrelax 02' 147 !########################################################## 148 !### 02. Compute the dimension of vectors (ndim=3*natom) 149 150 ndim=3*ab_mover%natom 151 nhist=hist%mxhist 152 153 if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim 154 155 !write(std_out,*) 'diisrelax 03' 156 !########################################################## 157 !### 03. Allocate the arrays 158 159 !Notice that the arrays could be allocated 160 !From a previous dataset with a different ndim 161 if(itime==1)then 162 ABI_SFREE(error) 163 ABI_SFREE(hessin) 164 165 ABI_MALLOC(error,(3, ab_mover%natom, nhist)) 166 ABI_MALLOC(hessin,(ndim,ndim)) 167 168 ident(:, :) = zero 169 ident(1, 1) = -one 170 ident(2, 2) = -one 171 ident(3, 3) = -one 172 call hessinit(ab_mover,hessin,ident,ndim,zero) 173 174 end if 175 176 !write(std_out,*) 'diisrelax 04' 177 !########################################################## 178 !### 04. Compute the shift in the history and the size 179 !### of the diisMatrix 180 181 !When itime > diismemory we need to shift the records 182 !in the history to obtain the right values, the variable 183 !'shift' contains the actual shift to be applied on each 184 !iteration 185 186 shift=max(0,itime-ab_mover%diismemory) 187 188 !Initially the diisMatrix grows with the iteration itime 189 !(itime+1) but when it arrives diismemory, the value of the 190 !matrix will be fixed on (diismemory+1) 191 192 if (ab_mover%diismemory>itime)then 193 diisSize=itime 194 else 195 diisSize=ab_mover%diismemory 196 end if 197 198 199 !write(std_out,*) 'diisrelax 05' 200 !########################################################## 201 !### 05. Obtain the present values from the history 202 203 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 204 205 strten(:)=hist%strten(:,hist%ihist) 206 etotal =hist%etot(hist%ihist) 207 208 if(zDEBUG)then 209 write (std_out,*) 'fcart:' 210 do kk=1,ab_mover%natom 211 write (std_out,*) hist%fcart(:,kk,hist%ihist) 212 end do 213 write (std_out,*) 'strten:' 214 write (std_out,*) strten(1:3),ch10,strten(4:6) 215 write (std_out,*) 'etotal:' 216 write (std_out,*) etotal 217 end if 218 219 !Need also history in cartesian coordinates 220 ABI_MALLOC(xcart_hist,(3,ab_mover%natom,diisSize)) 221 do ii=1,diisSize 222 call xred2xcart(ab_mover%natom,rprimd,xcart_hist(:,:,ii),hist%xred(:,:,ii+shift)) 223 end do 224 fcart_hist => hist%fcart(:,:,1+shift:diisSize+shift) 225 226 !write(std_out,*) 'diisrelax 06' 227 !########################################################## 228 !### 06. Precondition the error using the hessian matrix. 229 230 !Precondition the error using the hessian matrix. 231 !In the quadratic approximation, we have: 232 !e = H^-1.g, where e is the error vectors, H the hessian 233 !and g the gradient. 234 235 if(zDEBUG)then 236 write (std_out,*) 'Stored xcart:' 237 do ii = 1, diisSize, 1 238 write (std_out,*) 'ii,diisSize,shift',ii,diisSize,shift 239 do kk=1,ab_mover%natom 240 write (std_out,*) xcart_hist(:,kk,ii) 241 end do 242 end do 243 write (std_out,*) 'Stored fcart' 244 do ii = 1, diisSize, 1 245 write (std_out,*) 'ii,diisSize,shift',ii,diisSize,shift 246 do kk=1,ab_mover%natom 247 write (std_out,*) fcart_hist(:,kk,ii) 248 end do 249 end do 250 end if 251 252 do ii=1,diisSize,1 253 fcart_tmp(:)=RESHAPE( fcart_hist(:,:,ii), (/ ndim /) ) 254 ! * BLAS ROUTINE LEVEL 2 255 ! * DGEMV performs one of the matrix-vector operations 256 ! * 257 ! * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 258 ! * 259 ! * where alpha and beta are scalars, x and y are vectors and A is an 260 ! * m by n matrix. 261 262 ! Here we are computing: 263 ! error(ndim) := 1*hessin(ndim x ndim)*fcart(ndim) + 0*error(ndim) 264 ! 265 call DGEMV('N',ndim,ndim,one,hessin,& 266 & ndim,fcart_tmp,1,zero,error_tmp,1) 267 error(:,:,ii)=RESHAPE( error_tmp, (/ 3, ab_mover%natom /) ) 268 269 if(zDEBUG)then 270 write (std_out,*) 'Precondition ',ii 271 write (std_out,*) 'fcart_tmp:' 272 do kk=1,3*ab_mover%natom 273 write (std_out,*) fcart_tmp(kk) 274 end do 275 write (std_out,*) 'error:' 276 do kk=1,ab_mover%natom 277 write (std_out,*) error(:,kk,ii) 278 end do 279 write(std_out,*) 'Hessian matrix (hessin):',ndim,'x',ndim 280 do kk=1,ndim 281 do jj=1,ndim,3 282 if (jj+2<=ndim)then 283 write(std_out,'(I3,1p,3e22.14)') jj,hessin(jj:jj+2,kk) 284 else 285 write(std_out,'(I3,1p,3e22.14)') jj,hessin(jj:ndim,kk) 286 end if 287 end do 288 end do 289 end if 290 291 end do 292 293 if(zDEBUG)then 294 write (std_out,*) 'Computed error' 295 do ii = 1, diisSize, 1 296 write (std_out,*) 'ii,diisSize,shift',ii,diisSize,shift 297 do kk=1,ab_mover%natom 298 write (std_out,*) error(:,kk,ii+shift) 299 end do 300 end do 301 end if 302 303 !write(std_out,*) 'diisrelax 07' 304 !########################################################## 305 !### 07. Create the DIIS Matrix 306 307 ABI_MALLOC(diisMatrix,(diisSize + 1, diisSize + 1)) 308 diisMatrix(:,:) = zero 309 if(zDEBUG) write(std_out,*) "DIIS matrix", diisSize+1,'x',diisSize+1 310 do ii = 1, diisSize, 1 311 do jj = ii, diisSize, 1 312 diisMatrix(jj, ii) = ddot(ndim, error(1,1,ii),& 313 & 1, error(1,1,jj),1) 314 diisMatrix(ii, jj) = diisMatrix(jj, ii) 315 end do 316 diisMatrix(ii, diisSize + 1) = -one 317 diisMatrix(diisSize + 1, ii) = -one 318 if(zDEBUG) write(std_out,*) diisMatrix(1:diisSize + 1, ii) 319 end do 320 321 !write(std_out,*) 'diisrelax 08' 322 !########################################################## 323 !### 08. Solve the system using Lapack 324 325 ABI_MALLOC(diisCoeff,(diisSize + 1)) 326 diisCoeff(:) = zero 327 diisCoeff(diisSize + 1) = -one 328 if(zDEBUG) write(std_out,*) "B vector:", diisCoeff(1:diisSize + 1) 329 ABI_MALLOC(workMatrix,(diisSize + 1, diisSize + 1)) 330 ABI_MALLOC(workArray,((diisSize + 1) ** 2)) 331 ABI_MALLOC(ipiv,(diisSize + 1)) 332 !* DCOPY(N,DX,INCX,DY,INCY) 333 !* copies a vector, x, to a vector, y. 334 !* uses unrolled loops for increments equal to one. 335 call DCOPY((diisSize + 1) ** 2, diisMatrix(1:diisSize + 1, 1:diisSize + 1), 1, workMatrix, 1) 336 337 !* DSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, LWORK, INFO ) 338 !* 339 !* Purpose 340 !* ======= 341 !* 342 !* DSYSV computes the solution to a real system of linear equations 343 !* A * X = B, 344 !* where A is an N-by-N symmetric matrix and X and B are N-by-NRHS 345 !* matrices. 346 !* 347 !* The diagonal pivoting method is used to factor A as 348 !* A = U * D * U**T, if UPLO = 'U', or 349 !* A = L * D * L**T, if UPLO = 'L', 350 !* where U (or L) is a product of permutation and unit upper (lower) 351 !* triangular matrices, and D is symmetric and block diagonal with 352 !* 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then 353 !* used to solve the system of equations A * X = B. 354 !* 355 !* Arguments 356 !* ========= 357 !* 358 !* UPLO (input) CHARACTER*1 359 !* = 'U': Upper triangle of A is stored; 360 !* = 'L': Lower triangle of A is stored. 361 !* 362 !* N (input) INTEGER 363 !* The number of linear equations, i.e., the order of the 364 !* matrix A. N >= 0. 365 !* 366 !* NRHS (input) INTEGER 367 !* The number of right hand sides, i.e., the number of columns 368 !* of the matrix B. NRHS >= 0. 369 !* 370 !* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 371 !* On entry, the symmetric matrix A. If UPLO = 'U', the leading 372 !* N-by-N upper triangular part of A contains the upper 373 !* triangular part of the matrix A, and the strictly lower 374 !* triangular part of A is not referenced. If UPLO = 'L', the 375 !* leading N-by-N lower triangular part of A contains the lower 376 !* triangular part of the matrix A, and the strictly upper 377 !* triangular part of A is not referenced. 378 !* 379 !* On exit, if INFO = 0, the block diagonal matrix D and the 380 !* multipliers used to obtain the factor U or L from the 381 !* factorization A = U*D*U**T or A = L*D*L**T as computed by 382 !* DSYTRF. 383 !* 384 !* LDA (input) INTEGER 385 !* The leading dimension of the array A. LDA >= max(1,N). 386 !* 387 !* IPIV (output) INTEGER array, dimension (N) 388 !* Details of the interchanges and the block structure of D, as 389 !* determined by DSYTRF. If IPIV(k) > 0, then rows and columns 390 !* k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1 391 !* diagonal block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, 392 !* then rows and columns k-1 and -IPIV(k) were interchanged and 393 !* D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and 394 !* IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and 395 !* -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 396 !* diagonal block. 397 !* 398 !* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 399 !* On entry, the N-by-NRHS right hand side matrix B. 400 !* On exit, if INFO = 0, the N-by-NRHS solution matrix X. 401 !* 402 !* LDB (input) INTEGER 403 !* The leading dimension of the array B. LDB >= max(1,N). 404 !* 405 !* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 406 !* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 407 !* 408 !* LWORK (input) INTEGER 409 !* The length of WORK. LWORK >= 1, and for best performance 410 !* LWORK >= max(1,N*NB), where NB is the optimal blocksize for 411 !* DSYTRF. 412 !* 413 !* If LWORK = -1, then a workspace query is assumed; the routine 414 !* only calculates the optimal size of the WORK array, returns 415 !* this value as the first entry of the WORK array, and no error 416 !* message related to LWORK is issued by XERBLA. 417 !* 418 !* INFO (output) INTEGER 419 !* = 0: successful exit 420 !* < 0: if INFO = -i, the i-th argument had an illegal value 421 !* > 0: if INFO = i, D(i,i) is exactly zero. The factorization 422 !* has been completed, but the block diagonal matrix D is 423 !* exactly singular, so the solution could not be computed. 424 !* 425 !* ===================================================================== 426 427 if(zDEBUG)then 428 write(std_out,*) "(A*X=B) A Matrix:", diisSize + 1,'x',diisSize + 1 429 do ii=1,diisSize + 1 430 write(std_out,*) workMatrix(1:diisSize + 1,ii) 431 end do 432 433 write(std_out,*) "(A*X=B) B Vector:", diisSize + 1 434 do ii=1,diisSize + 1 435 write(std_out,*) diisCoeff(ii) 436 end do 437 end if 438 439 call DSYSV('L', diisSize + 1, 1, workMatrix, & 440 & diisSize + 1, ipiv, diisCoeff, diisSize + 1, & 441 & workArray, (diisSize + 1) ** 2, info) 442 443 if (info /= 0) then 444 write(std_out,*) "error solving DIIS matrix", info 445 do ii=1,diisSize + 1 446 write(std_out,*) workMatrix(1:diisSize + 1,ii) 447 end do 448 end if 449 450 if(zDEBUG)then 451 write(std_out,*) "(A*X=B) X Vector:",diisSize+1 452 suma=0.0 453 do ii=1,diisSize+1 454 write(std_out,*) ii,diisCoeff(ii) 455 suma=suma+diisCoeff(ii) 456 end do 457 458 suma=suma-diisCoeff(diisSize+1) 459 write(std_out,*) 'Sum of coefficients=',suma 460 end if 461 462 ABI_FREE(ipiv) 463 ABI_FREE(workArray) 464 ABI_FREE(workMatrix) 465 ABI_FREE(diisMatrix) 466 467 !write(std_out,*) 'diisrelax 09' 468 !########################################################## 469 !### 09. Build the new coordinates 470 471 !Build the new coordinates, to do it, we compute a new error e, 472 !using the linear coefficients (temporary store it in error) applied 473 !on previous gradient: e=H^-1(sum_i c_i.g_i) 474 xcart(:, :) = zero 475 error(:, :, diisSize) = zero 476 do ii = 1, diisSize, 1 477 xcart(:, :) = xcart(:, :) +& 478 & xcart_hist(:, :, ii) * diisCoeff(ii) 479 error(:, :, diisSize) = error(:, :, diisSize)+& 480 & fcart_hist(:, :, ii) * diisCoeff(ii) 481 482 if(zDEBUG)then 483 write (std_out,*) 'Building new coordinates (ii):',ii 484 write (std_out,*) 'diisCoeff(ii)',diisCoeff(ii) 485 write (std_out,*) 'xcart_hist(:, :, ii)' 486 do kk=1,ab_mover%natom 487 write (std_out,*) xcart_hist(:,kk,ii) 488 end do 489 write (std_out,*) 'fcart_hist(:, :, ii)' 490 do kk=1,ab_mover%natom 491 write (std_out,*) fcart_hist(:,kk,ii) 492 end do 493 write (std_out,*) 'xcart:' 494 do kk=1,ab_mover%natom 495 write (std_out,*) xcart(:,kk) 496 end do 497 write (std_out,*) 'error:' 498 do kk=1,ab_mover%natom 499 write (std_out,*) error(:,kk,diisSize) 500 end do 501 end if 502 503 end do 504 ABI_FREE(diisCoeff) 505 506 error_tmp(:)=RESHAPE( error(:,:,diisSize), (/ 3*ab_mover%natom /) ) 507 xcart_tmp(:)=RESHAPE( xcart(:,:), (/ 3*ab_mover%natom /) ) 508 509 !* BLAS ROUTINE LEVEL 2 510 !* DGEMV performs one of the matrix-vector operations 511 !* 512 !* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 513 !* 514 !* where alpha and beta are scalars, x and y are vectors and A is an 515 !* m by n matrix. 516 517 !Here we are computing: 518 !xcart_tmp(ndim) := -1*hessin(ndim x ndim)*error_tmp(ndim) + 1*xcart_tmp(ndim) 519 ! 520 call DGEMV('N', ndim, ndim, -one , hessin, & 521 & ndim, error_tmp, 1, one, xcart_tmp, 1) 522 xcart(:,:)=RESHAPE( xcart_tmp(:), (/ 3, ab_mover%natom /) ) 523 524 !write(std_out,*) 'diisrelax 10' 525 !########################################################## 526 !### 10. Update the hessian matrix using a BFGS algorithm. 527 if (itime > 1) then 528 call hessupdt(hessin, ab_mover%iatfix, ab_mover%natom, ndim, & 529 & reshape(xcart_hist(:,:,diisSize) , (/ ndim /)), & 530 & reshape(xcart_hist(:,:,diisSize-1), (/ ndim /)), & 531 & reshape(fcart_hist(:,:,diisSize) , (/ ndim /)), & 532 & reshape(fcart_hist(:,:,diisSize-1), (/ ndim /))) 533 end if 534 535 ABI_FREE(xcart_hist) 536 537 !write(std_out,*) 'diisrelax 11' 538 !########################################################## 539 !### 11. Update the history with the prediction 540 541 !Increase indexes 542 hist%ihist = abihist_findIndex(hist,+1) 543 544 !Compute xred from xcart and rprimd 545 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 546 547 !Fill the history with the variables 548 !xred, acell, rprimd, vel 549 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 550 ihist_prev = abihist_findIndex(hist,-1) 551 hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev) 552 553 if (.false.) write(std_out,*) ntime 554 555 end subroutine pred_diisrelax