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), ABI_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 : !!***
|
