Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_moldyn
2 : : !! NAME
3 : : !! m_pred_moldyn
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, SE)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_pred_moldyn
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_abimover
26 : : use m_abihist
27 : :
28 : : use m_geometry, only : xcart2xred, xred2xcart
29 : : use m_predtk, only : fdtion
30 : :
31 : : implicit none
32 : :
33 : : private
34 : : !!***
35 : :
36 : : public :: pred_moldyn
37 : : !!***
38 : :
39 : : contains
40 : : !!***
41 : :
42 : : !!****f* ABINIT/pred_moldyn
43 : : !! NAME
44 : : !! pred_moldyn
45 : : !!
46 : : !! FUNCTION
47 : : !! Ionmov predictor (1) Molecular dynamics
48 : : !!
49 : : !! Molecular dynamics, with or without viscous damping
50 : : !! This function should be called after the call to scfcv
51 : : !! Updates positions, velocities and forces
52 : : !!
53 : : !! INPUTS
54 : : !! ab_mover<type abimover>=Subset of dtset only related with
55 : : !! | movement of ions and acell, contains:
56 : : !! | dtion: Time step
57 : : !! ! natom: Number of atoms
58 : : !! | vis: viscosity
59 : : !! | iatfix: Index of atoms and directions fixed
60 : : !! | amass: Mass of ions
61 : : !! icycle: Index of the internal cycle inside a time step (itime)
62 : : !! itime: Index of time iteration
63 : : !! zDEBUG : if true print some debugging information
64 : : !!
65 : : !! OUTPUT
66 : : !!
67 : : !! SIDE EFFECTS
68 : : !! hist<type abihist>=Historical record of positions, forces,
69 : : !! stresses, cell and energies,
70 : : !!
71 : : !! ncycle: Number of cycles of a particular time step
72 : : !!
73 : : !! NOTES
74 : : !! * This routine is a predictor, it only produces new positions
75 : : !! to be computed in the next iteration, this routine should
76 : : !! produce not output at all
77 : : !! * ncycle changes from 4 for the first iteration (itime==1) to 1 for (itime>1)
78 : : !! * The arrays vec_tmp1 and vec_tmp2 are triky, they are use with
79 : : !! different meanings, during the initialization they contains
80 : : !! working positions and velocities that acumulated produce the
81 : : !! first positions of itime=1, for itime>1 they will contain
82 : : !! positions in 2 previous steps, those values are different
83 : : !! from the values store in the history, thats the reason why
84 : : !! we cannot simply use hist%xred to obtain those positions.
85 : : !!
86 : : !! SOURCE
87 : :
88 : 73 : subroutine pred_moldyn(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit)
89 : :
90 : : !Arguments ------------------------------------
91 : : !scalars
92 : : type(abimover),intent(in) :: ab_mover
93 : : type(abihist),intent(inout),target :: hist
94 : : integer,intent(in) :: icycle
95 : : integer,intent(inout) :: ncycle
96 : : integer,intent(in) :: itime
97 : : integer,intent(in) :: ntime
98 : : integer,intent(in) :: iexit
99 : : logical,intent(in) :: zDEBUG
100 : :
101 : : !Local variables-------------------------------
102 : : !scalars
103 : : integer :: kk,jj,ihist,ihist_next,ihist_prev,ihist_prev2
104 : : integer :: ihist_prev4,ihist_prev5
105 : : real(dp) :: aa,alfa,bb,cc,x0,xm,em,vis,dx,dv
106 : : real(dp) :: fcart,fprev,fprev2
107 : : real(dp) :: xc
108 : : real(dp) :: vel,vnow,xnow,vprev
109 : : real(dp),save :: hh,time
110 : : !arrays
111 : : real(dp) :: acell(3),rprimd(3,3)
112 : 6 : real(dp) :: xred(3,ab_mover%natom)
113 : 73 : real(dp),allocatable :: xcart(:,:),xcart_prev(:,:)
114 : : real(dp),save,allocatable :: vec_tmp1(:,:)
115 : : real(dp),save,allocatable :: vec_tmp2(:,:)
116 : 73 : real(dp), contiguous, pointer :: vel_cur(:,:),vel_next(:,:)
117 : 73 : real(dp),pointer :: fcart_cur(:,:),fcart_prev(:,:),fcart_prev2(:,:)
118 : :
119 : : !***************************************************************************
120 : : !Beginning of executable session
121 : : !***************************************************************************
122 : :
123 : 73 : if(iexit/=0)then
124 [ + - ]: 6 : if(allocated(vec_tmp1)) then
125 : 6 : ABI_FREE(vec_tmp1)
126 : : end if
127 [ + - ]: 6 : if(allocated(vec_tmp2)) then
128 : 6 : ABI_FREE(vec_tmp2)
129 : : end if
130 : : return
131 : : end if
132 : :
133 : 67 : vis= ab_mover%vis
134 : : !Just to avoid warnings of uninitialized variables
135 : : fprev=0.0_dp
136 : 67 : fprev=0.0_dp
137 : 67 : fprev2=0.0_dp
138 : 67 : vnow=0.0_dp
139 : 67 : vprev=0.0_dp
140 : 67 : xnow=0.0_dp
141 : :
142 : : !Those arrays contains intermediary results used with
143 : : !different meanings during the different time steps
144 : : !We need to preserv the allocation status, this is the
145 : : !reason to be 'SAVE'
146 [ + + + + ]: 67 : if (itime==1.and.icycle==1)then
147 [ - + ]: 6 : if(allocated(vec_tmp1)) then
148 : 0 : ABI_FREE(vec_tmp1)
149 : : end if
150 [ - + ]: 6 : if(allocated(vec_tmp2)) then
151 : 0 : ABI_FREE(vec_tmp2)
152 : : end if
153 [ - + + - : 18 : ABI_MALLOC(vec_tmp1,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
154 [ - + + - : 12 : ABI_MALLOC(vec_tmp2,(3,ab_mover%natom))
- + - + -
+ ]
155 : : end if
156 : :
157 : : !write(std_out,*) '00'
158 : : !##########################################################
159 : : !### 00. Copy from the history to the variables
160 : :
161 : 67 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
162 : :
163 [ - + + - : 201 : ABI_MALLOC(xcart,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
164 : 67 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
165 : :
166 [ + + ]: 67 : if (itime==1.or.itime==2)then
167 [ - + + - : 90 : ABI_MALLOC(xcart_prev,(3,ab_mover%natom))
+ - + - -
+ - + ]
168 : 30 : call xred2xcart(ab_mover%natom,rprimd,xcart_prev,hist%xred(:,:,1))
169 : : end if
170 : :
171 : 67 : ihist = abihist_findIndex(hist, 0)
172 : 67 : ihist_prev = abihist_findIndex(hist,-1)
173 : 67 : ihist_prev2 = abihist_findIndex(hist,-2)
174 : 67 : ihist_prev4 = abihist_findIndex(hist,-4)
175 : 67 : ihist_prev5 = abihist_findIndex(hist,-5)
176 : 67 : ihist_next = abihist_findIndex(hist,+1)
177 : :
178 : 67 : fcart_cur => hist%fcart(:,:,ihist)
179 [ + + ]: 67 : if (itime==2) fcart_prev => hist%fcart(:,:,ihist_prev4)
180 [ + + ]: 67 : if (itime==3) fcart_prev2 => hist%fcart(:,:,ihist_prev5)
181 [ + + + + ]: 67 : if (itime >2.or. icycle>=2)fcart_prev => hist%fcart(:,:,ihist_prev)
182 [ + + + + ]: 67 : if (itime >3.or. icycle>=3)fcart_prev2 => hist%fcart(:,:,ihist_prev2)
183 : :
184 : 67 : vel_cur => hist%vel(:,:,ihist)
185 : 67 : vel_next => hist%vel(:,:,ihist_next)
186 : :
187 : : !write(std_out,*) '01'
188 : : !##########################################################
189 : : !### 01. Get or compute the time step dtion
190 : :
191 [ + - ]: 67 : if (ab_mover%dtion>0)then
192 : 67 : hh = ab_mover%dtion
193 : : else
194 : 0 : hh=fdtion(ab_mover,itime,xcart,fcart_cur,vel_cur)
195 : : end if
196 : :
197 : : !write(std_out,*) '02'
198 : : !##########################################################
199 : : !### 02. For all atoms and directions
200 [ + + ]: 366 : do kk=1,ab_mover%natom
201 : 299 : em=ab_mover%amass(kk)
202 [ + + ]: 1263 : do jj=1,3
203 : :
204 : : ! write(std_out,*) '03'
205 : : ! ##########################################################
206 : : ! ### 03. Filling other values from history (forces and vel)
207 : 897 : fcart=fcart_cur(jj,kk)
208 : 897 : xc=xcart(jj,kk)
209 : 897 : vel=hist%vel(jj,kk,1)
210 : :
211 : : ! Previous values only after first iteration
212 [ + + + + ]: 897 : if (itime>=2.or.icycle>=2) then
213 : 813 : fprev=fcart_prev(jj,kk)
214 : 813 : vprev=hist%vel(jj,kk,hist%ihist)
215 : : end if
216 [ + + + + ]: 897 : if (itime>=3.or.icycle>=3) then
217 : 645 : fprev2=fcart_prev2(jj,kk)
218 : : end if
219 : :
220 [ + + ]: 897 : if (itime==2)then
221 : 84 : vec_tmp1(jj,kk)=xcart_prev(jj,kk)
222 : 84 : vec_tmp2(jj,kk)=xcart(jj,kk)
223 : : end if
224 : :
225 : : ! write(std_out,*) '04'
226 : : ! ##########################################################
227 : : ! ### 04. Take first the atoms that are not allowed to move along
228 : : ! ### this direction
229 : : ! ### Warning : implemented in cartesian coordinates
230 [ + + ]: 897 : if (ab_mover%iatfix(jj,kk)==1) then
231 : : ! Their positions will be the same as xcart
232 : 445 : xnow=xcart(jj,kk)
233 : : ! Their velocities are zero
234 : 445 : vnow=0.0_dp
235 : : else
236 : :
237 : : ! write(std_out,*) '05'
238 : : ! ##########################################################
239 : : ! ### 05. Initialization (itime==1):
240 : : ! ### 4 calls to obtain the forces are neeeded
241 : : ! ### The variables vec_tmp2 and vec_tmp1 from previous
242 : : ! ### calls are used in the following ones.
243 [ + + ]: 452 : if(itime==1)then
244 : 152 : x0=xcart_prev(jj,kk)
245 : :
246 : : ! Prepare the second cycle
247 [ + + ]: 152 : if(icycle==1)then
248 : 38 : dx=hh*vel
249 : 38 : dv=hh/em*(fcart-vis*vel)
250 : 38 : xnow=x0+.5_dp*dx
251 : 38 : vnow=vel+.5_dp*dv
252 : 38 : vec_tmp2(jj,kk)=xc+sixth*dx
253 : 38 : vec_tmp1(jj,kk)=vel+sixth*dv
254 [ + + ]: 114 : else if(icycle==2)then
255 : 38 : dx=hh*vprev
256 : 38 : dv=hh/em*(fcart-vis*vprev)
257 : 38 : xnow=x0+.5_dp*dx
258 : 38 : vnow=vel+.5_dp*dv
259 : 38 : vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx
260 : 38 : vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv
261 [ + + ]: 76 : else if(icycle==3)then
262 : 38 : dx=hh*vprev
263 : 38 : dv=hh/em*(fcart-vis*vprev)
264 : 38 : xnow=x0+dx
265 : 38 : vnow=vel+dv
266 : 38 : vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx
267 : 38 : vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv
268 [ + - ]: 38 : else if(icycle==4)then
269 : 38 : dx=hh*vprev
270 : 38 : dv=hh/em*(fcart-vis*vprev)
271 : 38 : xnow=vec_tmp2(jj,kk)+sixth*dx
272 : 38 : vnow=vec_tmp1(jj,kk)+sixth*dv
273 : : end if
274 : : else !(itime/=1)
275 : :
276 : : ! write(std_out,*) '06'
277 : : ! ##########################################################
278 : : ! ### 06. Change positions and velocities
279 : : ! ### These changes only applies for itime>2
280 [ + + ]: 300 : if (itime>2)then
281 : : ! Uses a corrector to have better value of xnow, and
282 : : ! derive vnow. Only update atoms position and
283 : : ! velocity along its allowed directions
284 : 262 : aa=fprev
285 : 262 : bb=(fcart-fprev2)/(2._dp*hh)
286 : 262 : cc=(fcart+fprev2-2._dp*fprev)/(2._dp*hh*hh)
287 : 262 : x0=vec_tmp2(jj,kk)
288 : 262 : xm=vec_tmp1(jj,kk)
289 [ - + ]: 262 : if(abs(vis)<=1.d-8)then
290 : : ! NON-DAMPED DYNAMICS (Post-Code)
291 : : xnow=2._dp*x0-xm+hh**2/em/12._dp*&
292 : 0 : & (fprev2+10._dp*fprev+fcart)
293 : : vnow=(bb*hh**2)/(3._dp*em)&
294 : : & +1.5_dp*aa*hh/em+&
295 : : & (5._dp/12._dp)*cc*hh**3/em&
296 : 0 : & +x0/hh-xm/hh
297 : : else
298 : : ! DAMPED DYNAMICS (Post-Code)
299 : 262 : alfa=exp(-vis*hh/em)
300 : : xnow=((-aa*hh*vis**2+0.5_dp*bb*hh**2*vis**2&
301 : : & -third*cc*hh**3*vis**2+em*bb*hh*vis&
302 : : & -em*cc*hh**2*vis-2._dp*em**2*cc*hh+x0*vis**3-xm*vis**3)*alfa&
303 : : & +aa*hh*vis**2-em*bb*hh*vis+third*cc*hh**3*vis**2&
304 : : & +2._dp*em**2*cc*hh+0.5D0*bb*hh**2*vis**2-em*cc*hh**2*vis+x0*vis**3)&
305 : 262 : & /vis**3
306 : : vnow=(em*aa*vis**2*alfa-em*aa*vis**2+bb*hh*vis**2*em*alfa&
307 : : & -bb*hh*vis**2*em+cc*hh**2*vis**2*em*alfa-cc*hh**2*vis**2*em&
308 : : & -em**2*bb*vis*alfa+em**2*bb*vis-2._dp*em**2*cc*hh*vis*alfa+&
309 : : & 2._dp*em**2*cc*hh*vis+2._dp*em**3*cc*alfa-2._dp*em**3*cc+&
310 : : & vis**3*alfa**2*aa*hh-0.5_dp*vis**3*alfa**2*bb*hh**2+&
311 : : & third*vis**3*alfa**2*cc*hh**3-vis**2*&
312 : : & alfa**2*em*bb*hh+vis**2*alfa**2*em*cc*hh**2+&
313 : : & 2._dp*vis*alfa**2*em**2*cc*hh-vis**4*alfa**2*x0+&
314 : 262 : & vis**4*alfa**2*xm)/vis**3/(alfa-1._dp)/em
315 : :
316 : : end if !if(abs(vis)<=1.d-8)
317 : :
318 : 262 : xc=xnow
319 : 262 : vec_tmp1(jj,kk)=vec_tmp2(jj,kk)
320 : 262 : vec_tmp2(jj,kk)=xnow
321 : : else
322 : : vnow=vprev
323 : : end if !if(itime>2)
324 : :
325 : : ! write(std_out,*) '07'
326 : : ! ##########################################################
327 : : ! ### 07. Correct positions
328 : : ! ### These changes only applies for itime>1
329 : :
330 [ - + ]: 300 : if(abs(vis)<=1.d-8)then
331 : : ! NON-DAMPED DYNAMICS (Pre-Code)
332 : : ! If the viscosity is too small, the equations become
333 : : ! ill conditioned due to rounding error so do regular
334 : : ! Verlet predictor Numerov corrector.
335 : 0 : x0=vec_tmp2(jj,kk)
336 : 0 : xm=vec_tmp1(jj,kk)
337 : : xnow=2._dp*x0-xm&
338 : 0 : & + hh**2/em*fcart
339 : : else
340 : : ! DAMPED DYNAMICS (Pre-Code)
341 : : ! These equations come from solving
342 : : ! m*d2x/dt2+vis*dx/dt=a+b*t+c*t**2
343 : : ! analytically under the boundary conditions that
344 : : ! x(0)=x0 and x(-h)=xm, and the following is the
345 : : ! expression for x(h). a, b and c are determined
346 : : ! from our knowledge of the driving forces.
347 : 300 : aa=fcart
348 : 300 : bb=(fcart-fprev)/hh
349 : 300 : x0=vec_tmp2(jj,kk)
350 : 300 : xm=vec_tmp1(jj,kk)
351 : 300 : alfa=exp(-vis*hh/em)
352 : : xnow=( (-aa*hh*vis**2 +0.5_dp*bb*hh**2*vis**2&
353 : : & +em*bb*hh*vis +x0*vis**3 -xm*vis**3)*alfa&
354 : : & +aa*hh*vis**2 -em*bb*hh*vis&
355 : 300 : & +0.5_dp*bb*hh**2*vis**2 +x0*vis**3)/vis**3
356 : : ! End of choice between initialisation, damped
357 : : ! dynamics and non-damped dynamics
358 : : end if
359 : :
360 : : end if !if(itime==1)
361 : :
362 : : end if !if(ab_mover%iatfix(jj,kk)==1)
363 : :
364 : : ! write(std_out,*) '08'
365 : : ! ##########################################################
366 : : ! ### 08. Update history
367 : :
368 : 897 : xcart(jj,kk)=xnow
369 : 1196 : vel_next(jj,kk)=vnow
370 : :
371 : : ! write(std_out,*) '09'
372 : : ! ##########################################################
373 : : ! ### 09. End loops of atoms and directions
374 : : end do ! jj=1,3
375 : : end do ! kk=1,ab_mover%natom
376 : :
377 : : !write(std_out,*) '10'
378 : : !##########################################################
379 : : !### 10. Filling history with the new values
380 : :
381 : 67 : hist%ihist = abihist_findIndex(hist,+1)
382 : :
383 : 67 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
384 : 67 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
385 : :
386 : : !Change ncycle for itime>1
387 [ + + ]: 67 : if (icycle==4) ncycle=1
388 : :
389 [ + + ]: 67 : if (itime==1)then
390 : 24 : time=0.0_dp
391 [ - + ]: 24 : if (ab_mover%dtion<0)then
392 : 0 : write(std_out,*) 'Time=',time
393 : : end if
394 : : end if
395 : 67 : time=time+hh
396 : 67 : hist%time(hist%ihist)=time
397 : :
398 [ + - ]: 67 : ABI_SFREE(xcart)
399 [ + + ]: 67 : ABI_SFREE(xcart_prev)
400 : :
401 : : if (.false.) write(std_out,*) ntime
402 : :
403 [ + + - + ]: 140 : end subroutine pred_moldyn
404 : : !!***
405 : :
406 : : end module m_pred_moldyn
407 : : !!***
|
