Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_diisrelax
2 : : !! NAME
3 : : !! m_pred_diisrelax
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JCC, SE)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
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
36 : : !!***
37 : :
38 : : public :: pred_diisrelax
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* ABINIT/pred_diisrelax
45 : : !! NAME
46 : : !! pred_diisrelax
47 : : !!
48 : : !! FUNCTION
49 : : !! Ionmov predictor (20) Direct inversion of the iterative subspace
50 : : !!
51 : : !! IONMOV 20:
52 : : !! Given a starting point xred that is a vector of length 3*natom
53 : : !! (reduced nuclei coordinates), and unit cell parameters (rprimd)
54 : : !! this routine uses the DIIS (direct inversion of the iterative
55 : : !! subspace) to minize the gradient (forces) on atoms. The preconditioning
56 : : !! used to compute errors from gradients is using an inversed hessian
57 : : !! matrix obtained by a BFGS algorithm.
58 : : !! This method is known to converge to the nearest point where gradients
59 : : !! vanish. This is efficient to refine positions around a saddle point
60 : : !! for instance.
61 : : !!
62 : : !! INPUTS
63 : : !! ab_mover <type(abimover)> : Datatype with all the information needed by the preditor
64 : : !! itime : Index of the present iteration
65 : : !! ntime : Maximal number of iterations
66 : : !! zDEBUG : if true print some debugging information
67 : : !!
68 : : !! OUTPUT
69 : : !!
70 : : !! SIDE EFFECTS
71 : : !! hist <type(abihist)> : History of positions,forces,acell, rprimd, stresses
72 : : !!
73 : : !! SOURCE
74 : :
75 : 0 : 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 : 0 : real(dp) :: xred(3,ab_mover%natom)
96 : 0 : 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 : 0 : integer, allocatable :: ipiv(:)
109 : 0 : real(dp) :: fcart_tmp(3*ab_mover%natom)
110 : 0 : real(dp) :: error_tmp(3*ab_mover%natom)
111 [ # # ]: 0 : real(dp) :: xcart_tmp(3*ab_mover%natom)
112 : : real(dp), allocatable,save :: error(:, :, :)
113 : 0 : real(dp), allocatable :: xcart_hist(:,:,:)
114 : 0 : real(dp), allocatable :: diisMatrix(:, :)
115 : 0 : real(dp), allocatable :: diisCoeff(:)
116 : 0 : real(dp), allocatable :: workArray(:)
117 : 0 : real(dp), allocatable :: workMatrix(:, :)
118 : 0 : real(dp),pointer :: fcart_hist(:,:,:)
119 : :
120 : : !***************************************************************************
121 : : !Beginning of executable session
122 : : !***************************************************************************
123 : :
124 : 0 : if(iexit/=0)then
125 : : ABI_SFREE(ipiv)
126 [ # # ]: 0 : ABI_SFREE(error)
127 : : ABI_SFREE(diisMatrix)
128 : : ABI_SFREE(diisCoeff)
129 : : ABI_SFREE(workArray)
130 : : ABI_SFREE(workMatrix)
131 [ # # ]: 0 : ABI_SFREE(hessin)
132 : 0 : return
133 : : end if
134 : :
135 : : !write(std_out,*) 'diisrelax 01'
136 : : !##########################################################
137 : : !### 01. Debugging and Verbose
138 : :
139 [ # # ]: 0 : if(zDEBUG)then
140 [ # # ]: 0 : write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),&
141 [ # # # # : 0 : & 'Debugging and Verbose for pred_diisrelax',('-',kk=1,37)
# # ]
142 : 0 : write(std_out,*) 'ionmov: ',20
143 : 0 : 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 : 0 : ndim=3*ab_mover%natom
151 : 0 : nhist=hist%mxhist
152 : :
153 [ # # ]: 0 : 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 [ # # ]: 0 : if(itime==1)then
162 [ # # ]: 0 : ABI_SFREE(error)
163 [ # # ]: 0 : ABI_SFREE(hessin)
164 : :
165 [ # # # # : 0 : ABI_MALLOC(error,(3, ab_mover%natom, nhist))
# # # # #
# # # # #
# # # # ]
166 [ # # # # : 0 : ABI_MALLOC(hessin,(ndim,ndim))
# # # # #
# # # # #
# # # # ]
167 : :
168 : 0 : ident(:, :) = zero
169 : 0 : ident(1, 1) = -one
170 : 0 : ident(2, 2) = -one
171 : 0 : ident(3, 3) = -one
172 : 0 : 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 : 0 : 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 [ # # ]: 0 : if (ab_mover%diismemory>itime)then
193 : 0 : diisSize=itime
194 : : else
195 : 0 : 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 : 0 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
204 : :
205 [ # # ]: 0 : strten(:)=hist%strten(:,hist%ihist)
206 : 0 : etotal =hist%etot(hist%ihist)
207 : :
208 [ # # ]: 0 : if(zDEBUG)then
209 : 0 : write (std_out,*) 'fcart:'
210 [ # # ]: 0 : do kk=1,ab_mover%natom
211 : 0 : write (std_out,*) hist%fcart(:,kk,hist%ihist)
212 : : end do
213 : 0 : write (std_out,*) 'strten:'
214 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
215 : 0 : write (std_out,*) 'etotal:'
216 : 0 : write (std_out,*) etotal
217 : : end if
218 : :
219 : : !Need also history in cartesian coordinates
220 [ # # # # : 0 : ABI_MALLOC(xcart_hist,(3,ab_mover%natom,diisSize))
# # # # #
# # # # #
# # ]
221 [ # # ]: 0 : do ii=1,diisSize
222 : 0 : call xred2xcart(ab_mover%natom,rprimd,xcart_hist(:,:,ii),hist%xred(:,:,ii+shift))
223 : : end do
224 : 0 : 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 [ # # ]: 0 : if(zDEBUG)then
236 : 0 : write (std_out,*) 'Stored xcart:'
237 [ # # ]: 0 : do ii = 1, diisSize, 1
238 : 0 : write (std_out,*) 'ii,diisSize,shift',ii,diisSize,shift
239 [ # # ]: 0 : do kk=1,ab_mover%natom
240 : 0 : write (std_out,*) xcart_hist(:,kk,ii)
241 : : end do
242 : : end do
243 : 0 : write (std_out,*) 'Stored fcart'
244 [ # # ]: 0 : do ii = 1, diisSize, 1
245 : 0 : write (std_out,*) 'ii,diisSize,shift',ii,diisSize,shift
246 [ # # ]: 0 : do kk=1,ab_mover%natom
247 : 0 : write (std_out,*) fcart_hist(:,kk,ii)
248 : : end do
249 : : end do
250 : : end if
251 : :
252 [ # # ]: 0 : do ii=1,diisSize,1
253 [ # # ]: 0 : 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 : 0 : & ndim,fcart_tmp,1,zero,error_tmp,1)
267 [ # # ]: 0 : error(:,:,ii)=RESHAPE( error_tmp, (/ 3, ab_mover%natom /) )
268 : :
269 [ # # ]: 0 : if(zDEBUG)then
270 : 0 : write (std_out,*) 'Precondition ',ii
271 : 0 : write (std_out,*) 'fcart_tmp:'
272 [ # # ]: 0 : do kk=1,3*ab_mover%natom
273 : 0 : write (std_out,*) fcart_tmp(kk)
274 : : end do
275 : 0 : write (std_out,*) 'error:'
276 [ # # ]: 0 : do kk=1,ab_mover%natom
277 : 0 : write (std_out,*) error(:,kk,ii)
278 : : end do
279 : 0 : write(std_out,*) 'Hessian matrix (hessin):',ndim,'x',ndim
280 [ # # ]: 0 : do kk=1,ndim
281 [ # # # # ]: 0 : do jj=1,ndim,3
282 [ # # ]: 0 : if (jj+2<=ndim)then
283 : 0 : write(std_out,'(I3,1p,3e22.14)') jj,hessin(jj:jj+2,kk)
284 : : else
285 : 0 : 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 [ # # ]: 0 : if(zDEBUG)then
294 : 0 : write (std_out,*) 'Computed error'
295 [ # # ]: 0 : do ii = 1, diisSize, 1
296 : 0 : write (std_out,*) 'ii,diisSize,shift',ii,diisSize,shift
297 [ # # ]: 0 : do kk=1,ab_mover%natom
298 : 0 : 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 [ # # # # : 0 : ABI_MALLOC(diisMatrix,(diisSize + 1, diisSize + 1))
# # # # #
# # # # #
# # # # ]
308 [ # # # # ]: 0 : diisMatrix(:,:) = zero
309 [ # # ]: 0 : if(zDEBUG) write(std_out,*) "DIIS matrix", diisSize+1,'x',diisSize+1
310 [ # # ]: 0 : do ii = 1, diisSize, 1
311 [ # # ]: 0 : do jj = ii, diisSize, 1
312 : : diisMatrix(jj, ii) = ddot(ndim, error(1,1,ii),&
313 : 0 : & 1, error(1,1,jj),1)
314 : 0 : diisMatrix(ii, jj) = diisMatrix(jj, ii)
315 : : end do
316 : 0 : diisMatrix(ii, diisSize + 1) = -one
317 : 0 : diisMatrix(diisSize + 1, ii) = -one
318 [ # # ]: 0 : 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 [ # # # # : 0 : ABI_MALLOC(diisCoeff,(diisSize + 1))
# # # # #
# # # ]
326 [ # # ]: 0 : diisCoeff(:) = zero
327 : 0 : diisCoeff(diisSize + 1) = -one
328 [ # # ]: 0 : if(zDEBUG) write(std_out,*) "B vector:", diisCoeff(1:diisSize + 1)
329 [ # # # # : 0 : ABI_MALLOC(workMatrix,(diisSize + 1, diisSize + 1))
# # # # #
# # # ]
330 [ # # # # : 0 : ABI_MALLOC(workArray,((diisSize + 1) ** 2))
# # # # #
# # # ]
331 [ # # # # : 0 : 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 [ # # # # : 0 : 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 [ # # ]: 0 : if(zDEBUG)then
428 : 0 : write(std_out,*) "(A*X=B) A Matrix:", diisSize + 1,'x',diisSize + 1
429 [ # # ]: 0 : do ii=1,diisSize + 1
430 : 0 : write(std_out,*) workMatrix(1:diisSize + 1,ii)
431 : : end do
432 : :
433 : 0 : write(std_out,*) "(A*X=B) B Vector:", diisSize + 1
434 [ # # ]: 0 : do ii=1,diisSize + 1
435 : 0 : 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 : 0 : & workArray, (diisSize + 1) ** 2, info)
442 : :
443 [ # # ]: 0 : if (info /= 0) then
444 : 0 : write(std_out,*) "error solving DIIS matrix", info
445 [ # # ]: 0 : do ii=1,diisSize + 1
446 : 0 : write(std_out,*) workMatrix(1:diisSize + 1,ii)
447 : : end do
448 : : end if
449 : :
450 [ # # ]: 0 : if(zDEBUG)then
451 : 0 : write(std_out,*) "(A*X=B) X Vector:",diisSize+1
452 : 0 : suma=0.0
453 [ # # ]: 0 : do ii=1,diisSize+1
454 : 0 : write(std_out,*) ii,diisCoeff(ii)
455 : 0 : suma=suma+diisCoeff(ii)
456 : : end do
457 : :
458 : 0 : suma=suma-diisCoeff(diisSize+1)
459 : 0 : write(std_out,*) 'Sum of coefficients=',suma
460 : : end if
461 : :
462 [ # # ]: 0 : ABI_FREE(ipiv)
463 [ # # ]: 0 : ABI_FREE(workArray)
464 [ # # ]: 0 : ABI_FREE(workMatrix)
465 [ # # ]: 0 : 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 [ # # # # ]: 0 : xcart(:, :) = zero
475 [ # # # # ]: 0 : error(:, :, diisSize) = zero
476 [ # # ]: 0 : do ii = 1, diisSize, 1
477 : : xcart(:, :) = xcart(:, :) +&
478 [ # # # # ]: 0 : & xcart_hist(:, :, ii) * diisCoeff(ii)
479 : : error(:, :, diisSize) = error(:, :, diisSize)+&
480 [ # # # # ]: 0 : & fcart_hist(:, :, ii) * diisCoeff(ii)
481 : :
482 [ # # ]: 0 : if(zDEBUG)then
483 : 0 : write (std_out,*) 'Building new coordinates (ii):',ii
484 : 0 : write (std_out,*) 'diisCoeff(ii)',diisCoeff(ii)
485 : 0 : write (std_out,*) 'xcart_hist(:, :, ii)'
486 [ # # ]: 0 : do kk=1,ab_mover%natom
487 : 0 : write (std_out,*) xcart_hist(:,kk,ii)
488 : : end do
489 : 0 : write (std_out,*) 'fcart_hist(:, :, ii)'
490 [ # # ]: 0 : do kk=1,ab_mover%natom
491 : 0 : write (std_out,*) fcart_hist(:,kk,ii)
492 : : end do
493 : 0 : write (std_out,*) 'xcart:'
494 [ # # ]: 0 : do kk=1,ab_mover%natom
495 : 0 : write (std_out,*) xcart(:,kk)
496 : : end do
497 : 0 : write (std_out,*) 'error:'
498 [ # # ]: 0 : do kk=1,ab_mover%natom
499 : 0 : write (std_out,*) error(:,kk,diisSize)
500 : : end do
501 : : end if
502 : :
503 : : end do
504 : 0 : ABI_FREE(diisCoeff)
505 : :
506 [ # # ]: 0 : error_tmp(:)=RESHAPE( error(:,:,diisSize), (/ 3*ab_mover%natom /) )
507 [ # # ]: 0 : 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 : 0 : & ndim, error_tmp, 1, one, xcart_tmp, 1)
522 [ # # ]: 0 : 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 [ # # ]: 0 : 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 [ # # # # : 0 : & reshape(fcart_hist(:,:,diisSize-1), (/ ndim /)))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
533 : : end if
534 : :
535 [ # # ]: 0 : 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 : 0 : hist%ihist = abihist_findIndex(hist,+1)
543 : :
544 : : !Compute xred from xcart and rprimd
545 : 0 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
546 : :
547 : : !Fill the history with the variables
548 : : !xred, acell, rprimd, vel
549 : 0 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
550 : 0 : ihist_prev = abihist_findIndex(hist,-1)
551 [ # # # # ]: 0 : hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
552 : :
553 : : if (.false.) write(std_out,*) ntime
554 : :
555 [ # # ]: 0 : end subroutine pred_diisrelax
556 : : !!***
557 : :
558 : : end module m_pred_diisrelax
559 : : !!***
|
