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