Line data Source code
1 : !!****m* ABINIT/m_pred_isokinetic
2 : !! NAME
3 : !! m_pred_isokinetic
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_isokinetic
23 :
24 : use m_abicore
25 : use defs_basis
26 : use m_abimover
27 : use m_abihist
28 :
29 : use m_numeric_tools, only : uniformrandom
30 : use m_geometry, only : xcart2xred, xred2xcart
31 :
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: pred_isokinetic
39 : !!***
40 :
41 : contains
42 : !!***
43 :
44 :
45 : !!****f* ABINIT/pred_isokinetic
46 : !! NAME
47 : !! pred_isokinetic
48 : !!
49 : !! FUNCTION
50 : !! Ionmov predictors (12) Isokinetic ensemble molecular dynamics
51 : !!
52 : !! IONMOV 12:
53 : !! Isokinetic ensemble molecular dynamics.
54 : !! The equation of motion of the ions in contact with a thermostat
55 : !! are solved with the algorithm proposed by Zhang [J. Chem. Phys. 106, 6102 (1997)] [[cite:Zhang1997]],
56 : !! as worked out by Minary et al, J. Chem. Phys. 188, 2510 (2003) [[cite:Minary2003]].
57 : !! The conservation of the kinetic energy is obtained within machine precision, at each step.
58 : !! Related parameters: the time step (dtion), the initial temperature (mdtemp(1)) if the velocities are not defined to start with.
59 : !!
60 : !! INPUTS
61 : !! ab_mover <type(abimover)> : Datatype with all the information needed by the preditor
62 : !! itime : Index of the present iteration
63 : !! ntime : Maximal number of iterations
64 : !! zDEBUG : if true print some debugging information
65 : !!
66 : !! OUTPUT
67 : !!
68 : !! SIDE EFFECTS
69 : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
70 : !! SOURCE
71 :
72 978 : subroutine pred_isokinetic(ab_mover,hist,itime,ntime,zDEBUG,iexit)
73 :
74 : !Arguments ------------------------------------
75 : !scalars
76 : integer,intent(in) :: itime
77 : integer,intent(in) :: ntime
78 : integer,intent(in) :: iexit
79 : logical,intent(in) :: zDEBUG
80 : type(abimover),intent(in) :: ab_mover
81 : type(abihist),intent(inout) :: hist
82 :
83 : !Local variables-------------------------------
84 : !scalars
85 : integer :: kk,iatom,idim,idum=5,nxyzatfree,ndegfreedom,nfirst,ifirst
86 : real(dp) :: a,as,b,sqb,s,s1,s2,scdot,sigma2,vtest,v2gauss
87 : real(dp),parameter :: v2tol=tol8
88 : real(dp) :: etotal,rescale_vel
89 : character(len=5000) :: message
90 : !arrays
91 : real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:)
92 :
93 : real(dp) :: acell(3),rprimd(3,3)
94 1956 : real(dp) :: fcart(3,ab_mover%natom)
95 1956 : real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
96 1956 : real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
97 978 : real(dp) :: vel(3,ab_mover%natom)
98 : real(dp) :: strten(6)
99 :
100 : !***************************************************************************
101 : !Beginning of executable session
102 : !***************************************************************************
103 :
104 : !DEBUG
105 : !write(std_out,*)' pred_isokinetic : enter '
106 : !stop
107 : !ENDDEBUG
108 :
109 978 : if(iexit/=0)then
110 33 : ABI_SFREE(fcart_m)
111 33 : ABI_SFREE(vel_nexthalf)
112 : return
113 : end if
114 :
115 : !write(std_out,*) 'isokinetic 01'
116 : !##########################################################
117 : !### 01. Debugging and Verbose
118 :
119 945 : if(zDEBUG)then
120 0 : write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),&
121 0 : & 'Debugging and Verbose for pred_isokinetic',('-',kk=1,37)
122 0 : write(std_out,*) 'ionmov: ',12
123 0 : write(std_out,*) 'itime: ',itime
124 : end if
125 :
126 : !write(std_out,*) 'isokinetic 02'
127 : !##########################################################
128 : !### 02. Allocate the vectors vin, vout and hessian matrix
129 : !### These arrays could be allocated from a previous
130 : !### dataset that exit before itime==ntime
131 :
132 945 : if(itime==1)then
133 33 : ABI_SFREE(fcart_m)
134 33 : ABI_SFREE(vel_nexthalf)
135 : end if
136 :
137 945 : if (.not.allocated(fcart_m)) then
138 99 : ABI_MALLOC(fcart_m,(3,ab_mover%natom))
139 : end if
140 945 : if (.not.allocated(vel_nexthalf)) then
141 99 : ABI_MALLOC(vel_nexthalf,(3,ab_mover%natom))
142 : end if
143 :
144 : !write(std_out,*) 'isokinetic 03'
145 : !##########################################################
146 : !### 03. Obtain the present values from the history
147 :
148 945 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
149 :
150 119905 : fcart(:,:)=hist%fcart(:,:,hist%ihist)
151 6615 : strten(:) =hist%strten(:,hist%ihist)
152 119905 : vel(:,:) =hist%vel(:,:,hist%ihist)
153 945 : etotal =hist%etot(hist%ihist)
154 :
155 945 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
156 :
157 945 : if(zDEBUG)then
158 0 : write (std_out,*) 'fcart:'
159 0 : do kk=1,ab_mover%natom
160 0 : write (std_out,*) fcart(:,kk)
161 : end do
162 0 : write (std_out,*) 'vel:'
163 0 : do kk=1,ab_mover%natom
164 0 : write (std_out,*) vel(:,kk)
165 : end do
166 0 : write (std_out,*) 'strten:'
167 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
168 0 : write (std_out,*) 'etotal:'
169 0 : write (std_out,*) etotal
170 : end if
171 :
172 : !Count the number of degrees of freedom, taking into account iatfix.
173 : !Also fix the velocity to zero for the fixed atoms
174 945 : nxyzatfree=0
175 30685 : do iatom=1,ab_mover%natom
176 119905 : do idim=1,3
177 118960 : if(ab_mover%iatfix(idim,iatom)==0)then
178 89184 : nxyzatfree=nxyzatfree+1
179 : else
180 36 : vel(idim,iatom)=zero
181 : endif
182 : enddo
183 : enddo
184 :
185 : !Now, the number of degrees of freedom is reduced by four because of the kinetic energy conservation
186 : !and because of the conservation of the total momentum for each dimension, in case no atom position is fixed for that dimension
187 : !(in the latter case, one degree of freedom has already been taken away)
188 : !This was not done until v8.9 of ABINIT ...
189 945 : ndegfreedom=nxyzatfree
190 945 : ndegfreedom=nxyzatfree-1 ! Kinetic energy conservation
191 3780 : do idim=1,3
192 93000 : if(sum(ab_mover%iatfix(idim,:))==0)then
193 2799 : ndegfreedom=ndegfreedom-1 ! Macroscopic momentum
194 : endif
195 : enddo
196 :
197 : !write(std_out,*) 'isokinetic 04'
198 : !##########################################################
199 : !### 04. Second half-velocity (Only after the first itime)
200 :
201 945 : if(itime>1) then
202 :
203 29888 : do iatom=1,ab_mover%natom
204 116816 : do idim=1,3
205 115904 : if(ab_mover%iatfix(idim,iatom)==0)then
206 86896 : fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom)
207 : else
208 32 : fcart_m(idim,iatom)=zero
209 : endif
210 : end do
211 : end do
212 :
213 : ! Computation of vel(:,:) at the next positions
214 : ! Computation of v2gauss, actually twice the kinetic energy.
215 : ! Called 2K, cf Eq. (A13) of [[cite:Minary2003]].
216 912 : v2gauss=0.0_dp
217 29888 : do iatom=1,ab_mover%natom
218 116816 : do idim=1,3
219 : v2gauss=v2gauss+&
220 : & vel_nexthalf(idim,iatom)*vel_nexthalf(idim,iatom)*&
221 115904 : & ab_mover%amass(iatom)
222 : end do
223 : end do
224 :
225 : ! Computation of a and b (4.13 of [[cite:Minary2003]])
226 912 : a=0.0_dp
227 912 : b=0.0_dp
228 29888 : do iatom=1,ab_mover%natom
229 116816 : do idim=1,3
230 86928 : a=a+fcart_m(idim,iatom)*vel_nexthalf(idim,iatom)*ab_mover%amass(iatom)
231 115904 : b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom)
232 : end do
233 : end do
234 912 : a=a/v2gauss
235 912 : b=b/v2gauss
236 :
237 :
238 : ! Computation of s and scdot
239 912 : sqb=sqrt(b)
240 912 : as=sqb*ab_mover%dtion/2.
241 : ! jmb
242 : if ( as > 300.0 ) as=300.0
243 912 : s1=cosh(as)
244 912 : s2=sinh(as)
245 912 : s=a*(s1-1.)/b+s2/sqb
246 912 : scdot=a*s2/sqb+s1
247 :
248 29888 : do iatom=1,ab_mover%natom
249 116816 : do idim=1,3
250 115904 : if(ab_mover%iatfix(idim,iatom)==0)then
251 86896 : vel(idim,iatom)=(vel_nexthalf(idim,iatom)+fcart_m(idim,iatom)*s)/scdot
252 : else
253 32 : vel(idim,iatom)=zero
254 : endif
255 : enddo
256 : enddo
257 :
258 912 : if (zDEBUG)then
259 0 : write(std_out,*) 'Computation of the second half-velocity'
260 0 : write(std_out,*) 'Cartesian forces per atomic mass (fcart_m):'
261 0 : do kk=1,ab_mover%natom
262 0 : write (std_out,*) fcart_m(:,kk)
263 : end do
264 0 : write(std_out,*) 'vel:'
265 0 : do kk=1,ab_mover%natom
266 0 : write (std_out,*) vel(:,kk)
267 : end do
268 0 : write(std_out,*) 'v2gauss:',v2gauss
269 0 : write(std_out,*) 'a:',a
270 0 : write(std_out,*) 'b:',b
271 0 : write(std_out,*) 's:',s
272 0 : write(std_out,*) 'scdot:',scdot
273 : end if
274 :
275 : end if ! (if itime>1)
276 :
277 : !write(std_out,*) 'isokinetic 05'
278 : !##########################################################
279 : !### 05. First half-time (First cycle the loop is double)
280 :
281 945 : if (itime==1) then
282 : nfirst=2
283 : else
284 912 : nfirst=1
285 : end if
286 :
287 1923 : do ifirst=1,nfirst
288 :
289 : ! Application of Gauss' principle of least constraint according to Fei Zhang's algorithm (J. Chem. Phys. 106, 1997, [[cite:Zhang1997]] p.6102)
290 :
291 : ! v2gauss is twice the kinetic energy
292 978 : v2gauss=0.0_dp
293 31482 : do iatom=1,ab_mover%natom
294 122994 : do idim=1,3
295 122016 : v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
296 : end do
297 : end do
298 :
299 : ! If there is no kinetic energy to start with ...
300 978 : if (v2gauss<=v2tol.and.itime==1) then
301 : ! Maxwell-Boltzman distribution
302 31 : v2gauss=zero
303 31 : vtest=zero
304 731 : do iatom=1,ab_mover%natom
305 2831 : do idim=1,3
306 2800 : if(ab_mover%iatfix(idim,iatom)==0)then
307 2097 : vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum))
308 2097 : vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum)))
309 : else
310 3 : vel(idim,iatom)=zero
311 : endif
312 : end do
313 : end do
314 :
315 : ! Get rid of center-of-mass velocity
316 731 : s1=sum(ab_mover%amass(:))
317 124 : do idim=1,3
318 2224 : if(sum(ab_mover%iatfix(idim,:))==0)then
319 2094 : s2=sum(ab_mover%amass(:)*vel(idim,:))
320 2094 : vel(idim,:)=vel(idim,:)-s2/s1
321 : endif
322 : end do
323 :
324 : ! Recompute v2gauss
325 731 : do iatom=1,ab_mover%natom
326 2831 : do idim=1,3
327 2100 : v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
328 2800 : vtest=vtest+vel(idim,iatom)/ndegfreedom
329 : end do
330 : end do
331 :
332 : ! Now rescale the velocities to give the exact temperature
333 31 : rescale_vel=sqrt(ndegfreedom*kb_HaK*ab_mover%mdtemp(1)/v2gauss)
334 2831 : vel(:,:)=vel(:,:)*rescale_vel
335 :
336 : ! Recompute v2gauss with the rescaled velocities
337 31 : v2gauss=zero
338 731 : do iatom=1,ab_mover%natom
339 2831 : do idim=1,3
340 2800 : v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
341 : end do
342 : end do
343 :
344 : ! Compute the variance and print
345 31 : sigma2=(v2gauss/ndegfreedom-ab_mover%amass(1)*vtest**2)/kb_HaK
346 :
347 : end if
348 :
349 31482 : do iatom=1,ab_mover%natom
350 122994 : do idim=1,3
351 122016 : if(ab_mover%iatfix(idim,iatom)==0)then
352 91472 : fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom)
353 : else
354 40 : fcart_m(idim,iatom)=zero
355 : endif
356 : end do
357 : end do
358 :
359 978 : if (zDEBUG)then
360 0 : write(std_out,*) 'Calculation first half-velocity '
361 0 : write (std_out,*) 'vel:'
362 0 : do kk=1,ab_mover%natom
363 0 : write (std_out,*) vel(:,kk)
364 : end do
365 0 : write (std_out,*) 'xcart:'
366 0 : do kk=1,ab_mover%natom
367 0 : write (std_out,*) xcart(:,kk)
368 : end do
369 0 : write (std_out,*) 'xred:'
370 0 : do kk=1,ab_mover%natom
371 0 : write (std_out,*) xred(:,kk)
372 : end do
373 0 : write (std_out,*) 'fcart_m'
374 0 : do kk=1,ab_mover%natom
375 0 : write (std_out,*) fcart_m(:,kk)
376 : end do
377 0 : write(std_out,*) 's2',s2
378 0 : write(std_out,*) 'v2gauss',v2gauss
379 0 : write(std_out,*) 'sigma2',sigma2
380 :
381 : write(message, '(a)' )&
382 0 : & ' --- Rescaling or initializing velocities to initial temperature'
383 0 : call wrtout(std_out,message,'COLL')
384 : write(message, '(a,d12.5,a,D12.5)' )&
385 0 : & ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1)
386 0 : call wrtout(std_out,message,'COLL')
387 : write(message, '(a,d12.5,a,D12.5)' )&
388 0 : & ' --- Effective temperature',v2gauss/(ndegfreedom*kb_HaK),' From variance', sigma2
389 0 : call wrtout(std_out,message,'COLL')
390 : end if
391 :
392 : ! Convert input xred (reduced coordinates) to xcart (cartesian)
393 978 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
394 :
395 1923 : if(itime==1.and.ifirst==1) then
396 33 : call wrtout(std_out,'if itime==1','COLL')
397 3089 : vel_nexthalf(:,:)=vel(:,:)
398 3089 : xcart_next(:,:)=xcart(:,:)
399 33 : call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
400 3089 : xred=xred_next
401 33 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
402 : end if
403 :
404 : end do
405 :
406 : !Computation of vel_nexthalf (4.16 of [[cite:Minary2003]])
407 : !Computation of a and b (4.13 of [[cite:Minary2003]])
408 945 : a=0.0_dp
409 945 : b=0.0_dp
410 30685 : do iatom=1,ab_mover%natom
411 119905 : do idim=1,3
412 89220 : a=a+fcart_m(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
413 118960 : b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom)
414 : end do
415 : end do
416 945 : a=a/v2gauss+tol20
417 945 : b=b/v2gauss+tol20
418 : !Computation of s and scdot
419 945 : sqb=sqrt(b)+tol20
420 945 : as=sqb*ab_mover%dtion/2.
421 : ! jmb
422 : if ( as > 300.0 ) as=300.0
423 945 : s1=cosh(as)
424 945 : s2=sinh(as)
425 945 : s=a*(s1-1.)/b+s2/sqb
426 945 : scdot=a*s2/sqb+s1
427 30685 : do iatom=1,ab_mover%natom
428 119905 : do idim=1,3
429 118960 : if(ab_mover%iatfix(idim,iatom)==0)then
430 89184 : vel_nexthalf(idim,iatom)=(vel(idim,iatom)+fcart_m(idim,iatom)*s)/scdot
431 : else
432 36 : vel_nexthalf(idim,iatom)=zero
433 : endif
434 : enddo
435 : enddo
436 :
437 : !Computation of the next positions
438 119905 : xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion
439 :
440 945 : if (zDEBUG)then
441 0 : write(std_out,*) 'a:',a
442 0 : write(std_out,*) 'b:',b
443 0 : write(std_out,*) 's:',s
444 0 : write(std_out,*) 'scdot:',scdot
445 : end if
446 :
447 : !Convert back to xred (reduced coordinates)
448 :
449 945 : call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
450 :
451 : !write(std_out,*) 'isokinetic 06'
452 : !##########################################################
453 : !### 06. Update the history with the prediction
454 :
455 119905 : xcart=xcart_next
456 119905 : xred=xred_next
457 :
458 : !increment the ihist
459 945 : hist%ihist = abihist_findIndex(hist,+1)
460 :
461 : !Fill the history with the variables
462 : !xred, acell, rprimd, vel
463 945 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
464 119905 : hist%vel(:,:,hist%ihist)=vel(:,:)
465 945 : hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
466 :
467 945 : if(zDEBUG)then
468 0 : write (std_out,*) 'fcart:'
469 0 : do kk=1,ab_mover%natom
470 0 : write (std_out,*) fcart(:,kk)
471 : end do
472 0 : write (std_out,*) 'vel:'
473 0 : do kk=1,ab_mover%natom
474 0 : write (std_out,*) vel(:,kk)
475 : end do
476 0 : write (std_out,*) 'strten:'
477 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
478 0 : write (std_out,*) 'etotal:'
479 0 : write (std_out,*) etotal
480 : end if
481 :
482 : if (.false.) write(std_out,*) ntime
483 :
484 : end subroutine pred_isokinetic
485 : !!***
486 :
487 : end module m_pred_isokinetic
488 : !!***
|
