Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_verlet
2 : : !! NAME
3 : : !! m_pred_verlet
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JCC, SE)
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_verlet
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_abimover
27 : : use m_abihist
28 : : use m_xfpack
29 : :
30 : : use m_geometry, only : mkrdim, xcart2xred, xred2xcart, fcart2gred, metric
31 : :
32 : : implicit none
33 : :
34 : : private
35 : : !!***
36 : :
37 : : public :: pred_verlet
38 : : !!***
39 : :
40 : : contains
41 : : !!***
42 : :
43 : : !!****f* ABINIT/pred_verlet
44 : : !! NAME
45 : : !! pred_verlet
46 : : !!
47 : : !! FUNCTION
48 : : !! Ionmov predictors (6 & 7) Verlet algorithm
49 : : !!
50 : : !! IONMOV 6:
51 : : !! Given a starting point xred that is a vector of length 3*natom
52 : : !! (reduced nuclei coordinates), a velocity vector (in cartesian
53 : : !! coordinates), and unit cell parameters (acell and rprimd -
54 : : !! without velocities in the present implementation), the Verlet
55 : : !! dynamics is performed, using the gradient of the energy
56 : : !! (atomic forces and stresses) as calculated by the routine scfcv.
57 : : !!
58 : : !! Some atoms can be kept fixed, while the propagation of unit cell
59 : : !! parameters is only performed if optcell/=0.
60 : : !! No more than ab_mover%ntime steps are performed.
61 : : !! The time step is governed by dtion, contained in ab_mover
62 : : !! (coming from dtset).
63 : : !! Returned quantities are xred, and eventually acell and rprimd (new ones!).
64 : : !!
65 : : !! IONMOV 7:
66 : : !! Block every atom for which the scalar product of velocity and
67 : : !! forces is negative, in order to reach the minimum.
68 : : !! The convergence requirement on the atomic forces, ab_mover%tolmxf,
69 : : !! allows an early exit.
70 : : !!
71 : : !! INPUTS
72 : : !! ab_mover <type(abimover)> : Datatype with all the information needed by the preditor
73 : : !! itime : Index of the present iteration
74 : : !! ntime : Maximal number of iterations
75 : : !! ionmov : (6 or 7) Specific kind of VERLET
76 : : !! zDEBUG : if true print some debugging information
77 : : !!
78 : : !! SIDE EFFECTS
79 : : !! hist <type(abihist)> : History of positions,forces
80 : : !! acell, rprimd, stresses
81 : : !!
82 : : !! SOURCE
83 : :
84 : 53 : subroutine pred_verlet(ab_mover,hist,ionmov,itime,ntime,zDEBUG,iexit)
85 : :
86 : : implicit none
87 : :
88 : : !Arguments ------------------------------------
89 : : !scalars
90 : : type(abimover),intent(in) :: ab_mover
91 : : type(abihist),intent(inout) :: hist
92 : : integer,intent(in) :: itime
93 : : integer,intent(in) :: ntime
94 : : integer,intent(in) :: ionmov
95 : : integer,intent(in) :: iexit
96 : : logical,intent(in) :: zDEBUG
97 : :
98 : : !Local variables-------------------------------
99 : : !scalars
100 : : integer :: ii,istopped,jj,kk,ndim,nstopped
101 : : real(dp) :: amass_tot,etotal,diag,gr_avg,ekin_corr,scprod,taylor,ucvol0
102 : : real(dp),save :: ucvol,ucvol_next
103 : : character(len=500) :: message
104 : : !arrays
105 : 106 : integer :: stopped(ab_mover%natom)
106 : 106 : real(dp) :: acell0(3),fcart(3,ab_mover%natom)
107 : 106 : real(dp) :: gred_corrected(3,ab_mover%natom)
108 : : real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3), strten(6)
109 : 106 : real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
110 : 106 : real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
111 : 53 : real(dp) :: vel(3,ab_mover%natom),vel_nexthalf(3,ab_mover%natom)
112 : : real(dp),save :: acell(3),acell_next(3)
113 : : real(dp),save :: rprimd(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3)
114 : : real(dp),allocatable,save :: hessin(:,:)
115 : : real(dp),allocatable,save :: vin(:),vin_prev(:),vin_next(:)
116 : : real(dp),allocatable,save :: vout(:),vout_prev(:),vel_prevhalf(:,:)
117 : :
118 : : !***************************************************************************
119 : : !Beginning of executable session
120 : : !***************************************************************************
121 : :
122 [ + + ]: 53 : if(iexit/=0)then
123 [ + - ]: 5 : ABI_SFREE(vin)
124 [ + - ]: 5 : ABI_SFREE(vin_next)
125 [ + - ]: 5 : ABI_SFREE(vout)
126 [ + - ]: 5 : ABI_SFREE(vin_prev)
127 [ + - ]: 5 : ABI_SFREE(vout_prev)
128 [ + - ]: 5 : ABI_SFREE(hessin)
129 [ + - ]: 5 : ABI_SFREE(vel_prevhalf)
130 : : return
131 : : end if
132 : :
133 : : !write(std_out,*) 'verlet 01'
134 : : !##########################################################
135 : : !### 01. Compute the dimension of vectors (ndim)
136 : :
137 : 48 : ndim=3*ab_mover%natom
138 [ - + ]: 48 : if(ab_mover%optcell==1) ndim=ndim+1
139 [ - + ]: 48 : if(ab_mover%optcell==2 .or.&
140 : 0 : & ab_mover%optcell==3) ndim=ndim+6
141 [ - + ]: 48 : if(ab_mover%optcell>=4) ndim=ndim+3
142 : :
143 : : !write(std_out,*) 'verlet 02'
144 : : !##########################################################
145 : : !### 02. Allocate the vectors vin, vout and hessian matrix
146 : :
147 : : !Notice that vin, vout, etc could be allocated
148 : : !From a previous dataset with a different ndim
149 [ + + ]: 48 : if(itime==1)then
150 [ - + ]: 5 : ABI_SFREE(vin)
151 [ - + ]: 5 : ABI_SFREE(vin_next)
152 [ - + ]: 5 : ABI_SFREE(vout)
153 [ - + ]: 5 : ABI_SFREE(vin_prev)
154 [ - + ]: 5 : ABI_SFREE(vout_prev)
155 [ - + ]: 5 : ABI_SFREE(hessin)
156 [ - + ]: 5 : ABI_SFREE(vel_prevhalf)
157 : :
158 [ - + + - : 15 : ABI_MALLOC(vin,(ndim))
+ - + - -
+ - + -
+ ]
159 [ - + + - : 10 : ABI_MALLOC(vin_next,(ndim))
- + - + -
+ ]
160 [ - + + - : 10 : ABI_MALLOC(vout,(ndim))
- + - + -
+ ]
161 [ - + + - : 10 : ABI_MALLOC(vin_prev,(ndim))
- + - + -
+ ]
162 [ - + + - : 10 : ABI_MALLOC(vout_prev,(ndim))
- + - + -
+ ]
163 [ - + + - : 20 : ABI_MALLOC(hessin,(ndim,ndim))
- + + - +
- + - - +
- + - + ]
164 [ - + + - : 15 : ABI_MALLOC(vel_prevhalf,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
165 : : end if
166 : :
167 : : !write(std_out,*) 'verlet 03'
168 : : !##########################################################
169 : : !### 03. Obtain the present values from the history
170 : :
171 : 48 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
172 : :
173 : 48 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
174 [ + + ]: 192 : do ii=1,3
175 [ + + ]: 624 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
176 : : end do
177 : :
178 : 48 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
179 [ + + + + ]: 432 : fcart(:,:) =hist%fcart(:,:,hist%ihist)
180 [ + + ]: 336 : strten(:) =hist%strten(:,hist%ihist)
181 [ + + + + ]: 432 : vel(:,:) =hist%vel(:,:,hist%ihist)
182 : 48 : etotal =hist%etot(hist%ihist)
183 : :
184 [ - + ]: 48 : if(zDEBUG)then
185 : 0 : write (std_out,*) 'fcart:'
186 [ # # ]: 0 : do kk=1,ab_mover%natom
187 : 0 : write (std_out,*) fcart(:,kk)
188 : : end do
189 : 0 : write (std_out,*) 'vel:'
190 [ # # ]: 0 : do kk=1,ab_mover%natom
191 : 0 : write (std_out,*) vel(:,kk)
192 : : end do
193 : 0 : write (std_out,*) 'strten:'
194 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
195 : 0 : write (std_out,*) 'etotal:'
196 : 0 : write (std_out,*) etotal
197 : : end if
198 : :
199 : 48 : acell0=acell ; ucvol0=ucvol
200 : :
201 : : !Get rid of mean force on whole unit cell, but only if no
202 : : !generalized constraints are in effect
203 : 48 : call fcart2gred(fcart,gred_corrected,rprimd,ab_mover%natom)
204 [ + - ]: 48 : if(ab_mover%nconeq==0)then
205 [ + + ]: 144 : amass_tot=sum(ab_mover%amass(:))
206 [ + + ]: 192 : do kk=1,3
207 [ + + + - ]: 192 : if (kk/=3.or.ab_mover%jellslab==0) then
208 [ + + ]: 432 : gr_avg=sum(gred_corrected(kk,:))/dble(ab_mover%natom)
209 [ + + ]: 432 : gred_corrected(kk,:)=gred_corrected(kk,:)-gr_avg*ab_mover%amass(:)/amass_tot
210 : : end if
211 : : end do
212 : : end if
213 : :
214 : : !write(std_out,*) 'verlet 04'
215 : : !##########################################################
216 : : !### 04. Fill the vectors vin and vout
217 : :
218 : : !Initialize input vectors : first vin, then vout
219 : : call xfpack_x2vin(acell, ab_mover%natom, ndim,&
220 : : & ab_mover%nsym, ab_mover%optcell, rprim, rprimd,&
221 [ + - - - : 48 : & ab_mover%symrel, ucvol, ucvol0, vin, xred)
- - - - -
- - + ]
222 : : call xfpack_f2vout(gred_corrected, ab_mover%natom, ndim,&
223 : : & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
224 [ + - - - : 48 : & vout)
- - - + ]
225 : :
226 : : !write(std_out,*) 'verlet 05'
227 : : !##########################################################
228 : : !### 05. Initialize or update the hessian matrix
229 : :
230 : : !Here, set up the matrix of transformation between forces and
231 : : !acceleration. Masses must be included here.
232 : : !Beside this feature, one could define
233 : : !a preconditioner, in which case it should
234 : : !be the inverse hessian, like in Broyden. This explains the
235 : : !name chosen for this transformation matrix. This would allow
236 : : !to find easily the optimal geometry with ionmov=7.
237 : : !The default, now implemented, corresponds to the identity matrix
238 : : !in cartesian coordinates, which makes use of metric tensor gmet
239 : : !in reduced coordinates.
240 : :
241 : : !Initialise the Hessian matrix using gmet
242 [ + + ]: 48 : if (itime==1)then
243 : : ! Initialize inverse hessian with identity matrix
244 : : ! in cartesian coordinates, which makes use of metric tensor gmet
245 : : ! in reduced coordinates.
246 [ + + + + ]: 215 : hessin(:,:)=zero
247 [ + + ]: 15 : do ii=1,ab_mover%natom
248 [ + + ]: 45 : do kk=1,3
249 [ + + ]: 130 : do jj=1,3
250 : : ! Warning : implemented in reduced coordinates
251 [ + - + - ]: 90 : if (ab_mover%iatfix(kk,ii)==0 .and.&
252 : 30 : & ab_mover%iatfix(jj,ii)==0 )then
253 : 90 : hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)/ab_mover%amass(ii)
254 : : end if
255 : : end do
256 : : end do
257 : : end do
258 [ - + ]: 5 : if(ab_mover%optcell/=0)then
259 : : ! These values might lead to too large changes in some cases ...
260 : 0 : diag=ab_mover%strprecon*30.0_dp/ucvol
261 [ # # ]: 0 : if(ab_mover%optcell==1) diag=diag/three
262 [ # # ]: 0 : do ii=3*ab_mover%natom+1,ndim
263 : 0 : hessin(ii,ii)=diag
264 : : end do
265 : : end if
266 : : end if
267 : :
268 : : !zDEBUG (vin,vout and hessin before prediction)
269 [ - + ]: 48 : if(zDEBUG)then
270 : 0 : write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
271 : 0 : write(std_out,*) 'vin:'
272 [ # # # # ]: 0 : do ii=1,ndim,3
273 [ # # ]: 0 : if (ii+2<=ndim)then
274 : 0 : write(std_out,*) ii,vin(ii:ii+2)
275 : : else
276 : 0 : write(std_out,*) ii,vin(ii:ndim)
277 : : end if
278 : : end do
279 : 0 : write(std_out,*) 'vout:'
280 [ # # # # ]: 0 : do ii=1,ndim,3
281 [ # # ]: 0 : if (ii+2<=ndim)then
282 : 0 : write(std_out,*) ii,vout(ii:ii+2)
283 : : else
284 : 0 : write(std_out,*) ii,vout(ii:ndim)
285 : : end if
286 : : end do
287 : 0 : write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
288 [ # # ]: 0 : do kk=1,ndim
289 [ # # # # ]: 0 : do jj=1,ndim,3
290 [ # # ]: 0 : if (jj+2<=ndim)then
291 : 0 : write(std_out,*) jj,hessin(jj:jj+2,kk)
292 : : else
293 : 0 : write(std_out,*) jj,hessin(jj:ndim,kk)
294 : : end if
295 : : end do
296 : : end do
297 : : end if
298 : :
299 : : !write(std_out,*) 'verlet 06'
300 : : !##########################################################
301 : : !### 06. Compute the next values
302 : :
303 : : !%%% VERLET ALGORITHM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
304 : :
305 [ - + ]: 48 : if(zDEBUG)then
306 : 0 : write(std_out,*) 'Shifted data ENTER'
307 : 0 : write(std_out,*) 'acell(:)',acell(:)
308 : 0 : write(std_out,*) 'rprimd(:,:)',rprimd(:,:)
309 : 0 : write(std_out,*) 'ucvol',ucvol
310 : 0 : write(std_out,*) 'vel_prevhalf(:,:)',vel_prevhalf(:,:)
311 : 0 : write(std_out,*) 'vin_prev(:)',vin_prev(:)
312 : 0 : write(std_out,*) 'vin(:)',vin(:)
313 : 0 : write(std_out,*) 'xcart(:,:)',xcart(:,:)
314 : 0 : write(std_out,*) 'xred(:,:)',xred(:,:)
315 : : end if
316 : :
317 : : !Compute next atomic coordinates and cell parameters, using
318 : : !Verlet algorithm
319 : :
320 : : !1. First propagate the position, without acceleration
321 [ + + ]: 48 : if(itime/=1)then
322 [ + + ]: 301 : vin_next(:)=2*vin(:)-vin_prev(:)
323 : : taylor=one
324 : : else
325 : : ! Initialisation : no vin_prev is available, but the ionic velocity
326 : : ! is available, in cartesian coordinates
327 : : ! Uses the velocity
328 [ + + + + ]: 45 : xcart_next(:,:)=xcart(:,:)+ab_mover%dtion*vel(:,:)
329 : : ! Convert back to xred_next (reduced coordinates)
330 : 5 : call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
331 : : ! Impose no change of acell, ucvol, rprim, and rprimd
332 : 5 : acell_next(:)=acell(:)
333 : 5 : ucvol_next=ucvol
334 : 5 : rprim_next(:,:)=rprim(:,:)
335 : 5 : rprimd_next(:,:)=rprimd(:,:)
336 : 5 : write(std_out,*) 'ucvol',ucvol
337 : : ! Store all these next values in vin_next
338 : : call xfpack_x2vin(acell_next,ab_mover%natom,&
339 : : & ndim,ab_mover%nsym,ab_mover%optcell,rprim_next,&
340 : : & rprimd,ab_mover%symrel,ucvol_next,ucvol0,&
341 [ + - - - : 5 : & vin_next,xred_next)
- - - - -
- - + ]
342 : 5 : taylor=half
343 : : end if
344 : :
345 : : !2. Now, take into account the acceleration
346 [ + + ]: 336 : do ii=1,ndim
347 : : ! Note the minus sign: the forces are minus the gradients,
348 : : ! contained in vout.
349 : : vin_next(:)=vin_next(:)-ab_mover%dtion**2*hessin(:,ii)*&
350 [ + + ]: 2064 : & vout(ii)*taylor
351 : : end do
352 : :
353 : : !3. Implement fixing of atoms : put back old values for fixed
354 : : !components
355 [ + + ]: 144 : do kk=1,ab_mover%natom
356 [ + + ]: 432 : do jj=1,3
357 : : ! Warning : implemented in reduced coordinates
358 [ - + ]: 384 : if (ab_mover%iatfix(jj,kk) == 1) then
359 : 0 : vin_next(jj+(kk-1)*3)=vin(jj+(kk-1)*3)
360 : : end if
361 : : end do
362 : : end do
363 : :
364 : : !4. Now, compute the velocity at the next half-step
365 : : !Get xred_next, and eventually acell_next, ucvol_next, rprim_next and
366 : : !rprimd_next, from vin_next
367 : : call xfpack_vin2x(acell_next,acell0,ab_mover%natom,ndim,&
368 : : & ab_mover%nsym,ab_mover%optcell,rprim_next,rprimd,&
369 [ + - - - : 48 : & ab_mover%symrel,ucvol_next,ucvol0,vin_next,xred_next)
- - - - -
- - + ]
370 [ - + ]: 48 : if(ab_mover%optcell/=0)then
371 : 0 : call mkrdim(acell_next,rprim_next,rprimd_next)
372 : 0 : call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next)
373 : : else
374 : 48 : rprimd_next(:,:)=rprimd(:,:)
375 : : end if
376 : : !Convert input xred_next (reduced coordinates) to
377 : : !xcart_next (cartesian)
378 : 48 : call xred2xcart(ab_mover%natom,rprimd_next,xcart_next,xred_next)
379 : : !Compute the velocity at half of the new step
380 [ + + + + ]: 432 : vel_nexthalf(:,:)=(xcart_next(:,:)-xcart(:,:))/ab_mover%dtion
381 : :
382 : : !5. If needed, compute the velocity at present position
383 [ + + ]: 48 : if(itime/=1)then
384 [ + + + + ]: 387 : vel(:,:)=(vel_nexthalf(:,:)+vel_prevhalf(:,:))*0.5_dp
385 : : end if
386 : :
387 : : !%%% VERLET ALGORITHM BLOCKING ATOMS %%%%%%%%%%%%%%%%%%%%%%
388 : :
389 : : !Here, stop the atoms for which the scalar product of velocity
390 : : !and force is negative, and recompute the kinetic energy.
391 [ + + ]: 48 : if(ionmov==7)then
392 [ + + ]: 87 : stopped(:)=0
393 [ + + ]: 87 : do ii=1,ab_mover%natom
394 : : scprod=fcart(1,ii)*vel(1,ii)+&
395 : : & fcart(2,ii)*vel(2,ii)+&
396 : 58 : & fcart(3,ii)*vel(3,ii)
397 [ + + + - ]: 87 : if(scprod<0.0_dp .and. itime/=1)then
398 : 26 : stopped(ii)=1
399 : 26 : write(std_out,*) 'Stopped atom',ii
400 : : ! Shift the velocities of the previous half-step and current
401 : : ! half-step, so that the acceleration is correct but the
402 : : ! present velocity vanishes.
403 [ + + ]: 104 : vel_prevhalf(:,ii)=vel_prevhalf(:,ii)-vel(:,ii)
404 [ + + ]: 104 : vel_nexthalf(:,ii)=vel_nexthalf(:,ii)-vel(:,ii)
405 [ + + ]: 104 : vel(:,ii)=0.0_dp
406 [ + + ]: 104 : xcart_next(:,ii)=xcart(:,ii)+ab_mover%dtion*vel_nexthalf(:,ii)
407 : : end if
408 : : end do
409 : :
410 [ - + ]: 29 : if(zDEBUG)then
411 : 0 : write (std_out,*) 'fcart:'
412 [ # # ]: 0 : do kk=1,ab_mover%natom
413 : 0 : write (std_out,*) fcart(:,kk)
414 : : end do
415 : 0 : write (std_out,*) 'vel_prevhalf:'
416 [ # # ]: 0 : do kk=1,ab_mover%natom
417 : 0 : write (std_out,*) vel_prevhalf(:,kk)
418 : : end do
419 : 0 : write (std_out,*) 'vel_nexthalf:'
420 [ # # ]: 0 : do kk=1,ab_mover%natom
421 : 0 : write (std_out,*) vel_nexthalf(:,kk)
422 : : end do
423 : 0 : write (std_out,*) 'vel:'
424 [ # # ]: 0 : do kk=1,ab_mover%natom
425 : 0 : write (std_out,*) vel(:,kk)
426 : : end do
427 : 0 : write (std_out,*) 'xcart_next:'
428 [ # # ]: 0 : do kk=1,ab_mover%natom
429 : 0 : write (std_out,*) xcart_next(:,kk)
430 : : end do
431 : : end if
432 : :
433 : : ! Establish a list of stopped atoms
434 [ + + ]: 87 : nstopped=sum(stopped(:))
435 : :
436 [ + + ]: 29 : if(nstopped/=0)then
437 : 13 : write(message,'(a)') ' List of stopped atoms (ionmov=7) :'
438 : 13 : call wrtout(ab_out,message,'COLL')
439 : 13 : istopped=1
440 [ + + ]: 39 : do ii=1,ab_mover%natom
441 [ + - ]: 39 : if(stopped(ii)==1)then
442 : 26 : stopped(istopped)=ii
443 : 26 : istopped=istopped+1
444 : : end if
445 : : end do
446 [ - + - + ]: 26 : do ii=1,nstopped,16
447 : 13 : write(message, '(16i4)' ) stopped(ii:min(ii+15,nstopped))
448 : 26 : call wrtout(ab_out,message,'COLL')
449 : : end do
450 : : ! Now, compute the corrected kinetic energy
451 : : ! Generate xred_next from xcart_next
452 : 13 : call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next)
453 : : ! Store xred_next, and eventual acell_next and rprim_next in vin
454 : : call xfpack_x2vin(acell_next,&
455 : : & ab_mover%natom,ndim,ab_mover%nsym,ab_mover%optcell,&
456 : : & rprim_next,rprimd,&
457 [ + - - - : 13 : & ab_mover%symrel,ucvol_next,ucvol0,vin_next,xred_next)
- - - - -
- - + ]
458 : :
459 : 13 : ekin_corr=0.0_dp
460 [ + + ]: 39 : do ii=1,ab_mover%natom
461 [ + + ]: 117 : do jj=1,3
462 : : ! Warning : the fixing of atomis is implemented in reduced
463 : : ! coordinates, so that this expression is wrong
464 : 104 : if (ab_mover%iatfix(jj,ii) == 0) then
465 : : ekin_corr=ekin_corr+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2
466 : : end if
467 : : end do
468 : : end do
469 : : ! End of test nstopped/=0
470 : : end if
471 : :
472 : : ! End of test ionmov==7
473 : : end if
474 : :
475 : : !write(std_out,*) 'verlet 07'
476 : : !##########################################################
477 : : !### 07. Shift the data from next values to the present
478 : :
479 : : !acell(:)=acell_next(:)
480 : : !rprim(:,:)=rprim_next(:,:)
481 : : !ucvol=ucvol_next
482 [ + + + + ]: 432 : vel_prevhalf(:,:)=vel_nexthalf(:,:)
483 [ + + ]: 336 : vin_prev(:)=vin(:)
484 [ + + ]: 336 : vin(:)=vin_next(:)
485 [ + + + + ]: 432 : xcart(:,:)=xcart_next(:,:)
486 [ + + + + ]: 432 : xred(:,:)=xred_next(:,:)
487 : :
488 : 48 : write(std_out,*) 'Shifted data EXIT'
489 : 48 : write(std_out,*) 'acell(:)',acell(:)
490 : 48 : write(std_out,*) 'rprim(:,:)',rprim(:,:)
491 : 48 : write(std_out,*) 'rprimd(:,:)',rprimd(:,:)
492 : 48 : write(std_out,*) 'ucvol',ucvol
493 : 48 : write(std_out,*) 'vel_prevhalf(:,:)',vel_prevhalf(:,:)
494 : 48 : write(std_out,*) 'vin_prev(:)',vin_prev(:)
495 : 48 : write(std_out,*) 'vin(:)',vin(:)
496 : 48 : write(std_out,*) 'xred(:,:)',xred(:,:)
497 : :
498 : :
499 : : !write(std_out,*) 'verlet 08'
500 : : !##########################################################
501 : : !### 08. Update the history with the prediction
502 : :
503 : : !Increase indexes
504 : 48 : hist%ihist=abihist_findIndex(hist,+1)
505 : :
506 : : !Fill the history with the variables
507 : : !xred, acell, rprimd, vel
508 : 48 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
509 [ + + + + ]: 432 : hist%vel(:,:,hist%ihist)=vel(:,:)
510 : 48 : hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
511 : :
512 [ - + ]: 48 : if(zDEBUG)then
513 : 0 : write (std_out,*) 'vel:'
514 [ # # ]: 0 : do kk=1,ab_mover%natom
515 : 0 : write (std_out,*) vel(:,kk)
516 : : end do
517 : : end if
518 : :
519 : : if (.false.) write(std_out,*) ntime
520 : :
521 : : end subroutine pred_verlet
522 : : !!***
523 : :
524 : : end module m_pred_verlet
525 : : !!***
|
