Branch data 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 : : !!***
|
