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