TABLE OF CONTENTS


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.

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 .
 For the initials of contributors,
 see ~abinit/doc/developers/contributors.txt .

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

OUTPUT

SIDE EFFECTS

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

NOTES

PARENTS

      mover

CHILDREN

      fcart2fred,hist2var,metric,mkrdim,var2hist,wrtout,xcart2xred
      xfpack_f2vout,xfpack_vin2x,xfpack_x2vin,xred2xcart

SOURCE

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