TABLE OF CONTENTS


ABINIT/m_pred_verlet [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_verlet

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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_verlet
23 
24  use defs_basis
25  use m_abicore
26  use m_abimover
27  use m_abihist
28  use m_xfpack
29 
30  use m_geometry,       only : mkrdim, xcart2xred, xred2xcart, fcart2gred, metric
31 
32  implicit none
33 
34  private

ABINIT/pred_verlet [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_verlet

FUNCTION

 Ionmov predictors (6 & 7) Verlet algorithm

 IONMOV 6:
 Given a starting point xred that is a vector of length 3*natom
 (reduced nuclei coordinates), a velocity vector (in cartesian
 coordinates), and unit cell parameters (acell and rprimd -
 without velocities in the present implementation), the Verlet
 dynamics is performed, using the gradient of the energy
 (atomic forces and stresses) as calculated by the routine scfcv.

 Some atoms can be kept fixed, while the propagation of unit cell
 parameters is only performed if optcell/=0.
 No more than ab_mover%ntime steps are performed.
 The time step is governed by dtion, contained in ab_mover
 (coming from dtset).
 Returned quantities are xred, and eventually acell and rprimd (new ones!).

 IONMOV 7:
 Block every atom for which the scalar product of velocity and
 forces is negative, in order to reach the minimum.
 The convergence requirement on the atomic forces, ab_mover%tolmxf,
 allows an early exit.

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
 ionmov : (6 or 7) Specific kind of VERLET
 zDEBUG : if true print some debugging information

SIDE EFFECTS

 hist <type(abihist)> : History of positions,forces
                               acell, rprimd, stresses

SOURCE

 84 subroutine pred_verlet(ab_mover,hist,ionmov,itime,ntime,zDEBUG,iexit)
 85 
 86  implicit none
 87 
 88 !Arguments ------------------------------------
 89 !scalars
 90  type(abimover),intent(in) :: ab_mover
 91  type(abihist),intent(inout) :: hist
 92  integer,intent(in) :: itime
 93  integer,intent(in) :: ntime
 94  integer,intent(in) :: ionmov
 95  integer,intent(in) :: iexit
 96  logical,intent(in) :: zDEBUG
 97 
 98 !Local variables-------------------------------
 99 !scalars
100  integer  :: ii,istopped,jj,kk,ndim,nstopped
101  real(dp) :: amass_tot,etotal,diag,gr_avg,ekin_corr,scprod,taylor,ucvol0
102  real(dp),save :: ucvol,ucvol_next
103  character(len=500) :: message
104 !arrays
105  integer  :: stopped(ab_mover%natom)
106  real(dp) :: acell0(3),fcart(3,ab_mover%natom)
107  real(dp) :: gred_corrected(3,ab_mover%natom)
108  real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3), strten(6)
109  real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
110  real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
111  real(dp) :: vel(3,ab_mover%natom),vel_nexthalf(3,ab_mover%natom)
112  real(dp),save :: acell(3),acell_next(3)
113  real(dp),save :: rprimd(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3)
114  real(dp),allocatable,save :: hessin(:,:)
115  real(dp),allocatable,save :: vin(:),vin_prev(:),vin_next(:)
116  real(dp),allocatable,save :: vout(:),vout_prev(:),vel_prevhalf(:,:)
117 
118 !***************************************************************************
119 !Beginning of executable session
120 !***************************************************************************
121 
122  if(iexit/=0)then
123     ABI_SFREE(vin)
124     ABI_SFREE(vin_next)
125     ABI_SFREE(vout)
126     ABI_SFREE(vin_prev)
127     ABI_SFREE(vout_prev)
128     ABI_SFREE(hessin)
129     ABI_SFREE(vel_prevhalf)
130    return
131  end if
132 
133 !write(std_out,*) 'verlet 01'
134 !##########################################################
135 !### 01. Compute the dimension of vectors (ndim)
136 
137  ndim=3*ab_mover%natom
138  if(ab_mover%optcell==1) ndim=ndim+1
139  if(ab_mover%optcell==2 .or.&
140 & ab_mover%optcell==3) ndim=ndim+6
141  if(ab_mover%optcell>=4) ndim=ndim+3
142 
143 !write(std_out,*) 'verlet 02'
144 !##########################################################
145 !### 02. Allocate the vectors vin, vout and hessian matrix
146 
147 !Notice that vin, vout, etc could be allocated
148 !From a previous dataset with a different ndim
149  if(itime==1)then
150    ABI_SFREE(vin)
151    ABI_SFREE(vin_next)
152    ABI_SFREE(vout)
153    ABI_SFREE(vin_prev)
154    ABI_SFREE(vout_prev)
155    ABI_SFREE(hessin)
156    ABI_SFREE(vel_prevhalf)
157 
158    ABI_MALLOC(vin,(ndim))
159    ABI_MALLOC(vin_next,(ndim))
160    ABI_MALLOC(vout,(ndim))
161    ABI_MALLOC(vin_prev,(ndim))
162    ABI_MALLOC(vout_prev,(ndim))
163    ABI_MALLOC(hessin,(ndim,ndim))
164    ABI_MALLOC(vel_prevhalf,(3,ab_mover%natom))
165  end if
166 
167 !write(std_out,*) 'verlet 03'
168 !##########################################################
169 !### 03. Obtain the present values from the history
170 
171  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
172 
173  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
174  do ii=1,3
175    rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
176  end do
177 
178  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
179  fcart(:,:)  =hist%fcart(:,:,hist%ihist)
180  strten(:)  =hist%strten(:,hist%ihist)
181  vel(:,:)   =hist%vel(:,:,hist%ihist)
182  etotal     =hist%etot(hist%ihist)
183 
184  if(zDEBUG)then
185    write (std_out,*) 'fcart:'
186    do kk=1,ab_mover%natom
187      write (std_out,*) fcart(:,kk)
188    end do
189    write (std_out,*) 'vel:'
190    do kk=1,ab_mover%natom
191      write (std_out,*) vel(:,kk)
192    end do
193    write (std_out,*) 'strten:'
194    write (std_out,*) strten(1:3),ch10,strten(4:6)
195    write (std_out,*) 'etotal:'
196    write (std_out,*) etotal
197  end if
198 
199  acell0=acell ; ucvol0=ucvol
200 
201 !Get rid of mean force on whole unit cell, but only if no
202 !generalized constraints are in effect
203  call fcart2gred(fcart,gred_corrected,rprimd,ab_mover%natom)
204  if(ab_mover%nconeq==0)then
205    amass_tot=sum(ab_mover%amass(:))
206    do kk=1,3
207      if (kk/=3.or.ab_mover%jellslab==0) then
208        gr_avg=sum(gred_corrected(kk,:))/dble(ab_mover%natom)
209        gred_corrected(kk,:)=gred_corrected(kk,:)-gr_avg*ab_mover%amass(:)/amass_tot
210      end if
211    end do
212  end if
213 
214 !write(std_out,*) 'verlet 04'
215 !##########################################################
216 !### 04. Fill the vectors vin and vout
217 
218 !Initialize input vectors : first vin, then vout
219  call xfpack_x2vin(acell, ab_mover%natom, ndim,&
220 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd,&
221 & ab_mover%symrel, ucvol, ucvol0, vin, xred)
222  call xfpack_f2vout(gred_corrected, ab_mover%natom, ndim,&
223 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
224 & vout)
225 
226 !write(std_out,*) 'verlet 05'
227 !##########################################################
228 !### 05. Initialize or update the hessian matrix
229 
230 !Here, set up the matrix of transformation between forces and
231 !acceleration. Masses must be included here.
232 !Beside this feature, one could define
233 !a preconditioner, in which case it should
234 !be the inverse hessian, like in Broyden. This explains the
235 !name chosen for this transformation matrix. This would allow
236 !to find easily the optimal geometry with ionmov=7.
237 !The default, now implemented, corresponds to the identity matrix
238 !in cartesian coordinates, which makes use of metric tensor gmet
239 !in reduced coordinates.
240 
241 !Initialise the Hessian matrix using gmet
242  if (itime==1)then
243 !  Initialize inverse hessian with identity matrix
244 !  in cartesian coordinates, which makes use of metric tensor gmet
245 !  in reduced coordinates.
246    hessin(:,:)=zero
247    do ii=1,ab_mover%natom
248      do kk=1,3
249        do jj=1,3
250 !        Warning : implemented in reduced coordinates
251          if (ab_mover%iatfix(kk,ii)==0 .and.&
252 &         ab_mover%iatfix(jj,ii)==0 )then
253            hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)/ab_mover%amass(ii)
254          end if
255        end do
256      end do
257    end do
258    if(ab_mover%optcell/=0)then
259 !    These values might lead to too large changes in some cases ...
260      diag=ab_mover%strprecon*30.0_dp/ucvol
261      if(ab_mover%optcell==1) diag=diag/three
262      do ii=3*ab_mover%natom+1,ndim
263        hessin(ii,ii)=diag
264      end do
265    end if
266  end if
267 
268 !zDEBUG (vin,vout and hessin before prediction)
269  if(zDEBUG)then
270    write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
271    write(std_out,*) 'vin:'
272    do ii=1,ndim,3
273      if (ii+2<=ndim)then
274        write(std_out,*) ii,vin(ii:ii+2)
275      else
276        write(std_out,*) ii,vin(ii:ndim)
277      end if
278    end do
279    write(std_out,*) 'vout:'
280    do ii=1,ndim,3
281      if (ii+2<=ndim)then
282        write(std_out,*) ii,vout(ii:ii+2)
283      else
284        write(std_out,*) ii,vout(ii:ndim)
285      end if
286    end do
287    write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
288    do kk=1,ndim
289      do jj=1,ndim,3
290        if (jj+2<=ndim)then
291          write(std_out,*) jj,hessin(jj:jj+2,kk)
292        else
293          write(std_out,*) jj,hessin(jj:ndim,kk)
294        end if
295      end do
296    end do
297  end if
298 
299 !write(std_out,*) 'verlet 06'
300 !##########################################################
301 !### 06. Compute the next values
302 
303 !%%% VERLET ALGORITHM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
304 
305  if(zDEBUG)then
306    write(std_out,*) 'Shifted data ENTER'
307    write(std_out,*) 'acell(:)',acell(:)
308    write(std_out,*) 'rprimd(:,:)',rprimd(:,:)
309    write(std_out,*) 'ucvol',ucvol
310    write(std_out,*) 'vel_prevhalf(:,:)',vel_prevhalf(:,:)
311    write(std_out,*) 'vin_prev(:)',vin_prev(:)
312    write(std_out,*) 'vin(:)',vin(:)
313    write(std_out,*) 'xcart(:,:)',xcart(:,:)
314    write(std_out,*) 'xred(:,:)',xred(:,:)
315  end if
316 
317 !Compute next atomic coordinates and cell parameters, using
318 !Verlet algorithm
319 
320 !1. First propagate the position, without acceleration
321  if(itime/=1)then
322    vin_next(:)=2*vin(:)-vin_prev(:)
323    taylor=one
324  else
325 !  Initialisation : no vin_prev is available, but the ionic velocity
326 !  is available, in cartesian coordinates
327 !  Uses the velocity
328    xcart_next(:,:)=xcart(:,:)+ab_mover%dtion*vel(:,:)
329 !  Convert back to xred_next (reduced coordinates)
330    call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
331 !  Impose no change of acell, ucvol, rprim, and rprimd
332    acell_next(:)=acell(:)
333    ucvol_next=ucvol
334    rprim_next(:,:)=rprim(:,:)
335    rprimd_next(:,:)=rprimd(:,:)
336    write(std_out,*) 'ucvol',ucvol
337 !  Store all these next values in vin_next
338    call xfpack_x2vin(acell_next,ab_mover%natom,&
339 &   ndim,ab_mover%nsym,ab_mover%optcell,rprim_next,&
340 &   rprimd,ab_mover%symrel,ucvol_next,ucvol0,&
341 &   vin_next,xred_next)
342    taylor=half
343  end if
344 
345 !2. Now, take into account the acceleration
346  do ii=1,ndim
347 !  Note the minus sign: the forces are minus the gradients,
348 !  contained in vout.
349    vin_next(:)=vin_next(:)-ab_mover%dtion**2*hessin(:,ii)*&
350 &   vout(ii)*taylor
351  end do
352 
353 !3. Implement fixing of atoms : put back old values for fixed
354 !components
355  do kk=1,ab_mover%natom
356    do jj=1,3
357 !    Warning : implemented in reduced coordinates
358      if (ab_mover%iatfix(jj,kk) == 1) then
359        vin_next(jj+(kk-1)*3)=vin(jj+(kk-1)*3)
360      end if
361    end do
362  end do
363 
364 !4. Now, compute the velocity at the next half-step
365 !Get xred_next, and eventually acell_next, ucvol_next, rprim_next and
366 !rprimd_next, from vin_next
367  call xfpack_vin2x(acell_next,acell0,ab_mover%natom,ndim,&
368 & ab_mover%nsym,ab_mover%optcell,rprim_next,rprimd,&
369 & ab_mover%symrel,ucvol_next,ucvol0,vin_next,xred_next)
370  if(ab_mover%optcell/=0)then
371    call mkrdim(acell_next,rprim_next,rprimd_next)
372    call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next)
373  else
374    rprimd_next(:,:)=rprimd(:,:)
375  end if
376 !Convert input xred_next (reduced coordinates) to
377 !xcart_next (cartesian)
378  call xred2xcart(ab_mover%natom,rprimd_next,xcart_next,xred_next)
379 !Compute the velocity at half of the new step
380  vel_nexthalf(:,:)=(xcart_next(:,:)-xcart(:,:))/ab_mover%dtion
381 
382 !5. If needed, compute the velocity at present position
383  if(itime/=1)then
384    vel(:,:)=(vel_nexthalf(:,:)+vel_prevhalf(:,:))*0.5_dp
385  end if
386 
387 !%%% VERLET ALGORITHM BLOCKING ATOMS %%%%%%%%%%%%%%%%%%%%%%
388 
389 !Here, stop the atoms for which the scalar product of velocity
390 !and force is negative, and recompute the kinetic energy.
391  if(ionmov==7)then
392    stopped(:)=0
393    do ii=1,ab_mover%natom
394      scprod=fcart(1,ii)*vel(1,ii)+&
395 &     fcart(2,ii)*vel(2,ii)+&
396 &     fcart(3,ii)*vel(3,ii)
397      if(scprod<0.0_dp .and. itime/=1)then
398        stopped(ii)=1
399        write(std_out,*) 'Stopped atom',ii
400 !      Shift the velocities of the previous half-step and current
401 !      half-step, so that the acceleration is correct but the
402 !      present velocity vanishes.
403        vel_prevhalf(:,ii)=vel_prevhalf(:,ii)-vel(:,ii)
404        vel_nexthalf(:,ii)=vel_nexthalf(:,ii)-vel(:,ii)
405        vel(:,ii)=0.0_dp
406        xcart_next(:,ii)=xcart(:,ii)+ab_mover%dtion*vel_nexthalf(:,ii)
407      end if
408    end do
409 
410    if(zDEBUG)then
411      write (std_out,*) 'fcart:'
412      do kk=1,ab_mover%natom
413        write (std_out,*) fcart(:,kk)
414      end do
415      write (std_out,*) 'vel_prevhalf:'
416      do kk=1,ab_mover%natom
417        write (std_out,*) vel_prevhalf(:,kk)
418      end do
419      write (std_out,*) 'vel_nexthalf:'
420      do kk=1,ab_mover%natom
421        write (std_out,*) vel_nexthalf(:,kk)
422      end do
423      write (std_out,*) 'vel:'
424      do kk=1,ab_mover%natom
425        write (std_out,*) vel(:,kk)
426      end do
427      write (std_out,*) 'xcart_next:'
428      do kk=1,ab_mover%natom
429        write (std_out,*) xcart_next(:,kk)
430      end do
431    end if
432 
433 !  Establish a list of stopped atoms
434    nstopped=sum(stopped(:))
435 
436    if(nstopped/=0)then
437      write(message,'(a)') ' List of stopped atoms (ionmov=7) :'
438      call wrtout(ab_out,message,'COLL')
439      istopped=1
440      do ii=1,ab_mover%natom
441        if(stopped(ii)==1)then
442          stopped(istopped)=ii
443          istopped=istopped+1
444        end if
445      end do
446      do ii=1,nstopped,16
447        write(message, '(16i4)' ) stopped(ii:min(ii+15,nstopped))
448        call wrtout(ab_out,message,'COLL')
449      end do
450 !    Now, compute the corrected kinetic energy
451 !    Generate xred_next from xcart_next
452      call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next)
453 !    Store xred_next, and eventual acell_next and rprim_next in vin
454      call xfpack_x2vin(acell_next,&
455 &     ab_mover%natom,ndim,ab_mover%nsym,ab_mover%optcell,&
456 &     rprim_next,rprimd,&
457 &     ab_mover%symrel,ucvol_next,ucvol0,vin_next,xred_next)
458 
459      ekin_corr=0.0_dp
460      do ii=1,ab_mover%natom
461        do jj=1,3
462 !        Warning : the fixing of atomis is implemented in reduced
463 !        coordinates, so that this expression is wrong
464          if (ab_mover%iatfix(jj,ii) == 0) then
465            ekin_corr=ekin_corr+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2
466          end if
467        end do
468      end do
469 !    End of test nstopped/=0
470    end if
471 
472 !  End of test ionmov==7
473  end if
474 
475 !write(std_out,*) 'verlet 07'
476 !##########################################################
477 !### 07. Shift the data from next values to the present
478 
479 !acell(:)=acell_next(:)
480 !rprim(:,:)=rprim_next(:,:)
481 !ucvol=ucvol_next
482  vel_prevhalf(:,:)=vel_nexthalf(:,:)
483  vin_prev(:)=vin(:)
484  vin(:)=vin_next(:)
485  xcart(:,:)=xcart_next(:,:)
486  xred(:,:)=xred_next(:,:)
487 
488  write(std_out,*) 'Shifted data EXIT'
489  write(std_out,*) 'acell(:)',acell(:)
490  write(std_out,*) 'rprim(:,:)',rprim(:,:)
491  write(std_out,*) 'rprimd(:,:)',rprimd(:,:)
492  write(std_out,*) 'ucvol',ucvol
493  write(std_out,*) 'vel_prevhalf(:,:)',vel_prevhalf(:,:)
494  write(std_out,*) 'vin_prev(:)',vin_prev(:)
495  write(std_out,*) 'vin(:)',vin(:)
496  write(std_out,*) 'xred(:,:)',xred(:,:)
497 
498 
499 !write(std_out,*) 'verlet 08'
500 !##########################################################
501 !### 08. Update the history with the prediction
502 
503 !Increase indexes
504  hist%ihist=abihist_findIndex(hist,+1)
505 
506 !Fill the history with the variables
507 !xred, acell, rprimd, vel
508  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
509  hist%vel(:,:,hist%ihist)=vel(:,:)
510  hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
511 
512  if(zDEBUG)then
513    write (std_out,*) 'vel:'
514    do kk=1,ab_mover%natom
515      write (std_out,*) vel(:,kk)
516    end do
517  end if
518 
519  if (.false.) write(std_out,*) ntime
520 
521 end subroutine pred_verlet