Line data Source code
1 : !!****m* ABINIT/m_pred_velverlet
2 : !! NAME
3 : !! m_pred_velverlet
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2017-2026 ABINIT group (SPr)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_pred_velverlet
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_abimover
28 : use m_abihist
29 :
30 : use m_geometry, only : xcart2xred, xred2xcart
31 :
32 : implicit none
33 :
34 : private
35 : !!***
36 :
37 : public :: pred_velverlet
38 : !!***
39 :
40 : contains
41 : !!***
42 :
43 : !!****f* ABINIT/pred_velverlet
44 : !! NAME
45 : !! pred_velverlet
46 : !!
47 : !! FUNCTION
48 : !! Velocity Verlet (VV) predictor of ionic positions (ionmov = 24).
49 : !! In constrast to Verlet algorithm, Velocity Verlet is a
50 : !! SYMPLECTIC integrator (for small enough step sizes "dtion",
51 : !! it better conserves the total energy and time-reversibility).
52 : !! VV is a second order integration scheme that requires a single
53 : !! evaluatoin of forces per time step. These properties make VV
54 : !! a good candidate integrator for use in Hybrid Monte Carlo simulation scheme.
55 : !!
56 : !! INPUTS
57 : !! ab_mover = Data structure containing information about
58 : !! input variables related to MD, e.g dtion, masses, etc.
59 : !! hist = history of ionic positions, forces,
60 : !! itime = index of current time step
61 : !! ntime = total number of time steps
62 : !! zDEBUG = flag indicating whether to print debug info
63 : !! iexit = flag indicating finilization of mover loop
64 : !! hmcflag = optional argument indicating whether the predictor is called from the Hybrid Monte Carlo (HMC) routine
65 : !! icycle = if hmcflag==1, then icycle providing information about number of HMC cycle is needed
66 : !! ncycle = if hmcflag==1, then ncycle provides the total number of cycles within one HMC iteration
67 : !!
68 : !! OUTPUT
69 : !! hist = history of ionic positions, forces etc. is updated
70 : !!
71 : !! SIDE EFFECTS
72 : !!
73 : !! NOTES
74 : !!
75 : !! This routine can be used either to simulate NVE molecular dynamics (ionmov = 24) or
76 : !! is called from pred_hmc routine (ionmov = 25) to perform updates of ionic positions
77 : !! in Hybrid Monte Carlo iterations.
78 : !!
79 : !! SOURCE
80 :
81 10 : subroutine pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,hmcflag,icycle,ncycle)
82 :
83 : !Arguments ------------------------------------
84 : type(abimover),intent(in) :: ab_mover
85 : type(abihist),intent(inout) :: hist
86 : integer,intent(in) :: itime
87 : integer,intent(in) :: ntime
88 : integer,intent(in) :: iexit
89 : logical,intent(in) :: zDEBUG
90 : integer,intent(in),optional :: hmcflag
91 : integer,intent(in),optional :: icycle
92 : integer,intent(in),optional :: ncycle
93 :
94 : !Local variables-------------------------------
95 :
96 : integer :: ii,jj ! dummy integers for loop indexes
97 : real(dp) :: epot,ekin !,ekin_tmp ! potential (electronic), kinetic (ionic) energies
98 20 : real(dp) :: xcart(3,ab_mover%natom) ! Cartesian coordinates of all ions
99 20 : real(dp) :: xred(3,ab_mover%natom) ! reduced coordinates of all ions
100 20 : real(dp) :: vel(3,ab_mover%natom) ! ionic velocities in Cartesian coordinates
101 1 : real(dp) :: fcart(3,ab_mover%natom),gred(3,ab_mover%natom) ! cartesian forces, and gradient in reduced coordinates
102 : !real(dp) :: factor ! factor, indicating change of time step at last iteration
103 : integer :: hmcflag_
104 : integer :: icycle_
105 : integer :: ncycle_
106 :
107 : real(dp) :: acell(3) ! lattice parameters
108 : real(dp) :: rprimd(3,3) ! lattice vectors
109 : real(dp),allocatable,save :: vel_prev(:,:) ! velocities at the end of each time step (half time step ahead of coordinates)
110 :
111 : !***************************************************************************
112 : !Beginning of executable session
113 : !***************************************************************************
114 :
115 : DBG_ENTER("COLL")
116 :
117 : ABI_UNUSED((/ntime/))
118 :
119 10 : hmcflag_=0
120 10 : if(present(hmcflag))then
121 0 : hmcflag_=hmcflag
122 : end if
123 :
124 10 : icycle_ =0
125 10 : if(present(icycle))then
126 0 : icycle_=icycle
127 : end if
128 :
129 10 : ncycle_ =0
130 : if(present(ncycle))then
131 : ncycle_=ncycle
132 : end if
133 :
134 :
135 10 : if(iexit/=0)then
136 1 : ABI_SFREE(vel_prev)
137 : return
138 : end if
139 :
140 9 : if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then
141 1 : ABI_SFREE(vel_prev)
142 3 : ABI_MALLOC(vel_prev,(3,ab_mover%natom))
143 : end if
144 :
145 : ! Start preparation for velocity verlet, get information about current ionic positions, forces, velocities, etc.
146 :
147 9 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
148 9 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
149 81 : fcart(:,:) = hist%fcart(:,:,hist%ihist) ! forces in Cartesian coordinates
150 81 : vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality
151 9 : epot = hist%etot(hist%ihist) ! electronic sub-system energy, not needed
152 9 : ekin = hist%ekin(hist%ihist) ! kinetic energy, not needed
153 :
154 :
155 9 : if(zDEBUG)then
156 0 : write (std_out,*) 'velverlet step',itime
157 0 : write (std_out,*) 'fcart:'
158 0 : do ii=1,3
159 0 : write (std_out,*) fcart(ii,:)
160 : end do
161 0 : write (std_out,*) 'gred:'
162 0 : do ii=1,3
163 0 : write (std_out,*) gred(ii,:)
164 : end do
165 0 : write (std_out,*) 'xcart:'
166 0 : do ii=1,3
167 0 : write (std_out,*) xcart(ii,:)
168 : end do
169 0 : write (std_out,*) 'vel:'
170 0 : do ii=1,3
171 0 : write (std_out,*) vel(ii,:)
172 : end do
173 : end if
174 :
175 :
176 9 : if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then
177 :
178 : !the following breakdown of single time step in two halfs is needed for initialization.
179 : !half step velocities "vel_prev" are saved to be used in the next iteration
180 : !the velocities "vel" are only used to estimate kinetic energy at correct time instances
181 3 : do ii=1,ab_mover%natom ! propagate velocities half time step forward
182 9 : do jj=1,3
183 8 : vel_prev(jj,ii) = vel(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
184 : end do
185 : end do
186 :
187 : ! propagate velocities half time step forward
188 3 : do ii=1,ab_mover%natom
189 9 : do jj=1,3
190 8 : vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
191 : end do
192 : end do
193 : ! use half-step behind velocity values to propagate coordinates one time step forward!!!!
194 3 : do ii=1,ab_mover%natom
195 9 : do jj=1,3
196 8 : xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii)
197 : end do
198 : end do
199 : ! now, at this 1st iteration, "vel_prev" correspond to a time instance half-step behind
200 : ! that of "xcart"
201 :
202 : else
203 :
204 : !at this moment "vel_prev" is behind "xcart" by half of a time step
205 : !(saved from the previous iteration) and these are the velocity values to be propagated
206 : !using forces that are evaluated at the same time instance as xcart
207 24 : do ii=1,ab_mover%natom ! propagate velocities one time step forward
208 72 : do jj=1,3
209 64 : vel_prev(jj,ii) = vel_prev(jj,ii) + ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
210 : end do
211 : end do
212 : !now, the "vel_prev" velocities are half of a time step ahead and can be used to propagate xcart
213 :
214 : !if((hmcflag_==0.and.itime==ntime-1).or.(hmcflag_==1.and.icycle_==ncycle_-1))then
215 : ! factor=0.5_dp
216 : !else
217 : ! factor=one
218 : !end if
219 :
220 24 : do ii=1,ab_mover%natom ! propagate coordinates
221 72 : do jj=1,3
222 64 : xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii)
223 : !xcart(jj,ii) = xcart(jj,ii) + factor*ab_mover%dtion*vel_prev(jj,ii)
224 : end do
225 : end do
226 : !to estimate kinetic energy at the same time instance as the potential (electronic sub-system) energy
227 : !propagate "vel" another half-time forward (these values are not to be used in the next time-step)
228 24 : do ii=1,ab_mover%natom ! propagate velocities half time step forward
229 72 : do jj=1,3
230 64 : vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
231 : !vel(jj,ii) = vel_prev(jj,ii) +(factor-0.5_dp) * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii)
232 : end do
233 : end do
234 :
235 : ekin=0.0
236 : do ii=1,ab_mover%natom
237 : do jj=1,3
238 : ekin=ekin+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2
239 : end do
240 : end do
241 : !ekin_tmp=0.0
242 : !do ii=1,ab_mover%natom
243 : ! do jj=1,3
244 : ! ekin_tmp=ekin_tmp+0.5_dp*ab_mover%amass(ii)*vel_prev(jj,ii)**2
245 : ! end do
246 : !end do
247 : !write(238,*) itime,icycle,ekin_tmp,ekin,epot,factor
248 :
249 : end if
250 :
251 : !Convert new xcart to xred to set correct output values
252 : !Update the history with the new coordinates, velocities, etc.
253 :
254 : !Increase indexes
255 9 : hist%ihist=abihist_findIndex(hist,+1)
256 :
257 9 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
258 9 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
259 :
260 81 : hist%vel(:,:,hist%ihist)=vel(:,:)
261 9 : hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
262 :
263 : DBG_EXIT("COLL")
264 :
265 : end subroutine pred_velverlet
266 : !!***
267 :
268 : end module m_pred_velverlet
269 : !!***
|
