Line data Source code
1 : !!****m* ABINIT/m_pimd_nosehoover
2 : !! NAME
3 : !! m_pimd_nosehoover
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2011-2026 ABINIT group (GG,MT)
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_pimd_nosehoover
23 :
24 : use defs_basis
25 : use m_pimd
26 : use m_abicore
27 :
28 : use m_geometry, only : xcart2xred, xred2xcart
29 :
30 : implicit none
31 :
32 : private
33 : !!***
34 :
35 : public :: pimd_nosehoover_npt
36 : public :: pimd_nosehoover_nvt
37 : !!***
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* ABINIT/pimd_nosehoover_npt
43 : !! NAME
44 : !! pimd_nosehoover_npt
45 : !!
46 : !! FUNCTION
47 : !! Predicts new positions in Path Integral Molecular Dynamics using Nose-Hoover in the NPT ensemble.
48 : !! Given the positions at time t and t-dtion, an estimation of the velocities at time t,
49 : !! the forces and an estimation of the stress at time t, and an estimation of the cell at time t,
50 : !! computes in the Path Integral Molecular Dynamics framework the new positions at time t+dtion,
51 : !! computes self-consistently the velocities, the stress and the cell at time t and produces
52 : !! an estimation of the velocities, stress and new cell at time t+dtion
53 : !! No change of acell and rprim at present.
54 : !!
55 : !! INPUTS
56 : !! etotal(trotter)=electronic total energy for all images
57 : !! itimimage=number of the current time for image propagation (itimimage+1 is to be predicted here)
58 : !! natom=dimension of vel_timimage and xred_timimage
59 : !! pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
60 : !! prtvolimg=printing volume
61 : !! rprimd(3,3)=dimensionless unit cell vectors (common to all images) at time t (present time step)
62 : !! rprimd_prev(3,3)=dimensionless unit cell vectors (common to all images) at time t-dt (previous time step)
63 : !! stressin(3,3,trotter)=electronic stress tensor for each image
64 : !! trotter=Trotter number (total number of images)
65 : !! volume=voume of unit cell (common to all images)
66 : !! xred(3,natom,trotter)=reduced coordinates of atoms for all images at time t (present time step)
67 : !! xred_prev(3,natom,trotter)=reduced coordinates of atoms for all images at time t-dt (previous time step)
68 : !!
69 : !! OUTPUT
70 : !! rprimd_next(3,3)=dimensionless unit cell vectors (common to all images) at time t+dt (next time step)
71 : !! xred_next(3,natom,trotter)=reduced coordinates of atoms for all images at time t+dt (next time step)
72 : !!
73 : !! SIDE EFFECTS
74 : !! forces(3,natom,trotter)=forces over atoms for all images
75 : !! at input, electronic forces
76 : !! at output, electronic forces + Langevin contribution
77 : !! vel(3,natom,trotter)/vel_next(3,natom,trotter)/=velocies of atoms for all images
78 : !! at input, vel(itimimage) = values of the estimated vel at t
79 : !! at output, vel(itimimage) = values of the exact vel at t
80 : !! vel_next(itimimage) = values of the estimated vel at t+dt
81 : !! vel_cell(3,3)/vel_cell_next(3,3)=time derivative of cell parameters
82 : !! at input, vel_cell(itimimage) = values of the estimated cell vel at t
83 : !! at output, vel_cell(itimimage) = values of the exact cell vel at t
84 : !! vel_cell_next(itimimage) = values of the estimated cell vel at t+dt
85 : !!
86 : !! SOURCE
87 :
88 0 : subroutine pimd_nosehoover_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
89 0 : & rprimd,rprimd_next,rprimd_prev,stressin,trotter,vel,vel_next,vel_cell,&
90 0 : & vel_cell_next,volume,xred,xred_next,xred_prev)
91 :
92 : !Arguments ------------------------------------
93 : !scalars
94 : integer,intent(in) :: itimimage,natom,prtvolimg,trotter
95 : real(dp),intent(in) :: volume
96 : type(pimd_type),intent(in) :: pimd_param
97 : !arrays
98 : real(dp),intent(in) :: etotal(trotter),rprimd(3,3),rprimd_prev(3,3),stressin(3,3,trotter)
99 : real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
100 : real(dp),intent(out) :: rprimd_next(3,3),xred_next(3,natom,trotter)
101 : real(dp),intent(out) :: vel_next(3,natom,trotter),vel_cell_next(3,3)
102 : real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter),vel_cell(3,3)
103 :
104 : !Local variables-------------------------------
105 : !Options
106 : real(dp),parameter :: tolerance=tol7 ! SCF tolerance
107 : !scalars
108 : integer :: idum=-5
109 : integer :: constraint,iimage,irestart,ndof,nnos,pitransform,prtstress
110 : real(dp) :: dtion,eharm,eharm2,epot,initemp,kt,temperature1,temperature2,thermtemp
111 : !arrays
112 : real(dp) :: constraint_output(2),ddh(3,3),stress_pimd(3,3,3)
113 0 : real(dp),allocatable :: dzeta(:,:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:)
114 0 : real(dp),allocatable :: inertmass(:),masseff(:,:),qmass(:),quantummass(:),springeff(:,:)
115 0 : real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:),zeta(:)
116 :
117 : ! *************************************************************************
118 :
119 : !############# Initializations ###########################
120 :
121 : !Allocation of local arrays
122 0 : ABI_MALLOC(xcart,(3,natom,trotter))
123 0 : ABI_MALLOC(xcart_prev,(3,natom,trotter))
124 0 : ABI_MALLOC(xcart_next,(3,natom,trotter))
125 0 : ABI_MALLOC(forces_orig,(3,natom,trotter))
126 0 : ABI_MALLOC(forces_pimd,(3,natom,trotter))
127 0 : ABI_MALLOC(inertmass,(natom))
128 0 : ABI_MALLOC(quantummass,(natom))
129 :
130 : !Fill in the local variables
131 0 : ndof=3*natom*trotter
132 0 : quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass
133 0 : inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
134 0 : initemp=pimd_param%mdtemp(1);thermtemp=pimd_param%mdtemp(2)
135 0 : dtion=pimd_param%dtion;pitransform=pimd_param%pitransform
136 0 : kt=thermtemp*kb_HaK
137 0 : forces_orig=forces
138 :
139 : !Allocation/initialization of local variables used for Nose-Hoover chains
140 : !Associated variables:
141 : !nnos = number of thermostats
142 : !dzeta(3,natom,trotter,nnos) = variables of thermostats, in (atomic time unit)^(-1)
143 : !qmass(nnos) = masses of thermostats
144 : !specific to PIMD: pitransform = coordinate transformation (0:no; 1:normal mode; 2:staging)
145 0 : nnos=pimd_param%nnos
146 0 : ABI_MALLOC(qmass,(nnos))
147 0 : ABI_MALLOC(zeta,(nnos))
148 0 : ABI_MALLOC(dzeta,(3,natom,trotter,nnos))
149 0 : qmass(1:nnos)=pimd_param%qmass(1:nnos)
150 0 : zeta=zero;dzeta=zero
151 :
152 : !Compute cartesian coordinates
153 0 : do iimage=1,trotter
154 0 : call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage))
155 0 : call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
156 : end do
157 :
158 : !Determine if it is a restart or not
159 : !If this is a calculation from scratch,generate random distribution of velocities
160 0 : irestart=1;if (itimimage==1) irestart=pimd_is_restart(masseff,vel,vel_cell)
161 :
162 : !Initialize derivatives
163 0 : if (mod(irestart,10)==0) then
164 0 : call pimd_initvel(idum,masseff,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
165 : end if
166 : !vel_cell does not depend on Trotter...
167 0 : ddh=vel_cell(:,:);if (irestart<10) ddh=zero
168 :
169 : !Compute temperature at t
170 0 : temperature1=pimd_temperature(masseff,vel)
171 :
172 : !################## Images evolution #####################
173 :
174 : !This is temporary
175 0 : xcart_next=zero
176 0 : rprimd_next=rprimd_prev
177 0 : temperature2=pimd_temperature(masseff,vel)
178 :
179 : !Compute contributions to energy
180 0 : call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,springeff,trotter,xcart)
181 :
182 : !Compute stress tensor at t from virial theorem
183 0 : call pimd_stresses(masseff,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
184 :
185 :
186 : !############# Final operations ############################
187 :
188 : !Print messages
189 0 : prtstress=1
190 : call pimd_print(constraint,constraint_output,&
191 : & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
192 : & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
193 : & stress_pimd,temperature2,&
194 0 : & pimd_param%traj_unit,trotter,vel,ddh,xcart,xred)
195 :
196 : !Come back to reduced coordinates
197 0 : do iimage=1,trotter
198 0 : call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage))
199 : end do
200 :
201 : !Return cell velocities (does not depend on Trotter)
202 0 : vel_cell(:,:)=ddh(:,:)
203 :
204 : !Free memory
205 0 : ABI_FREE(xcart)
206 0 : ABI_FREE(xcart_prev)
207 0 : ABI_FREE(xcart_next)
208 0 : ABI_FREE(forces_orig)
209 0 : ABI_FREE(forces_pimd)
210 0 : ABI_FREE(inertmass)
211 0 : ABI_FREE(quantummass)
212 0 : ABI_FREE(masseff)
213 : ABI_FREE(springeff)
214 : ABI_FREE(qmass)
215 : ABI_FREE(dzeta)
216 : ABI_FREE(zeta)
217 :
218 0 : end subroutine pimd_nosehoover_npt
219 : !!***
220 :
221 : !!****f* ABINIT/pimd_nosehoover_nvt
222 : !! NAME
223 : !! pimd_nosehoover_nvt
224 : !!
225 : !! FUNCTION
226 : !! Predicts new positions in Path Integral Molecular Dynamics using Nose-Hoover in the NVT ensemble.
227 : !! Given the positions at time t and t-dtion, an estimation of the velocities at time t,
228 : !! the forces and an estimation of the stress at time t, and an estimation of the cell at time t,
229 : !! computes in the Path Integral Molecular Dynamics framework the new positions at time t+dtion,
230 : !! computes self-consistently the velocities, the stress and the cell at time t and produces
231 : !! an estimation of the velocities, stress and new cell at time t+dtion
232 : !! No change of acell and rprim at present.
233 : !!
234 : !! INPUTS
235 : !! etotal(trotter)=electronic total energy for all images
236 : !! itimimage=number of the current time for image propagation (itimimage+1 is to be predicted here)
237 : !! natom=dimension of vel_timimage and xred_timimage
238 : !! pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
239 : !! prtvolimg=printing volume
240 : !! rprimd(3,3)=dimensionless unit cell vectors (common to all images)
241 : !! stressin(3,3,trotter)=electronic stress tensor for each image
242 : !! trotter=Trotter number (total number of images)
243 : !! volume=voume of unit cell (common to all images)
244 : !! xred(3,natom,trotter)=reduced coordinates of atoms for all images at time t (present time step)
245 : !! xred_prev(3,natom,trotter)=reduced coordinates of atoms for all images at time t-dt (previous time step)
246 : !!
247 : !! OUTPUT
248 : !! xred_next(3,natom,trotter)=reduced coordinates of atoms for all images at time t+dt (next time step)
249 : !!
250 : !! SIDE EFFECTS
251 : !! forces(3,natom,trotter)=forces over atoms for all images
252 : !! at input, electronic forces
253 : !! at output, electronic forces + quantum spring contribution
254 : !! vel(3,natom,trotter)/vel_next(3,natom,trotter)/=velocies of atoms for all images
255 : !! at input, vel(itimimage) = values of the estimated vel at t
256 : !! at output, vel(itimimage) = values of the exact vel at t
257 : !! vel_next(itimimage) = values of the estimated vel at t+dt
258 : !!
259 : !! NOTES
260 : !! Thermization by Nose-Hoover chains according to
261 : !! Martyna, Klein, Tuckerman, J. Chem. Phys. 97, 2635 (1992) [[cite:Martyna1992]]
262 : !! Tuckerman, Marx, Klein, Parrinello, J. Chem. Phys. 104, 5579 (1996) [[cite:Tuckerman1996]]
263 : !!
264 : !! SOURCE
265 :
266 15 : subroutine pimd_nosehoover_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
267 15 : & rprimd,stressin,trotter,vel,vel_next,volume,xred,xred_next,xred_prev)
268 :
269 : !Arguments ------------------------------------
270 : !scalars
271 : integer,intent(in) :: itimimage,natom,prtvolimg,trotter
272 : real(dp),intent(in) :: volume
273 : type(pimd_type),intent(inout) :: pimd_param
274 : !arrays
275 : real(dp),intent(in) :: etotal(trotter),rprimd(3,3),stressin(3,3,trotter)
276 : real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
277 : real(dp),intent(out) :: xred_next(3,natom,trotter),vel_next(3,natom,trotter)
278 : real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter)
279 :
280 : !Local variables-------------------------------
281 : !Options
282 : real(dp),parameter :: tolerance=tol9 ! SCF tolerance
283 : !scalars
284 : integer :: idum=-5
285 : integer :: iimage,irestart,ndof,nnos,pitransform,prtstress
286 : real(dp) :: dtion,eharm,eharm2,epot,initemp,kt
287 : real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol
288 : character(len=500) :: msg
289 : !arrays
290 30 : real(dp) :: constraint_output(2),spring_prim(natom),stress_pimd(3,3,3),vel_cell(3,3)
291 15 : real(dp),allocatable :: forces_orig(:,:,:),forces_pimd(:,:,:)
292 15 : real(dp),allocatable :: inertmass(:),mass(:,:),qmass(:),quantummass(:),spring(:,:)
293 15 : real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:)
294 15 : real(dp),allocatable :: dzeta(:,:,:,:),zeta_prev(:,:,:,:),zeta(:,:,:,:)
295 15 : real(dp),allocatable :: zeta_next(:,:,:,:)
296 :
297 : ! *************************************************************************
298 :
299 : !############# Initializations ###########################
300 :
301 : !Allocation of local arrays
302 60 : ABI_MALLOC(xcart,(3,natom,trotter))
303 45 : ABI_MALLOC(xcart_prev,(3,natom,trotter))
304 45 : ABI_MALLOC(xcart_next,(3,natom,trotter))
305 45 : ABI_MALLOC(forces_orig,(3,natom,trotter))
306 45 : ABI_MALLOC(forces_pimd,(3,natom,trotter))
307 45 : ABI_MALLOC(inertmass,(natom))
308 30 : ABI_MALLOC(quantummass,(natom))
309 :
310 : !Fill in the local variables
311 15 : ndof=3*natom*trotter
312 15 : pitransform=pimd_param%pitransform
313 75 : quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass
314 75 : inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
315 30 : if(pitransform==1) inertmass=quantummass !compulsory for good definition of normal mode masses
316 30 : if(pitransform==2) inertmass=quantummass !compulsory for good definition of staging masses
317 15 : initemp=pimd_param%mdtemp(1);thermtemp=pimd_param%mdtemp(2)
318 15 : dtion=pimd_param%dtion
319 15 : kt=thermtemp*kb_HaK
320 840 : forces_orig=forces
321 :
322 : !Allocation/initialization of local variables used for Nose-Hoover chains
323 : !Associated variables:
324 : !nnos = number of thermostats
325 : !zeta,zeta_next,zeta_prev(3,natom,trotter,nnos) = variables of thermostats, dzeta, its time derivative
326 : !qmass(nnos) = masses of thermostats
327 : !specific to PIMD: pitransform = coordinate transformation (0:no; 1:normal mode; 2:staging)
328 15 : nnos=pimd_param%nnos
329 45 : ABI_MALLOC(qmass,(nnos))
330 75 : ABI_MALLOC(zeta_prev,(3,natom,trotter,nnos))
331 60 : ABI_MALLOC(zeta,(3,natom,trotter,nnos))
332 60 : ABI_MALLOC(zeta_next,(3,natom,trotter,nnos))
333 60 : ABI_MALLOC(dzeta,(3,natom,trotter,nnos))
334 : !initialization
335 90 : qmass(1:nnos)=pimd_param%qmass(1:nnos)
336 4140 : zeta_prev(:,:,:,:)=pimd_param%zeta_prev(:,:,:,:)
337 4140 : zeta(:,:,:,:) =pimd_param%zeta(:,:,:,:)
338 4140 : dzeta(:,:,:,:) =pimd_param%dzeta(:,:,:,:) !unuseful to initialize zeta_next
339 :
340 : !Masses and spring constants (according to pitransform)
341 5 : select case(pitransform)
342 : case(0)
343 10 : ABI_MALLOC(mass,(natom,1))
344 10 : ABI_MALLOC(spring,(natom,1))
345 : case(1,2)
346 40 : ABI_MALLOC(mass,(natom,trotter))
347 45 : ABI_MALLOC(spring,(natom,trotter))
348 : end select
349 45 : spring_prim(:)=quantummass(:)*dble(trotter)*kt*kt
350 :
351 15 : call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,pitransform,trotter)
352 :
353 : !Recommended value of Nose mass
354 15 : write(msg,'(2a,f9.2,3a)') ch10,&
355 15 : & ' Recommended value of Nose mass is',one/(dble(trotter)*kt),' (atomic units)',ch10,&
356 30 : & '(see Tuckerman et al, J. Chem. Phys. 104, 5579 (1996))' ! [[cite:Tuckerman1996]]
357 15 : call wrtout(std_out,msg,'COLL')
358 :
359 : !Compute cartesian coordinates
360 105 : do iimage=1,trotter
361 90 : call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage))
362 105 : call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
363 : end do
364 :
365 : !Determine if it is a restart or not
366 : !If this is a calculation from scratch,generate random distribution of velocities
367 15 : irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel)
368 15 : if (irestart==0) then
369 3 : call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
370 : end if
371 :
372 : !Compute temperature at t
373 15 : temperature1=pimd_temperature(mass,vel)
374 :
375 : !################## Images evolution #####################
376 :
377 : !Transform the coordinates and forces (according to pitransform)
378 15 : call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
379 15 : call pimd_force_transform(forces,1,natom,pitransform,trotter) !compute staging forces
380 15 : call pimd_forces(forces,natom,spring,pitransform,trotter,xcart)
381 15 : call pimd_nosehoover_forces(dzeta,forces,forces_pimd,mass,natom,nnos,trotter,vel)
382 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
383 15 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
384 :
385 : !Compute atomic positions at t+dt
386 15 : if (itimimage<=1) then
387 :
388 : ! === 1st time step: single Taylor algorithm
389 : ! Predict positions
390 : call pimd_predict_taylor(dtion,forces_pimd,mass,natom,trotter,&
391 3 : & vel,xcart,xcart_next)
392 :
393 : ! Compute new temperature
394 3 : temperature2=pimd_temperature(mass,vel)
395 :
396 : ! Propagate the thermostat variables
397 : call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,&
398 3 : & thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform)
399 :
400 831 : dzeta=(zeta_next-zeta)/dtion
401 :
402 : else
403 :
404 : ! === Other time steps: Verlet algorithm + SC cycle
405 : ! Predict positions
406 12 : call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
407 : call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
408 12 : & xcart,xcart_next,xcart_prev)
409 : ! Propagate the thermostat variables
410 : call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,&
411 12 : & thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform)
412 : ! Self-consistent loop
413 12 : temperature2=pimd_temperature(mass,vel)
414 12 : temp2_prev=temperature2; tol=one
415 129 : do while (tol>tolerance)
416 : ! Recompute a (better) estimation of the velocity at time step t
417 6435 : vel = (xcart_next - xcart_prev) / (two*dtion)
418 32409 : dzeta=(zeta_next - zeta_prev) / (two*dtion)
419 117 : temperature2=pimd_temperature(mass,vel)
420 : ! Reestimate the force
421 117 : call pimd_nosehoover_forces(dzeta,forces,forces_pimd,mass,natom,nnos,trotter,vel)
422 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
423 117 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
424 : ! Compute new positions
425 : call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
426 117 : & xcart,xcart_next,xcart_prev)
427 : ! Propagate the thermostat variables
428 : call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,&
429 117 : & thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform)
430 : ! Compute variation of temperature (to check convergence of SC loop)
431 117 : tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev)
432 117 : temp2_prev=temperature2
433 : end do ! End self-consistent loop
434 :
435 : end if ! itimimage==1
436 :
437 : !Come back to primitive coordinates and velocities
438 15 : call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
439 15 : call pimd_coord_transform(xcart ,-1,natom,pitransform,trotter)
440 15 : call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
441 :
442 : !Compute contributions to energy
443 15 : call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring_prim,trotter,xcart)
444 :
445 : !Compute stress tensor at t from virial theorem
446 15 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
447 600 : stress_pimd=-stress_pimd ! Translate pressure to stress
448 :
449 : !############# Final operations ############################
450 :
451 : !Print messages
452 15 : vel_cell=zero;prtstress=1;if (prtvolimg>=2) prtstress=0
453 : call pimd_print(pimd_param%constraint,constraint_output,&
454 : & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
455 : & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
456 : & stress_pimd,temperature2,&
457 15 : & pimd_param%traj_unit,trotter,vel,vel_cell,xcart,xred)
458 :
459 : !If possible, estimate the velocities at t+dt
460 : !if (itimimage>1) then
461 15 : call pimd_predict_vel(dtion,itimimage,natom,trotter,0,xcart,xcart_next,xcart_prev,vel_next)
462 : !end if
463 :
464 : !Come back to reduced coordinates
465 105 : do iimage=1,trotter
466 105 : call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage))
467 : end do
468 :
469 : !update thermostat variables
470 4155 : dzeta = (three*zeta_next - four*zeta + zeta_prev)/(two * dtion)
471 4155 : zeta_prev=zeta
472 4155 : zeta=zeta_next
473 4140 : pimd_param%zeta_prev(:,:,:,:)=zeta_prev(:,:,:,:)
474 4140 : pimd_param%zeta(:,:,:,:) =zeta(:,:,:,:)
475 4140 : pimd_param%dzeta(:,:,:,:) =dzeta(:,:,:,:)
476 :
477 : !Free memory
478 15 : ABI_FREE(xcart)
479 15 : ABI_FREE(xcart_prev)
480 15 : ABI_FREE(xcart_next)
481 15 : ABI_FREE(forces_orig)
482 15 : ABI_FREE(forces_pimd)
483 15 : ABI_FREE(inertmass)
484 15 : ABI_FREE(quantummass)
485 15 : ABI_FREE(mass)
486 15 : ABI_FREE(spring)
487 15 : ABI_FREE(qmass)
488 15 : ABI_FREE(zeta_prev)
489 15 : ABI_FREE(zeta)
490 15 : ABI_FREE(zeta_next)
491 15 : ABI_FREE(dzeta)
492 :
493 15 : end subroutine pimd_nosehoover_nvt
494 : !!***
495 :
496 : end module m_pimd_nosehoover
497 : !!***
|
