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