Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_nose
2 : : !! NAME
3 : : !! m_pred_nose
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_nose
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_abimover
27 : : use m_abihist
28 : :
29 : : use m_numeric_tools, only : uniformrandom
30 : : use m_geometry, only : xcart2xred, xred2xcart, metric
31 : :
32 : : implicit none
33 : :
34 : : private
35 : : !!***
36 : :
37 : : public :: pred_nose
38 : : !!***
39 : :
40 : : contains
41 : : !!***
42 : :
43 : : !!****f* ABINIT/pred_nose
44 : : !! NAME
45 : : !! pred_nose
46 : : !!
47 : : !! FUNCTION
48 : : !! Ionmov predictors (8) Verlet algorithm with a nose-hoover thermostat
49 : : !!
50 : : !! IONMOV 8:
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),
55 : : !! the Verlet dynamics is performed, using the gradient of the
56 : : !! energy (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 "ntime" steps are performed.
61 : : !! The time step is governed by dtion (contained in dtset)
62 : : !! Returned quantities are xred, and eventually acell and rprimd
63 : : !! (new ones!).
64 : : !!
65 : : !! See ionmov=6, but with a nose-hoover thermostat
66 : : !! Velocity verlet algorithm : Swope et al JCP 76 (1982) 637
67 : : !!
68 : : !! INPUTS
69 : : !! ab_mover <type(abimover)> : Datatype with all the information
70 : : !! needed by the preditor
71 : : !! itime : Index of the present iteration
72 : : !! ntime : Maximal number of iterations
73 : : !! zDEBUG : if true print some debugging information
74 : : !!
75 : : !! SIDE EFFECTS
76 : : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
77 : : !!
78 : : !! SOURCE
79 : :
80 : 20 : subroutine pred_nose(ab_mover,hist,itime,ntime,zDEBUG,iexit)
81 : :
82 : : !Arguments ------------------------------------
83 : : !scalars
84 : : type(abimover),intent(in) :: ab_mover
85 : : type(abihist),intent(inout) :: hist
86 : : integer,intent(in) :: itime
87 : : integer,intent(in) :: ntime
88 : : integer,intent(in) :: iexit
89 : : logical,intent(in) :: zDEBUG
90 : :
91 : : !Local variables-------------------------------
92 : : !scalars
93 : : integer :: ii,jj,kk
94 : : integer :: idum=-5
95 : : real(dp),parameter :: v2tol=tol8,nosetol=tol10
96 : : real(dp) :: delxi,xio,ktemp,rescale_vel
97 : : real(dp) :: dnose,v2nose,xin_nose
98 : : real(dp),save :: xi_nose,fsnose,snose
99 : : real(dp) :: gnose
100 : : real(dp) :: ucvol,ucvol_next
101 : : real(dp) :: etotal
102 : : logical :: ready
103 : :
104 : : !arrays
105 : : real(dp) :: acell(3),acell_next(3)
106 : : real(dp) :: rprimd(3,3),rprimd_next(3,3)
107 : : real(dp) :: gprimd(3,3)
108 : : real(dp) :: gmet(3,3)
109 : : real(dp) :: rmet(3,3)
110 : 40 : real(dp) :: fcart(3,ab_mover%natom)
111 : 40 : real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
112 : 40 : real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
113 : 40 : real(dp) :: vel(3,ab_mover%natom),vel_temp(3,ab_mover%natom)
114 : 40 : real(dp) :: finose(3,ab_mover%natom),binose(3,ab_mover%natom)
115 : 2 : real(dp) :: vonose(3,ab_mover%natom),hnose(3,ab_mover%natom)
116 : : real(dp),allocatable,save :: fcart_m(:,:),fcart_mold(:,:)
117 : : real(dp) :: strten(6)
118 : : character(len=500) :: message
119 : :
120 : : !***************************************************************************
121 : : !Beginning of executable session
122 : : !***************************************************************************
123 : :
124 [ + + ]: 20 : if(iexit/=0)then
125 [ + - ]: 2 : ABI_SFREE(fcart_m)
126 [ + - ]: 2 : ABI_SFREE(fcart_mold)
127 : : return
128 : : end if
129 : :
130 : : !write(std_out,*) 'nose 01'
131 : : !##########################################################
132 : : !### 01. Allocate the arrays fcart_m and fcart_mold
133 : :
134 [ + + ]: 18 : if(itime==1)then
135 [ - + ]: 2 : ABI_SFREE(fcart_m)
136 [ - + ]: 2 : ABI_SFREE(fcart_mold)
137 : : end if
138 : :
139 [ + + ]: 18 : if(.not.allocated(fcart_m)) then
140 [ - + + - : 6 : ABI_MALLOC(fcart_m,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
141 : : end if
142 [ + + ]: 18 : if(.not.allocated(fcart_mold)) then
143 [ - + + - : 6 : ABI_MALLOC(fcart_mold,(3,ab_mover%natom))
+ - + - -
+ - + -
+ ]
144 : : end if
145 : :
146 : : !write(std_out,*) 'nose 02'
147 : : !##########################################################
148 : : !### 02. Obtain the present values from the history
149 : :
150 : 18 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
151 : 18 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
152 : :
153 [ + + + + ]: 198 : fcart(:,:)=hist%fcart(:,:,hist%ihist)
154 : : strten(:)=hist%strten(:,hist%ihist)
155 [ + + + + ]: 198 : vel(:,:)=hist%vel(:,:,hist%ihist)
156 : 18 : etotal=hist%etot(hist%ihist)
157 : :
158 : 18 : write(std_out,*) 'RPRIMD'
159 [ + + ]: 72 : do ii=1,3
160 : 72 : write(std_out,*) rprimd(:,ii)
161 : : end do
162 : 18 : write(std_out,*) 'RMET'
163 [ + + ]: 72 : do ii=1,3
164 : 72 : write(std_out,*) rmet(ii,:)
165 : : end do
166 : :
167 : : !write(std_out,*) 'nose 03'
168 : : !##########################################################
169 : : !### 03. Fill the vectors vin and vout
170 : :
171 : : !write(std_out,*) 'nose 04'
172 : : !##########################################################
173 : : !### 04. Initialize or update the hessian matrix
174 : :
175 : : !write(std_out,*) 'nose 05'
176 : : !##########################################################
177 : : !### 05. Compute the next values
178 : :
179 : : !The temperature is linear between initial and final values
180 : : !It is here converted from Kelvin to Hartree (kb_HaK)
181 : 18 : ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK
182 : :
183 : : !%%% NOSE DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184 : :
185 : : acell_next(:)=acell(:)
186 : 18 : ucvol_next=ucvol
187 : : rprimd_next(:,:)=rprimd(:,:)
188 : :
189 [ + + ]: 18 : if(itime==1)then
190 : 2 : snose=0.0_dp
191 : 2 : xi_nose=0.0_dp
192 : : ! Compute twice the kinetic energy of the system, called v2nose
193 : 2 : v2nose=0.0_dp
194 [ + + ]: 7 : do kk=1,ab_mover%natom
195 [ + + ]: 22 : do jj=1,3
196 : 20 : v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
197 : : end do
198 : : end do
199 [ - + ]: 2 : if (zDEBUG)then
200 : 0 : write(std_out,*) 'itime ntime KTEMP=',itime-1,ntime-1,ktemp
201 : 0 : write(std_out,*) 'V2NOSE=',v2nose
202 : 0 : write (std_out,*) 'VEL'
203 [ # # ]: 0 : do kk=1,ab_mover%natom
204 : 0 : write (std_out,*) vel(:,kk)
205 : : end do
206 : : end if
207 : :
208 : : ! If there is no kinetic energy, use a random initial velocity
209 [ + - ]: 2 : if (v2nose<=v2tol) then
210 : 2 : v2nose=0.0_dp
211 [ + + ]: 7 : do kk=1,ab_mover%natom
212 [ + + ]: 22 : do jj=1,3
213 : : ! Uniform random returns a uniform random deviate between 0.0
214 : : ! and 1.0
215 : : ! if it were always 0 or 1, then the following expression
216 : : ! would give the requested temperature
217 : : vel(jj,kk)=(1.0_dp-2.0_dp*uniformrandom(idum))*&
218 : 15 : & sqrt( (ab_mover%mdtemp(1)) * kb_HaK / ab_mover%amass(kk) )
219 : : ! Recompute v2nose
220 : 15 : v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
221 [ - + ]: 20 : if (zDEBUG)then
222 : 0 : write(std_out,*) 'jj kk vel(jj,kk)=',jj,kk,vel(jj,kk)
223 : 0 : write(std_out,*) 'jj kk V2NOSE=',jj,kk,v2nose
224 : : end if
225 : : end do
226 : : end do
227 : : end if
228 : 2 : write(std_out,*) 'V2NOSE=',v2nose
229 : :
230 : : ! Now, rescale the velocities to give the proper temperature
231 : 2 : rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose)
232 : 2 : write(std_out,*) 'RESCALE_VEL=',rescale_vel
233 [ + + + + ]: 22 : vel(:,:)=vel(:,:)*rescale_vel
234 : : ! Recompute v2nose with the rescaled velocities
235 : 2 : v2nose=0.0_dp
236 [ + + ]: 7 : do kk=1,ab_mover%natom
237 [ + + ]: 22 : do jj=1,3
238 : 20 : v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
239 : : end do
240 : : end do
241 : : write(message, '(a)' )&
242 : 2 : & ' Rescaling or initializing velocities to initial temperature'
243 : 2 : call wrtout(std_out,message,'COLL')
244 : 2 : call wrtout(std_out,message,'COLL')
245 : : write(message, '(2(a,es22.14))' )&
246 : 2 : & ' --- Scaling factor : ',rescale_vel,&
247 : 4 : & ' Asked T (K) ',ab_mover%mdtemp(1)
248 : 2 : call wrtout(std_out,message,'COLL')
249 : 2 : call wrtout(std_out,message,'COLL')
250 : : write(message, '(a,es22.14)' )&
251 : 2 : & ' --- Effective temperature',v2nose/(3.0_dp*ab_mover%natom*kb_HaK)
252 : 2 : call wrtout(std_out,message,'COLL')
253 : 2 : call wrtout(std_out,message,'COLL')
254 : : end if
255 : :
256 [ + + ]: 63 : do kk=1,ab_mover%natom
257 [ + + ]: 198 : fcart_m(:,kk)=fcart(:,kk)/ab_mover%amass(kk)
258 : : end do
259 : :
260 : : !First step of velocity verlet algorithm
261 : 18 : gnose=3*ab_mover%natom
262 : :
263 : : !Convert input xred (reduced coordinates) to xcart (cartesian)
264 : 18 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
265 : :
266 : : !Calculate nose-hoover force on atoms
267 : : !If first iteration, no old force are available, so use present
268 : : !forces
269 [ + + + + : 38 : if (itime==1) fcart_mold(:,:)=fcart_m(:,:)
+ + ]
270 : :
271 [ - + ]: 18 : if (zDEBUG)then
272 : 0 : write (std_out,*) 'FCART_MOLD'
273 [ # # ]: 0 : do kk=1,ab_mover%natom
274 : 0 : write (std_out,*) fcart_mold(:,kk)
275 : : end do
276 : 0 : write (std_out,*) 'FCART_M'
277 [ # # ]: 0 : do kk=1,ab_mover%natom
278 : 0 : write (std_out,*) fcart_m(:,kk)
279 : : end do
280 : : end if
281 : :
282 [ + + + + ]: 198 : finose(:,:)=fcart_mold(:,:)-xi_nose*vel(:,:)
283 [ + + + + ]: 198 : xcart(:,:)=xcart(:,:)+ab_mover%dtion*(vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp)
284 : :
285 : : !Convert back to xred (reduced coordinates)
286 : 18 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
287 : :
288 [ - + ]: 18 : if (zDEBUG)then
289 : 0 : write (std_out,*) 'VEL'
290 [ # # ]: 0 : do kk=1,ab_mover%natom
291 : 0 : write (std_out,*) vel(:,kk)
292 : : end do
293 : : end if
294 : :
295 : : !Calculate v2nose
296 : 18 : v2nose=0.0_dp
297 [ + + ]: 63 : do kk=1,ab_mover%natom
298 [ + + ]: 198 : do jj=1,3
299 : 180 : v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
300 : : end do
301 : : end do
302 [ + + + + ]: 198 : vel(:,:)=vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp
303 : :
304 [ - + ]: 18 : if (zDEBUG)then
305 : 0 : write(std_out,*) 'NOSE BEFORE'
306 : 0 : write(std_out,*) 'FSNOSE=',fsnose
307 : 0 : write(std_out,*) 'SNOSE=',snose
308 : 0 : write(std_out,*) 'XI_NOSE=',xi_nose
309 : 0 : write (std_out,*) 'VEL'
310 [ # # ]: 0 : do kk=1,ab_mover%natom
311 : 0 : write (std_out,*) vel(:,kk)
312 : : end do
313 : 0 : write (std_out,*) 'NOSEINERT',ab_mover%noseinert
314 : : end if
315 : :
316 : : !Update thermostat
317 : 18 : fsnose=(v2nose-gnose*ktemp)/ab_mover%noseinert
318 : 18 : snose=snose+ab_mover%dtion*(xi_nose+ab_mover%dtion*fsnose/2.0_dp)
319 : 18 : xi_nose=xi_nose+ab_mover%dtion*fsnose/2.0_dp
320 [ - + ]: 18 : if (zDEBUG)then
321 : 0 : write(std_out,*) 'NOSE AFTER'
322 : 0 : write(std_out,*) 'FSNOSE=',fsnose
323 : 0 : write(std_out,*) 'SNOSE=',snose
324 : 0 : write(std_out,*) 'XI_NOSE=',xi_nose
325 : 0 : write (std_out,*) 'VEL'
326 [ # # ]: 0 : do kk=1,ab_mover%natom
327 : 0 : write (std_out,*) vel(:,kk)
328 : : end do
329 : : end if
330 : :
331 : : !Second step of the velocity Verlet algorithm, uses the 'new forces'
332 : : !Calculate v2nose
333 : 18 : v2nose=0.0_dp
334 [ + + ]: 63 : do kk=1,ab_mover%natom
335 [ + + ]: 198 : do jj=1,3
336 : 180 : v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
337 : : end do
338 : : end do
339 [ + + + + ]: 198 : vel_temp(:,:)=vel(:,:)
340 : :
341 [ - + ]: 18 : if (zDEBUG)then
342 : 0 : write(std_out,*) 'V2NOSE=',v2nose
343 : 0 : write (std_out,*) 'VEL'
344 [ # # ]: 0 : do kk=1,ab_mover%natom
345 : 0 : write (std_out,*) vel(:,kk)
346 : : end do
347 : 0 : write (std_out,*) 'Starting Newton Raphson'
348 : : end if
349 : :
350 : 18 : xin_nose=xi_nose
351 : :
352 : : !Start Newton-Raphson loop
353 : 18 : ready=.false.
354 [ + + ]: 88 : do while (.not.ready)
355 : : xio=xin_nose
356 : 244 : delxi=0.0D0
357 [ + + + + ]: 766 : vonose(:,:)=vel_temp(:,:)
358 : : hnose(:,:)=-ab_mover%dtion/2.0_dp*(fcart_m(:,:)-xio*vonose(:,:))-&
359 [ + + + + ]: 766 : & (vel(:,:)-vonose(:,:))
360 [ + + ]: 244 : do kk=1,ab_mover%natom
361 [ + + ]: 766 : do jj=1,3
362 : 522 : binose(jj,kk)=vonose(jj,kk)*ab_mover%dtion/ab_mover%noseinert*ab_mover%amass(kk) ! a verifier
363 : 696 : delxi=delxi+hnose(jj,kk)*binose(jj,kk)
364 : : end do
365 : : end do
366 : 70 : dnose=-(xio*ab_mover%dtion/2.0D0+1.0D0)
367 : : delxi=delxi-dnose*((-v2nose+gnose*ktemp)*ab_mover%dtion/2.0_dp/ &
368 : 70 : & ab_mover%noseinert-(xi_nose-xio))
369 : 70 : delxi=delxi/(-ab_mover%dtion*ab_mover%dtion/2.0_dp*v2nose/ab_mover%noseinert+dnose)
370 : :
371 : : ! hzeronose=-(xio-xi_nose-(v2nose-gnose*ktemp)
372 : : ! *dtion/(2.0_dp*ab_mover%noseinert) )
373 : : ! cibinose=-v2nose*dtion*dtion/(2.0_dp*ab_mover%noseinert)
374 : : ! delxi=(delxi+hzeronose*dnose)/(dnose+cibinose)
375 : :
376 : : ! DEBUG
377 : : ! write(message, '(a,es22.14)' )' after delxi',delxi
378 : : ! call wrtout(std_out,message,'COLL')
379 : : ! call wrtout(std_out,message,'COLL')
380 : : ! ENDDEBUG
381 : 70 : v2nose=0.0_dp
382 : :
383 : : vel_temp(:,:)=vel_temp(:,:)+&
384 [ + + + + ]: 766 : & (hnose+ab_mover%dtion/2.0_dp*vonose(:,:)*delxi)/dnose
385 [ + + ]: 244 : do kk=1,ab_mover%natom
386 [ + + ]: 766 : do jj=1,3
387 : : v2nose=v2nose+vel_temp(jj,kk)*&
388 : 696 : & vel_temp(jj,kk)*ab_mover%amass(kk)
389 : : end do
390 : : end do
391 : : ! New guess for xi
392 : 70 : xin_nose=xio+delxi
393 : :
394 : : ! zDEBUG
395 : : ! write(message, '(a,es22.14)' )' v2nose=',v2nose
396 : : ! call wrtout(std_out,message,'COLL')
397 : : ! call wrtout(std_out,message,'COLL')
398 : : ! ENDDEBUG
399 : :
400 : 70 : ready=.true.
401 : : ! Test for convergence
402 : 70 : kk=0
403 : 70 : jj=1
404 [ - + + + : 302 : do while((kk<=ab_mover%natom).and.(jj<=3).and.ready)
+ + ]
405 : 214 : kk=kk+1
406 [ + + ]: 214 : if (kk>ab_mover%natom) then
407 : 57 : kk=1
408 : 57 : jj=jj+1
409 : : end if
410 [ + - + + ]: 284 : if ((kk<=ab_mover%natom) .and.(jj<=3)) then
411 [ - + ]: 195 : if (abs(vel_temp(jj,kk))<1.0d-50)&
412 : 0 : & vel_temp(jj,kk)=1.0d-50
413 [ + + ]: 195 : if (abs((vel_temp(jj,kk)-vonose(jj,kk))&
414 : 51 : & /vel_temp(jj,kk))>nosetol) ready=.false.
415 : : else
416 : : if (xin_nose<1.0d-50) xin_nose=1.0d-50
417 [ + + ]: 19 : if (abs((xin_nose-xio)/xin_nose)>nosetol) ready=.false.
418 : : end if
419 : : end do ! end of while
420 : :
421 : : ! Enddo ready
422 : : end do
423 : :
424 : : !Update velocities to converged value
425 [ + + + + ]: 198 : vel(:,:)=vel_temp(:,:)
426 : 18 : write(message, '(a,es14.7)' )' converged velocities for T=',ktemp
427 : 18 : call wrtout(std_out,message,'COLL')
428 : :
429 [ - + ]: 18 : if (zDEBUG)then
430 : 0 : write (std_out,*) 'Final Values for NOSE'
431 : 0 : write (std_out,*) 'VEL'
432 [ # # ]: 0 : do kk=1,ab_mover%natom
433 : 0 : write (std_out,*) vel(:,kk)
434 : : end do
435 : 0 : write (std_out,*) 'XCART'
436 [ # # ]: 0 : do kk=1,ab_mover%natom
437 : 0 : write (std_out,*) xcart(:,kk)
438 : : end do
439 : : end if
440 : :
441 : : !Update thermostat
442 : 18 : xi_nose=xin_nose
443 [ + + + + ]: 198 : xcart_next(:,:)=xcart(:,:)
444 : : !Convert back to xred_next (reduced coordinates)
445 : 18 : call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
446 : : !Store 'new force' as 'old force'
447 [ + + + + ]: 198 : fcart_mold(:,:)=fcart_m(:,:)
448 : :
449 : : !write(std_out,*) 'nose 06'
450 : : !##########################################################
451 : : !### 06. Update the history with the prediction
452 : :
453 : : !Increase indexes
454 : 18 : hist%ihist=abihist_findIndex(hist,+1)
455 : :
456 : 18 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
457 [ + + + + ]: 198 : hist%vel(:,:,hist%ihist)=vel(:,:)
458 : 18 : hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
459 : :
460 : : end subroutine pred_nose
461 : : !!***
462 : :
463 : : end module m_pred_nose
464 : : !!***
|
