TABLE OF CONTENTS


ABINIT/m_pimd_nosehoover [ Modules ]

[ Top ] [ Modules ]

NAME

   m_pimd_nosehoover

FUNCTION

COPYRIGHT

  Copyright (C) 2011-2018 ABINIT group (GG,MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_pimd_nosehoover
28 
29  use defs_basis
30  use m_pimd
31  use m_abicore
32 
33  use m_geometry,  only : xcart2xred, xred2xcart
34 
35  implicit none
36 
37  private

ABINIT/pimd_nosehoover_npt [ Functions ]

[ Top ] [ Functions ]

NAME

 pimd_nosehoover_npt

FUNCTION

 Predicts new positions in Path Integral Molecular Dynamics using Nose-Hoover in the NPT ensemble.
 Given the positions at time t and t-dtion, an estimation of the velocities at time t,
 the forces and an estimation of the stress at time t, and an estimation of the cell at time t,
 computes in the Path Integral Molecular Dynamics framework the new positions at time t+dtion,
 computes self-consistently the velocities, the stress and the cell at time t and produces
 an estimation of the velocities, stress and new cell at time t+dtion
 No change of acell and rprim at present.

INPUTS

  etotal(trotter)=electronic total energy for all images
  itimimage=number of the current time for image propagation (itimimage+1 is to be predicted here)
  natom=dimension of vel_timimage and xred_timimage
  pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
  prtvolimg=printing volume
  rprimd(3,3)=dimensionless unit cell vectors (common to all images) at time t (present time step)
  rprimd_prev(3,3)=dimensionless unit cell vectors (common to all images) at time t-dt (previous time step)
  stressin(3,3,trotter)=electronic stress tensor for each image
  trotter=Trotter number (total number of images)
  volume=voume of unit cell (common to all images)
  xred(3,natom,trotter)=reduced coordinates of atoms for all images at time t (present time step)
  xred_prev(3,natom,trotter)=reduced coordinates of atoms for all images at time t-dt (previous time step)

OUTPUT

  rprimd_next(3,3)=dimensionless unit cell vectors (common to all images) at time t+dt (next time step)
  xred_next(3,natom,trotter)=reduced coordinates of atoms for all images at time t+dt (next time step)

SIDE EFFECTS

  forces(3,natom,trotter)=forces over atoms for all images
    at input,  electronic forces
    at output, electronic forces + Langevin contribution
  vel(3,natom,trotter)=velocies of atoms for all images
    at input,  values at time t
    at output, values at time t+dt
  vel_cell(3,3,trotter)=time derivative of cell parameters
    at input,  values at time t
    at output, values at time t+dt

PARENTS

      predict_pimd

CHILDREN

      pimd_energies,pimd_initvel,pimd_print,pimd_stresses,xcart2xred
      xred2xcart

SOURCE

 98 subroutine pimd_nosehoover_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
 99 &          rprimd,rprimd_next,rprimd_prev,stressin,trotter,vel,vel_cell,&
100 &          volume,xred,xred_next,xred_prev)
101 
102 
103 !This section has been created automatically by the script Abilint (TD).
104 !Do not modify the following lines by hand.
105 #undef ABI_FUNC
106 #define ABI_FUNC 'pimd_nosehoover_npt'
107 !End of the abilint section
108 
109  implicit none
110 
111 !Arguments ------------------------------------
112 !scalars
113  integer,intent(in) :: itimimage,natom,prtvolimg,trotter
114  real(dp),intent(in) :: volume
115  type(pimd_type),intent(in) :: pimd_param
116 !arrays
117  real(dp),intent(in) :: etotal(trotter),rprimd(3,3),rprimd_prev(3,3),stressin(3,3,trotter)
118  real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
119  real(dp),intent(out) :: rprimd_next(3,3),xred_next(3,natom,trotter)
120  real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter),vel_cell(3,3,trotter)
121 
122 !Local variables-------------------------------
123 !Options
124  real(dp),parameter :: tolerance=tol7 ! SCF tolerance
125 !scalars
126  integer :: idum=-5
127  integer :: constraint,iimage,irestart,ndof,nnos,pitransform,prtstress
128  real(dp) :: dtion,eharm,eharm2,epot,initemp,kt,temperature1,temperature2,thermtemp
129 !arrays
130  real(dp) :: constraint_output(2),ddh(3,3),stress_pimd(3,3,3)
131  real(dp),allocatable :: dzeta(:,:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:)
132  real(dp),allocatable :: inertmass(:),masseff(:,:),qmass(:),quantummass(:),springeff(:,:)
133  real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:),zeta(:)
134 
135 ! *************************************************************************
136 
137 !############# Initializations ###########################
138 
139 !Allocation of local arrays
140  ABI_ALLOCATE(xcart,(3,natom,trotter))
141  ABI_ALLOCATE(xcart_prev,(3,natom,trotter))
142  ABI_ALLOCATE(xcart_next,(3,natom,trotter))
143  ABI_ALLOCATE(forces_orig,(3,natom,trotter))
144  ABI_ALLOCATE(forces_pimd,(3,natom,trotter))
145  ABI_ALLOCATE(inertmass,(natom))
146  ABI_ALLOCATE(quantummass,(natom))
147 
148 !Fill in the local variables
149  ndof=3*natom*trotter
150  quantummass(1:natom)=pimd_param%amu   (pimd_param%typat(1:natom))*amu_emass
151  inertmass  (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
152  initemp=pimd_param%mdtemp(1);thermtemp=pimd_param%mdtemp(2)
153  dtion=pimd_param%dtion;pitransform=pimd_param%pitransform
154  kt=thermtemp*kb_HaK
155  forces_orig=forces
156 
157 !Allocation/initialization of local variables used for Nose-Hoover chains
158 !Associated variables:
159 !nnos = number of thermostats
160 !dzeta(3,natom,trotter,nnos) = variables of thermostats, in (atomic time unit)^(-1)
161 !qmass(nnos) = masses of thermostats
162 !specific to PIMD: pitransform = coordinate transformation (0:no; 1:normal mode; 2:staging)
163  nnos=pimd_param%nnos
164  ABI_ALLOCATE(qmass,(nnos))
165  ABI_ALLOCATE(zeta,(nnos))
166  ABI_ALLOCATE(dzeta,(3,natom,trotter,nnos))
167  qmass(1:nnos)=pimd_param%qmass(1:nnos)
168  zeta=zero;dzeta=zero
169 
170 !Compute cartesian coordinates
171  do iimage=1,trotter
172    call xred2xcart(natom,rprimd,xcart     (:,:,iimage),xred(:,:,iimage))
173    call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
174  end do
175 
176 !Determine if it is a restart or not
177 !If this is a calculation from scratch,generate random distribution of velocities
178  irestart=1;if (itimimage==1) irestart=pimd_is_restart(masseff,vel,vel_cell)
179 
180 !Initialize derivatives
181  if (mod(irestart,10)==0) then
182    call pimd_initvel(idum,masseff,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
183  end if
184 !vel_cell does not depend on Trotter...
185  ddh=vel_cell(:,:,1);if (irestart<10) ddh=zero
186 
187 !Compute temperature at t
188  temperature1=pimd_temperature(masseff,vel)
189 
190 !################## Images evolution #####################
191 
192 !This is temporary
193  xcart_next=zero
194  rprimd_next=rprimd_prev
195  temperature2=pimd_temperature(masseff,vel)
196 
197 !Compute contributions to energy
198  call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,springeff,trotter,xcart)
199 
200 !Compute stress tensor at t from virial theorem
201  call pimd_stresses(masseff,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
202 
203 
204 !############# Final operations ############################
205 
206 !Print messages
207  prtstress=1
208  call pimd_print(constraint,constraint_output,&
209 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
210 & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
211 & stress_pimd,temperature1,temperature2,&
212 & pimd_param%traj_unit,trotter,vel,ddh,xcart,xred)
213 
214 !Come back to reduced coordinates
215  do iimage=1,trotter
216    call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage))
217  end do
218 
219 !Return cell velocities (does not depend on Trotter)
220  do iimage=1,trotter
221    vel_cell(:,:,iimage)=ddh(:,:)
222  end do
223 
224 !Free memory
225  ABI_DEALLOCATE(xcart)
226  ABI_DEALLOCATE(xcart_prev)
227  ABI_DEALLOCATE(xcart_next)
228  ABI_DEALLOCATE(forces_orig)
229  ABI_DEALLOCATE(forces_pimd)
230  ABI_DEALLOCATE(inertmass)
231  ABI_DEALLOCATE(quantummass)
232  ABI_DEALLOCATE(masseff)
233  ABI_DEALLOCATE(springeff)
234  ABI_DEALLOCATE(qmass)
235  ABI_DEALLOCATE(dzeta)
236  ABI_DEALLOCATE(zeta)
237 
238 end subroutine pimd_nosehoover_npt

ABINIT/pimd_nosehoover_nvt [ Functions ]

[ Top ] [ Functions ]

NAME

 pimd_nosehoover_nvt

FUNCTION

 Predicts new positions in Path Integral Molecular Dynamics using Nose-Hoover in the NVT ensemble.
 Given the positions at time t and t-dtion, an estimation of the velocities at time t,
 the forces and an estimation of the stress at time t, and an estimation of the cell at time t,
 computes in the Path Integral Molecular Dynamics framework the new positions at time t+dtion,
 computes self-consistently the velocities, the stress and the cell at time t and produces
 an estimation of the velocities, stress and new cell at time t+dtion
 No change of acell and rprim at present.

INPUTS

  etotal(trotter)=electronic total energy for all images
  itimimage=number of the current time for image propagation (itimimage+1 is to be predicted here)
  natom=dimension of vel_timimage and xred_timimage
  pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
  prtvolimg=printing volume
  rprimd(3,3)=dimensionless unit cell vectors (common to all images)
  stressin(3,3,trotter)=electronic stress tensor for each image
  trotter=Trotter number (total number of images)
  volume=voume of unit cell (common to all images)
  xred(3,natom,trotter)=reduced coordinates of atoms for all images at time t (present time step)
  xred_prev(3,natom,trotter)=reduced coordinates of atoms for all images at time t-dt (previous time step)

OUTPUT

  xred_next(3,natom,trotter)=reduced coordinates of atoms for all images at time t+dt (next time step)

SIDE EFFECTS

  forces(3,natom,trotter)=forces over atoms for all images
    at input,  electronic forces
    at output, electronic forces + quantum spring contribution
  vel(3,natom,trotter)=velocies of atoms for all images
    at input,  values at time t
    at output, values at time t+dt

NOTES

  Thermization by Nose-Hoover chains according to
  Martyna, Klein, Tuckerman, J. Chem. Phys. 97, 2635 (1992) [[cite:Martyna1992]]
  Tuckerman, Marx, Klein, Parrinello, J. Chem. Phys. 104, 5579 (1996) [[cite:Tuckerman1996]]

PARENTS

      predict_pimd

CHILDREN

      pimd_apply_constraint,pimd_coord_transform,pimd_energies
      pimd_force_transform,pimd_forces,pimd_initvel,pimd_mass_spring
      pimd_nosehoover_forces,pimd_nosehoover_propagate,pimd_predict_taylor
      pimd_predict_verlet,pimd_print,pimd_stresses,wrtout,xcart2xred
      xred2xcart

SOURCE

295 subroutine pimd_nosehoover_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
296 &                              rprimd,stressin,trotter,vel,volume,xred,xred_next,xred_prev)
297 
298 
299 !This section has been created automatically by the script Abilint (TD).
300 !Do not modify the following lines by hand.
301 #undef ABI_FUNC
302 #define ABI_FUNC 'pimd_nosehoover_nvt'
303 !End of the abilint section
304 
305  implicit none
306 
307 !Arguments ------------------------------------
308 !scalars
309  integer,intent(in) :: itimimage,natom,prtvolimg,trotter
310  real(dp),intent(in) :: volume
311  type(pimd_type),intent(inout) :: pimd_param
312 !arrays
313  real(dp),intent(in) :: etotal(trotter),rprimd(3,3),stressin(3,3,trotter)
314  real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
315  real(dp),intent(out) :: xred_next(3,natom,trotter)
316  real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter)
317 
318 !Local variables-------------------------------
319 !Options
320  real(dp),parameter :: tolerance=tol9 ! SCF tolerance
321 !scalars
322  integer :: idum=-5
323  integer :: iimage,irestart,ndof,nnos,pitransform,prtstress
324  real(dp) :: dtion,eharm,eharm2,epot,initemp,kt
325  real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol
326  character(len=500) :: msg
327 !arrays
328  real(dp) :: constraint_output(2),spring_prim(natom),stress_pimd(3,3,3),vel_cell(3,3)
329  real(dp),allocatable :: forces_orig(:,:,:),forces_pimd(:,:,:)
330  real(dp),allocatable :: inertmass(:),mass(:,:),qmass(:),quantummass(:),spring(:,:)
331  real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:)
332  real(dp),allocatable :: dzeta(:,:,:,:),zeta_prev(:,:,:,:),zeta(:,:,:,:)
333  real(dp),allocatable :: zeta_next(:,:,:,:)
334 
335 ! *************************************************************************
336 
337 !############# Initializations ###########################
338 
339 !Allocation of local arrays
340  ABI_ALLOCATE(xcart,(3,natom,trotter))
341  ABI_ALLOCATE(xcart_prev,(3,natom,trotter))
342  ABI_ALLOCATE(xcart_next,(3,natom,trotter))
343  ABI_ALLOCATE(forces_orig,(3,natom,trotter))
344  ABI_ALLOCATE(forces_pimd,(3,natom,trotter))
345  ABI_ALLOCATE(inertmass,(natom))
346  ABI_ALLOCATE(quantummass,(natom))
347 
348 !Fill in the local variables
349  ndof=3*natom*trotter
350  quantummass(1:natom)=pimd_param%amu   (pimd_param%typat(1:natom))*amu_emass
351  inertmass  (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
352  if(pitransform==1) inertmass=quantummass !compulsory for good definition of normal mode masses
353  if(pitransform==2) inertmass=quantummass !compulsory for good definition of staging masses
354  initemp=pimd_param%mdtemp(1);thermtemp=pimd_param%mdtemp(2)
355  dtion=pimd_param%dtion;pitransform=pimd_param%pitransform
356  kt=thermtemp*kb_HaK
357  forces_orig=forces
358 
359 !Allocation/initialization of local variables used for Nose-Hoover chains
360 !Associated variables:
361 !nnos = number of thermostats
362 !zeta,zeta_next,zeta_prev(3,natom,trotter,nnos) = variables of thermostats, dzeta, its time derivative
363 !qmass(nnos) = masses of thermostats
364 !specific to PIMD: pitransform = coordinate transformation (0:no; 1:normal mode; 2:staging)
365  nnos=pimd_param%nnos
366  ABI_ALLOCATE(qmass,(nnos))
367  ABI_ALLOCATE(zeta_prev,(3,natom,trotter,nnos))
368  ABI_ALLOCATE(zeta,(3,natom,trotter,nnos))
369  ABI_ALLOCATE(zeta_next,(3,natom,trotter,nnos))
370  ABI_ALLOCATE(dzeta,(3,natom,trotter,nnos))
371 !initialization
372  qmass(1:nnos)=pimd_param%qmass(1:nnos)
373  zeta_prev(:,:,:,:)=pimd_param%zeta_prev(:,:,:,:)
374  zeta(:,:,:,:)     =pimd_param%zeta(:,:,:,:)
375  dzeta(:,:,:,:)    =pimd_param%dzeta(:,:,:,:) !unuseful to initialize zeta_next
376 
377 !Masses and spring constants (according to pitransform)
378  select case(pitransform)
379  case(0)
380    ABI_ALLOCATE(mass,(natom,1))
381    ABI_ALLOCATE(spring,(natom,1))
382  case(1,2)
383    ABI_ALLOCATE(mass,(natom,trotter))
384    ABI_ALLOCATE(spring,(natom,trotter))
385  end select
386  spring_prim(:)=quantummass(:)*dble(trotter)*kt*kt
387 
388  call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,pitransform,trotter)
389 
390 !Recommended value of Nose mass
391  write(msg,'(2a,f9.2,3a)') ch10,&
392 & ' Recommended value of Nose mass is',one/(dble(trotter)*kt),' (atomic units)',ch10,&
393 & '(see Tuckerman et al, J. Chem. Phys. 104, 5579 (1996))' ! [[cite:Tuckerman1996]]
394  call wrtout(std_out,msg,'COLL')
395 
396 !Compute cartesian coordinates
397  do iimage=1,trotter
398    call xred2xcart(natom,rprimd,xcart     (:,:,iimage),xred(:,:,iimage))
399    call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
400  end do
401 
402 !Determine if it is a restart or not
403 !If this is a calculation from scratch,generate random distribution of velocities
404  irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel)
405  if (irestart==0) then
406    call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
407  end if
408 
409 !Compute temperature at t
410  temperature1=pimd_temperature(mass,vel)
411 
412 !################## Images evolution #####################
413 
414 !Transform the coordinates and forces (according to pitransform)
415  call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
416  call pimd_force_transform(forces,1,natom,pitransform,trotter) !compute staging forces
417  call pimd_forces(forces,natom,spring,pitransform,trotter,xcart)
418  call pimd_nosehoover_forces(dzeta,forces,forces_pimd,mass,natom,nnos,trotter,vel)
419  call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
420 & mass,natom,trotter,pimd_param%wtatcon,xcart)
421 
422 !Compute atomic positions at t+dt
423  if (itimimage<=1) then
424 
425 !  === 1st time step: single Taylor algorithm
426 !  Predict positions
427    call pimd_predict_taylor(dtion,forces_pimd,mass,natom,trotter,&
428 &   vel,xcart,xcart_next)
429 
430 !  Estimate the velocities at t+dt/2
431    vel=(xcart_next-xcart)/dtion
432 
433 !  Compute new temperature
434    temperature2=pimd_temperature(mass,vel)
435 
436 !  Propagate the thermostat variables
437    call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,&
438 &   thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform)
439 
440    dzeta=(zeta_next-zeta)/dtion
441 
442  else
443 
444 !  === Other time steps: Verlet algorithm + SC cycle
445 !  Predict positions
446    call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
447    call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
448 &   xcart,xcart_next,xcart_prev)
449 !  Propagate the thermostat variables
450    call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,&
451 &   thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform)
452 !  Self-consistent loop
453    temperature2=pimd_temperature(mass,vel)
454    temp2_prev=temperature2; tol=one
455    do while (tol>tolerance)
456 !    Recompute a (better) estimation of the velocity at time step t
457      vel = (xcart_next - xcart_prev) / (two*dtion)
458      dzeta=(zeta_next  - zeta_prev)  / (two*dtion)
459      temperature2=pimd_temperature(mass,vel)
460 !    Reestimate the force
461      call pimd_nosehoover_forces(dzeta,forces,forces_pimd,mass,natom,nnos,trotter,vel)
462      call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
463 &     mass,natom,trotter,pimd_param%wtatcon,xcart)
464 !    Compute new positions
465      call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
466 &     xcart,xcart_next,xcart_prev)
467 !    Propagate the thermostat variables
468      call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,&
469 &     thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform)
470 !    Compute variation of temperature (to check convergence of SC loop)
471      tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev)
472      temp2_prev=temperature2
473    end do ! End self-consistent loop
474 
475  end if ! itimimage==1
476 
477 !Come back to primitive coordinates and velocities
478  call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
479  call pimd_coord_transform(xcart     ,-1,natom,pitransform,trotter)
480  call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
481 
482 !Compute contributions to energy
483  call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring_prim,trotter,xcart)
484 
485 !Compute stress tensor at t from virial theorem
486  call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
487  stress_pimd=-stress_pimd ! Translate pressure to stress
488 
489 !############# Final operations ############################
490 
491 !Print messages
492  vel_cell=zero;prtstress=1;if (prtvolimg>=2) prtstress=0
493  call pimd_print(pimd_param%constraint,constraint_output,&
494 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
495 & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
496 & stress_pimd,temperature1,temperature2,&
497 & pimd_param%traj_unit,trotter,vel,vel_cell,xcart,xred)
498 
499 !If possible, estimate the (transformed) velocities at t+dt
500  if (itimimage>1) then
501    call pimd_coord_transform(xcart_next,1,natom,pitransform,trotter)
502    call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
503    call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
504    vel = (three*xcart_next - four*xcart + xcart_prev)/(two * dtion)
505    call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
506    call pimd_coord_transform(xcart,-1,natom,pitransform,trotter)
507    call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
508  end if
509 
510 !Come back to reduced coordinates
511  do iimage=1,trotter
512    call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage))
513  end do
514 
515 !update thermostat variables
516  dzeta = (three*zeta_next - four*zeta + zeta_prev)/(two * dtion)
517  zeta_prev=zeta
518  zeta=zeta_next
519  pimd_param%zeta_prev(:,:,:,:)=zeta_prev(:,:,:,:)
520  pimd_param%zeta(:,:,:,:)     =zeta(:,:,:,:)
521  pimd_param%dzeta(:,:,:,:)    =dzeta(:,:,:,:)
522 
523 !Free memory
524  ABI_DEALLOCATE(xcart)
525  ABI_DEALLOCATE(xcart_prev)
526  ABI_DEALLOCATE(xcart_next)
527  ABI_DEALLOCATE(forces_orig)
528  ABI_DEALLOCATE(forces_pimd)
529  ABI_DEALLOCATE(inertmass)
530  ABI_DEALLOCATE(quantummass)
531  ABI_DEALLOCATE(mass)
532  ABI_DEALLOCATE(spring)
533  ABI_DEALLOCATE(qmass)
534  ABI_DEALLOCATE(zeta_prev)
535  ABI_DEALLOCATE(zeta)
536  ABI_DEALLOCATE(zeta_next)
537  ABI_DEALLOCATE(dzeta)
538 
539 end subroutine pimd_nosehoover_nvt