TABLE OF CONTENTS


ABINIT/m_pred_velverlet [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_velverlet

FUNCTION

COPYRIGHT

  Copyright (C) 2017-2018 ABINIT group (SPr)
  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_velverlet
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_abimover
33  use m_abihist
34 
35  use m_geometry,  only : xcart2xred, xred2xcart
36 
37  implicit none
38 
39  private

ABINIT/pred_velverlet [ Functions ]

[ Top ] [ Functions ]

NAME

  pred_velverlet

FUNCTION

  Velocity Verlet (VV) predictor of ionic positions (ionmov = 24).
  In constrast to Verlet algorithm, Velocity Verlet is a
  SYMPLECTIC integrator (for small enough step sizes "dtion",
  it better conserves the total energy and time-reversibility).
  VV is a second order integration scheme that requires a single
  evaluatoin of forces per time step. These properties make VV
  a good candidate integrator for use in Hybrid Monte Carlo simulation scheme.

INPUTS

  ab_mover =  Data structure containing information about
              input variables related to MD, e.g dtion, masses, etc.
  hist     =  history of ionic positions, forces,
  itime    =  index of current time step
  ntime    =  total number of time steps
  zDEBUG   =  flag indicating whether to print debug info
  iexit    =  flag indicating finilization of mover loop
  hmcflag  =  optional argument indicating whether the predictor is called from the Hybrid Monte Carlo (HMC) routine
  icycle   =  if hmcflag==1, then icycle providing information about number of HMC cycle is needed
  ncycle   =  if hmcflag==1, then ncycle provides the total number of cycles within one HMC iteration

OUTPUT

  hist =  history of ionic positions, forces etc. is updated

SIDE EFFECTS

NOTES

 This routine can be used either to simulate NVE molecular dynamics (ionmov = 24) or
 is called from pred_hmc routine (ionmov = 25) to perform updates of ionic positions
 in Hybrid Monte Carlo iterations.

PARENTS

      mover,pred_hmc

CHILDREN

      hist2var,var2hist,xcart2xred,xred2xcart

SOURCE

 92 subroutine pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,hmcflag,icycle,ncycle)
 93 
 94 
 95 !This section has been created automatically by the script Abilint (TD).
 96 !Do not modify the following lines by hand.
 97 #undef ABI_FUNC
 98 #define ABI_FUNC 'pred_velverlet'
 99 !End of the abilint section
100 
101  implicit none
102 
103 !Arguments ------------------------------------
104  type(abimover),intent(in)   :: ab_mover
105  type(abihist),intent(inout) :: hist
106  integer,intent(in) :: itime
107  integer,intent(in) :: ntime
108  integer,intent(in) :: iexit
109  logical,intent(in) :: zDEBUG
110  integer,intent(in),optional :: hmcflag
111  integer,intent(in),optional :: icycle
112  integer,intent(in),optional :: ncycle
113 
114 !Local variables-------------------------------
115 
116  integer  :: ii,jj                                                              ! dummy integers for loop indexes
117  real(dp) :: epot,ekin,ekin_tmp                                                 ! potential (electronic), kinetic (ionic) energies
118  real(dp) :: xcart(3,ab_mover%natom)                                            ! Cartesian coordinates of all ions
119  real(dp) :: xred(3,ab_mover%natom)                                             ! reduced coordinates of all ions
120  real(dp) :: vel(3,ab_mover%natom)                                              ! ionic velocities in Cartesian coordinates
121  real(dp) :: fcart(3,ab_mover%natom),fred(3,ab_mover%natom)                     ! forces, Cartesian and reduced coordinates
122  real(dp) :: factor                                                             ! factor, indicating change of time step at last iteration
123  integer :: hmcflag_
124  integer :: icycle_
125  integer :: ncycle_
126 
127  real(dp) :: acell(3)                                                           ! lattice parameters
128  real(dp) :: rprimd(3,3)                                                        ! lattice vectors
129  real(dp),allocatable,save :: vel_prev(:,:)                                     ! velocities at the end of each time step (half time step ahead of coordinates)
130 
131 !***************************************************************************
132 !Beginning of executable session
133 !***************************************************************************
134 
135  DBG_ENTER("COLL")
136 
137 ! if (option/=1 .and. option/=2 ) then
138 !   write(msg,'(3a,i0)')&
139 !&   'The argument option should be 1 or 2,',ch10,&
140 !&   'however, option=',option
141 !   MSG_BUG(msg)
142 ! end if
143 !
144 ! if (sizein<1) then
145 !   write(msg,'(3a,i0)')&
146 !&   'The argument sizein should be a positive number,',ch10,&
147 !&   'however, sizein=',sizein
148 !   MSG_ERROR(msg)
149 ! end if
150 
151  DBG_EXIT("COLL")
152 
153 
154  hmcflag_=0
155  if(present(hmcflag))then
156    hmcflag_=hmcflag
157  end if
158 
159  icycle_ =0
160  if(present(icycle))then
161    icycle_=icycle
162  end if
163 
164  ncycle_ =0
165  if(present(ncycle))then
166    ncycle_=ncycle
167  end if
168 
169 
170  if(iexit/=0)then
171    if (allocated(vel_prev))  then
172      ABI_DEALLOCATE(vel_prev)
173    end if
174    return
175  end if
176 
177  if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then
178    if (allocated(vel_prev))  then
179      ABI_DEALLOCATE(vel_prev)
180    end if
181    ABI_ALLOCATE(vel_prev,(3,ab_mover%natom))
182  end if
183 
184 
185 
186  ! Start preparation for velocity verlet, get information about current ionic positions, forces, velocities, etc.
187 
188  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
189  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
190  fcart(:,:) = hist%fcart(:,:,hist%ihist)              ! forces in Cartesian coordinates
191  vel(:,:)   = hist%vel(:,:,hist%ihist)                ! velocities of all ions, not needed in reality
192  epot       = hist%etot(hist%ihist)                   ! electronic sub-system energy, not needed
193  ekin       = hist%ekin(hist%ihist)                   ! kinetic energy, not needed
194 
195 
196  if(zDEBUG)then
197    write (std_out,*) 'velverlet step',itime
198    write (std_out,*) 'fcart:'
199    do ii=1,3
200      write (std_out,*) fcart(ii,:)
201    end do
202    write (std_out,*) 'fred:'
203    do ii=1,3
204      write (std_out,*) fred(ii,:)
205    end do
206    write (std_out,*) 'xcart:'
207    do ii=1,3
208      write (std_out,*) xcart(ii,:)
209    end do
210    write (std_out,*) 'vel:'
211    do ii=1,3
212      write (std_out,*) vel(ii,:)
213    end do
214  end if
215 
216 
217  if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then
218 
219    !the following breakdown of single time step in two halfs is needed for initialization.
220    !half step velocities "vel_prev" are saved to be used in the next iteration
221    !the velocities "vel" are only used to estimate kinetic energy at correct time instances
222    do ii=1,ab_mover%natom ! propagate velocities half time step forward
223      do jj=1,3
224        vel_prev(jj,ii) = vel(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
225      end do
226    end do
227 
228    ! propagate velocities half time step forward
229    do ii=1,ab_mover%natom
230      do jj=1,3
231        vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
232      end do
233    end do
234    ! use half-step behind velocity values to propagate coordinates one time step forward!!!!
235    do ii=1,ab_mover%natom
236      do jj=1,3
237        xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii)
238      end do
239    end do
240    ! now, at this 1st iteration, "vel_prev" correspond to a time instance half-step behind
241    ! that of "xcart"
242 
243  else
244 
245    !at this moment "vel_prev" is behind "xcart" by half of a time step
246    !(saved from the previous iteration) and these are the velocity values to be propagated
247    !using forces that are evaluated at the same time instance as xcart
248    do ii=1,ab_mover%natom ! propagate velocities one time step forward
249      do jj=1,3
250        vel_prev(jj,ii) = vel_prev(jj,ii) + ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
251      end do
252    end do
253    !now, the "vel_prev" velocities are half of a time step ahead and can be used to propagate xcart
254 
255    !if((hmcflag_==0.and.itime==ntime-1).or.(hmcflag_==1.and.icycle_==ncycle_-1))then
256    !  factor=0.5_dp
257    !else
258    !  factor=one
259    !end if
260 
261    do ii=1,ab_mover%natom ! propagate coordinates
262      do jj=1,3
263        xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii)
264       !xcart(jj,ii) = xcart(jj,ii) + factor*ab_mover%dtion*vel_prev(jj,ii)
265      end do
266    end do
267    !to estimate kinetic energy at the same time instance as the potential (electronic sub-system) energy
268    !propagate "vel" another half-time forward (these values are not to be used in the next time-step)
269    do ii=1,ab_mover%natom ! propagate velocities half time step forward
270      do jj=1,3
271        vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
272       !vel(jj,ii) = vel_prev(jj,ii) +(factor-0.5_dp) * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
273      end do
274    end do
275 
276    ekin=0.0
277    do ii=1,ab_mover%natom
278      do jj=1,3
279        ekin=ekin+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2
280      end do
281    end do
282    !ekin_tmp=0.0
283    !do ii=1,ab_mover%natom
284    !  do jj=1,3
285    !    ekin_tmp=ekin_tmp+0.5_dp*ab_mover%amass(ii)*vel_prev(jj,ii)**2
286    !  end do
287    !end do
288    !write(238,*) itime,icycle,ekin_tmp,ekin,epot,factor
289 
290  end if
291 
292  !Convert new xcart to xred to set correct output values
293  !Update the history with the new coordinates, velocities, etc.
294 
295  !Increase indexes
296  hist%ihist=abihist_findIndex(hist,+1)
297 
298  call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
299  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
300 
301  hist%vel(:,:,hist%ihist)=vel(:,:)
302  hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
303 
304 end subroutine pred_velverlet