TABLE OF CONTENTS


ABINIT/m_pred_diisrelax [ Modules ]

[ Top ] [ 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 ]

[ Top ] [ 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