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