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 : !!***
|
