TABLE OF CONTENTS


ABINIT/m_pred_verlet [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_verlet

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_verlet
28 
29  use defs_basis
30  use m_abicore
31  use m_abimover
32  use m_abihist
33  use m_xfpack
34 
35  use m_geometry,       only : mkrdim, xcart2xred, xred2xcart, fcart2fred, metric
36 
37  implicit none
38 
39  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

PARENTS

      mover

CHILDREN

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

SOURCE

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