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