TABLE OF CONTENTS


ABINIT/m_pimd_langevin [ Modules ]

[ Top ] [ Modules ]

NAME

 m_pimd_langevin

FUNCTION

COPYRIGHT

  Copyright (C) 2011-2022 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_pimd_langevin
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_pimd
28  use m_random_zbq
29 
30  use m_symtk,     only : matr3inv
31  use m_geometry,  only : xcart2xred, xred2xcart
32 
33 
34  implicit none
35 
36  private

ABINIT/pimd_langevin_npt [ Functions ]

[ Top ] [ Functions ]

NAME

 pimd_langevin_npt

FUNCTION

 Predicts new positions in Path Integral Molecular Dynamics using Langevin thermostat 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

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=volume 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 + 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
  vel_cell(3,3)=time derivative of cell parameters
    at input,  values at time t
    at output, values at time t+dt

NOTES

  Here follows PIMD in the NPT ensemble within the Langevin barostat algorithm
  of Quigley and Probert: J. Chem. Phys. 120, 11432 (2004) [[cite:Quigley2004]]
  and Comput. Phys. Comm. 169, 322 (2005) [[cite:Quigley2005]]

SOURCE

 94 subroutine pimd_langevin_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
 95 &          rprimd,rprimd_next,rprimd_prev,stressin,trotter,vel,vel_cell,&
 96 &          volume,xred,xred_next,xred_prev)
 97 
 98  implicit none
 99 
100 !Arguments ------------------------------------
101 !scalars
102  integer,intent(in) :: itimimage,natom,prtvolimg,trotter
103  real(dp),intent(in) :: volume
104  type(pimd_type),intent(in) :: pimd_param
105 !arrays
106  real(dp),intent(in) :: etotal(trotter),rprimd(3,3),rprimd_prev(3,3),stressin(3,3,trotter)
107  real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
108  real(dp),intent(out) :: rprimd_next(3,3),xred_next(3,natom,trotter)
109  real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter),vel_cell(3,3)
110 
111 !Local variables-------------------------------
112 !Options
113 !        Option for the Langevin algorithm correction
114  integer,parameter :: ilangevin=0
115 !        The following option forces the total of forces to be zero
116 !         It prevents the translation of the center of mass
117 !         If it is zero, no constraint on mass center is applied
118  integer,parameter :: zeroforce=1
119 !        Tolerance for the SC cycle
120  real(dp),parameter :: tolerance=tol9
121 
122 !scalars
123  integer :: idum=-5
124  integer :: constraint,iatom,ii,iimage,irestart,jj,ndof,prtstress
125  real(dp) :: dtion,eharm,eharm2,epot,friction,frictionbar,initemp,kt,rescale_temp,scalebar
126  real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol,tracepg,wg
127 !arrays
128  real, parameter :: identity(3,3)=reshape((/(one,(zero,ii=1,3),jj=1,2),one/),(/3,3/))
129  real(dp) :: aleabar(3,3),constraint_output(2),ddh(3,3),diffstress(3,3)
130  real(dp) :: dstrhh(3,3),fg(3,3),invrprimd(3,3)
131  real(dp) :: langev_bar(3,3),pg(3,3),pgdh(3,3),stress_pimd(3,3,3),strtarget(6),tmp(3,3)
132  real(dp),allocatable :: alea(:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:),forces_pimd_red(:,:)
133  real(dp),allocatable :: fsup(:,:),hxredpoint(:,:,:),inertmass(:),langev(:,:)
134  real(dp),allocatable ::  mass(:,:),quantummass(:),spring(:,:)
135  real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:)
136  real(dp),allocatable :: xredpoint(:,:,:)
137 
138 ! *************************************************************************
139 
140 !############# Initializations ###########################
141 
142 !Allocation of local arrays
143  ABI_MALLOC(xcart,(3,natom,trotter))
144  ABI_MALLOC(xcart_prev,(3,natom,trotter))
145  ABI_MALLOC(xcart_next,(3,natom,trotter))
146  ABI_MALLOC(forces_orig,(3,natom,trotter))
147  ABI_MALLOC(forces_pimd,(3,natom,trotter))
148  ABI_MALLOC(inertmass,(natom))
149  ABI_MALLOC(quantummass,(natom))
150 
151 !Fill in the local variables
152  ndof=3*natom*trotter
153  rescale_temp=one;if(zeroforce==1)rescale_temp=dble(ndof)/dble(ndof-3)
154  quantummass(1:natom)=pimd_param%amu   (pimd_param%typat(1:natom))*amu_emass
155  inertmass  (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
156  initemp=pimd_param%mdtemp(1)/rescale_temp; thermtemp=pimd_param%mdtemp(2)
157  dtion=pimd_param%dtion
158  kt=thermtemp*kb_HaK
159  friction=pimd_param%vis
160  wg=pimd_param%bmass
161  strtarget(:)=pimd_param%strtarget(:) ! imposed stress tensor
162  frictionbar=pimd_param%friction      ! friction coeff of barostat
163  scalebar=sqrt(two*frictionbar*wg*kt/dtion)
164  forces_orig=forces
165  constraint=0
166 
167 !Masses and spring constants
168  ABI_MALLOC(mass,(natom,1))
169  ABI_MALLOC(spring,(natom,1))
170  call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,0,trotter)
171 
172 !Initialize random forces
173  ABI_MALLOC(alea,(3,natom,trotter))
174  ABI_MALLOC(langev,(natom,trotter))
175  langev(:,1)=sqrt(two*friction*inertmass(:)*kt/dtion)
176  if(ilangevin==1)then
177    langev(:,1)=langev(:,1)*sqrt(one-(friction*dtion/(two*inertmass(:))))
178  end if
179 
180 !Random number generator initialization
181  if(itimimage<=1) then
182    call pimd_langevin_random_init(pimd_param%irandom,idum)
183  end if
184 
185 !Compute cartesian coordinates
186  do iimage=1,trotter
187    call xred2xcart(natom,rprimd,xcart     (:,:,iimage),xred(:,:,iimage))
188    call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
189  end do
190 
191 !Determine if it is a restart or not
192 !If this is a calculation from scratch,generate random distribution of velocities
193  irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel,vel_cell)
194 
195 !Initialize derivatives
196  if (mod(irestart,10)==0) then
197    call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
198  end if
199 !vel_cell does not depend on Trotter...
200  ddh(:,:)=vel_cell(:,:);if (irestart<10) ddh=zero
201 
202  if (itimimage<=1) then
203 
204 !  ========================= FIRST TIME STEP =======================================
205 
206    ABI_MALLOC(hxredpoint,(3,natom,trotter))
207    ABI_MALLOC(xredpoint,(3,natom,trotter))
208    ABI_MALLOC(forces_pimd_red,(3,natom))
209    ABI_MALLOC(fsup,(3,natom))
210 
211    do iimage=1,trotter
212      hxredpoint(:,:,iimage)=vel(:,:,iimage) - matmul(ddh(:,:),xred(:,:,iimage))
213    end do
214    call matr3inv(rprimd,invrprimd)
215    do iimage=1,trotter
216      xredpoint(:,:,iimage)=matmul(invrprimd(:,:),hxredpoint(:,:,iimage))
217    end do
218 
219 !  Generate random numbers
220    call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce)
221 
222 !  Compute PIMD and Langevin contributions to forces
223    call pimd_forces(forces,natom,spring,0,trotter,xcart)
224    call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint)
225    call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
226 &   mass,natom,trotter,pimd_param%wtatcon,xcart)
227    tmp=matmul(invrprimd,ddh)
228    pg=wg*matmul(ddh,invrprimd)
229    tracepg=pg(1,1)+pg(2,2)+pg(3,3)
230 
231 !  Taylor algorithm
232    do iimage=1,trotter
233      call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red)
234      fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage))
235      do iatom=1,natom
236        xred_next(:,iatom,iimage)=xred(:,iatom,iimage)+dtion*xredpoint(:,iatom,iimage) + &
237 &       half*(  &
238 &       forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom)- &
239 &       (tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) &
240 &       )*dtion*dtion/inertmass(iatom)
241      end do
242    end do
243 
244 !  predict rprimd at time t+dt from taylor algorithm
245    call pimd_langevin_random_bar(aleabar,pimd_param%irandom,idum)
246    langev_bar=matmul(aleabar,rprimd)*scalebar
247    call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart)
248    diffstress=pimd_diff_stress(stress_pimd,strtarget)
249    dstrhh=matmul(diffstress,rprimd)
250    pgdh=matmul(pg,ddh)
251    temperature1=pimd_temperature(mass,hxredpoint)*rescale_temp
252    fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature1*kb_HaK*rprimd(:,:)
253 
254    rprimd_next(:,:)=rprimd(:,:) + dtion*ddh(:,:) + half*(  &
255 &   fg(:,:)-wg*frictionbar*ddh(:,:)+langev_bar  &
256 &   )*dtion*dtion/wg
257 
258 !  Recompute xcart_next
259    do iimage=1,trotter
260      call xred2xcart(natom,rprimd_next,xcart_next(:,:,iimage),xred_next(:,:,iimage))
261    end do
262 
263 !  Compute stress tensor at t from virial theorem
264    call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
265 
266 !  Translate from pressure to stress for print
267    stress_pimd=-stress_pimd
268 
269 !  Compute temperature at current step
270    temperature1=pimd_temperature(mass,vel)*rescale_temp
271 
272 !  Estimate the velocities at t+dt
273    do iimage=1,trotter
274      do iatom=1,natom
275        xredpoint(:,iatom,iimage)=xredpoint(:,iatom,iimage)+(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom)- &
276             (tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))))*dtion/inertmass(iatom)
277      end do
278    end do
279    ddh(:,:)=ddh(:,:)+( fg(:,:)-wg*frictionbar*ddh(:,:)+langev_bar(:,:) )*dtion/wg
280    do iimage=1,trotter
281      call xred2xcart(natom,rprimd_next,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
282    end do
283    do iimage=1,trotter
284      do iatom=1,natom
285        vel(:,iatom,iimage)=hxredpoint(:,iatom,iimage)+matmul(ddh(:,:),xred(:,iatom,iimage))
286      end do
287    end do
288 
289 !  Compute new temperature
290    temperature2=pimd_temperature(mass,vel)*rescale_temp
291 
292    vel=xredpoint !scaled velocities transmitted to step 2
293 
294    ABI_FREE(hxredpoint)
295    ABI_FREE(xredpoint)
296    ABI_FREE(forces_pimd_red)
297    ABI_FREE(fsup)
298 
299  else
300 
301 !  ========================= OTHER TIME STEPS ======================================
302 
303 !  Additional allocations
304    ABI_MALLOC(hxredpoint,(3,natom,trotter))
305    ABI_MALLOC(xredpoint,(3,natom,trotter))
306    ABI_MALLOC(forces_pimd_red,(3,natom))
307    ABI_MALLOC(fsup,(3,natom))
308 
309    ddh=vel_cell(:,:)
310 
311    do iimage=1,trotter
312      hxredpoint(:,:,iimage)=vel(:,:,iimage) - matmul(ddh(:,:),xred(:,:,iimage))
313    end do
314 
315 !  first estimation of ddh, pg and its trace:
316    call matr3inv(rprimd,invrprimd)
317    pg=wg*matmul(ddh,invrprimd)
318    tracepg=pg(1,1)+pg(2,2)+pg(3,3)
319 
320 !  Momenta hxredpoint = H ds/dt: estimation
321    do iimage=1,trotter
322      hxredpoint(:,:,iimage)=matmul(rprimd,vel(:,:,iimage))
323    end do
324 
325 !  Compute temperature at t
326    temperature1=pimd_temperature(mass,hxredpoint)*rescale_temp
327 
328 !  Generate random numbers
329    call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce)
330 
331 !  Generate random numbers for the barostat
332    call pimd_langevin_random_bar(aleabar,pimd_param%irandom,idum)
333    langev_bar=matmul(aleabar,rprimd)*scalebar
334 
335 !  Compute PIMD and Langevin contributions to forces
336    call pimd_forces(forces,natom,spring,0,trotter,xcart)
337    call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint)
338    call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
339 &   mass,natom,trotter,pimd_param%wtatcon,xcart)
340 
341 !  Compute difference between instantaneous stress and imposed stress (barostat)
342    call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart)
343 
344    diffstress=pimd_diff_stress(stress_pimd,strtarget)
345 
346 !  Compute "force" on supercell vectors
347    dstrhh=matmul(diffstress,rprimd)
348    pgdh=matmul(pg,ddh)
349    fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature1*kb_HaK*rprimd(:,:)
350 
351 !  Evolve the supercell (fist estimation)
352    rprimd_next=two*rprimd-rprimd_prev+(fg-wg*frictionbar*ddh+langev_bar)*dtion*dtion/wg
353 
354 !  Evolve atomic positions (first estimation)
355    tmp=matmul(invrprimd,ddh)
356    do iimage=1,trotter
357      call xcart2xred(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
358      fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage))
359      call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red)
360 
361      do iatom=1,natom
362        xred_next(:,iatom,iimage)= &
363 &       two*xred(:,iatom,iimage) - xred_prev(:,iatom,iimage) &
364 &       +(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom) &
365 &       -tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) &
366 &       *dtion*dtion/inertmass(iatom)
367      end do
368    end do
369 
370 !  Self-consistent loop
371    temperature2=pimd_temperature(mass,xredpoint)*rescale_temp
372    temp2_prev=temperature2; tol=one
373 
374    do while (tol>tolerance)
375 !    Reestimate dH/dt at t
376      ddh(:,:)=(rprimd_next(:,:)-rprimd_prev(:,:))/(two*dtion)
377 
378 !    Reestimate the scaled velocities at t
379      do iimage=1,trotter
380        xredpoint(:,:,iimage)=(xred_next(:,:,iimage)-xred_prev(:,:,iimage))/(two*dtion)
381        call xred2xcart(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
382      end do
383 !    Reestimate the forces
384      call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint)
385      call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
386 &     mass,natom,trotter,pimd_param%wtatcon,xcart)
387 !    Compute variation of temperature (to check convergence of SC loop)
388      temperature2=pimd_temperature(mass,xredpoint)*rescale_temp
389      tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev)
390      temp2_prev=temperature2
391 !    Recompute the temperature
392      temperature2=pimd_temperature(mass,hxredpoint)*rescale_temp
393 !    Recompute pg
394      pg=wg*matmul(ddh,invrprimd)
395      tracepg=pg(1,1)+pg(2,2)+pg(3,3)
396 !    Recompute difference between instantaneous stress and imposed stress (barostat)
397      call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart)
398 
399      diffstress=pimd_diff_stress(stress_pimd,strtarget)
400 
401 !    Recompute "force" on supercell vectors
402      dstrhh=matmul(diffstress,rprimd)
403      pgdh=matmul(pg,ddh)
404      fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature2*kb_HaK*rprimd(:,:)
405 !    Evolve the supercell (better estimation)
406      rprimd_next=two*rprimd-rprimd_prev+(fg-wg*frictionbar*ddh+langev_bar)*dtion*dtion/wg
407 
408 !    Evolve atomic positions (better estimation):
409      tmp=matmul(invrprimd,ddh)
410      do iimage=1,trotter
411        call xcart2xred(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
412        fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage))
413        call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red)
414        do iatom=1,natom
415          xred_next(:,iatom,iimage)= &
416 &         two*xred(:,iatom,iimage) - xred_prev(:,iatom,iimage) &
417 &         +(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom) &
418 &         -(tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof)))) &
419 &         *dtion*dtion/inertmass(iatom)
420        end do
421      end do
422    end do ! End self-consistent loop
423 
424 !  Computation of true temperature from true velocities at t
425    do iimage=1,trotter
426      vel(:,:,iimage)=hxredpoint(:,:,iimage)+matmul(ddh,xred(:,:,iimage))
427    end do
428    temperature2=pimd_temperature(mass,vel)*rescale_temp
429 
430 !  Computation of the real stress tensor at t
431    call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
432 
433 !  translate from pressure to stress
434    stress_pimd=-stress_pimd
435 
436 !  Deallocations (Verlet algo)
437    ABI_FREE(xredpoint)
438    ABI_FREE(hxredpoint)
439    ABI_FREE(forces_pimd_red)
440    ABI_FREE(fsup)
441 
442  end if ! itimimage==1
443 
444 !############# Final operations ############################
445 
446 !Compute contributions to energy
447  call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring,trotter,xcart)
448 
449 !Print messages
450  prtstress=1
451  call pimd_print(constraint,constraint_output,&
452 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
453 & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
454 & stress_pimd,temperature1,temperature2,&
455 & pimd_param%traj_unit,trotter,vel,ddh,xcart,xred)
456 
457  if (itimimage>1) then
458 !  Estimation of ds/dt and ddh at t+dt
459    vel = (three*xred_next   - four*xred   + xred_prev)/(two * dtion)
460    ddh = (three*rprimd_next - four*rprimd + rprimd_prev)/(two * dtion)
461  end if
462 
463 !Return cell velocities (does not depend on Trotter)
464  vel_cell(:,:)=ddh(:,:)
465 
466 !Free memory
467  ABI_FREE(xcart)
468  ABI_FREE(xcart_prev)
469  ABI_FREE(xcart_next)
470  ABI_FREE(forces_orig)
471  ABI_FREE(forces_pimd)
472  ABI_FREE(inertmass)
473  ABI_FREE(quantummass)
474  ABI_FREE(mass)
475  ABI_FREE(spring)
476  ABI_FREE(alea)
477  ABI_FREE(langev)
478 
479 end subroutine pimd_langevin_npt

ABINIT/pimd_langevin_nvt [ Functions ]

[ Top ] [ Functions ]

NAME

 pimd_langevin_nvt

FUNCTION

 Predicts new positions in Path Integral Molecular Dynamics using Langevin thermostat 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

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=volume 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

   See Quigley,Probert, JCP 120, 11432 (2004) [[cite:Quigley2004]], part III

SOURCE

523 subroutine pimd_langevin_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
524 &                            rprimd,stressin,trotter,vel,volume,xred,xred_next,xred_prev)
525 
526  implicit none
527 
528 !Arguments ------------------------------------
529 !scalars
530  integer,intent(in) :: itimimage,natom,prtvolimg,trotter
531  real(dp),intent(in) :: volume
532  type(pimd_type),intent(in) :: pimd_param
533 !arrays
534  real(dp),intent(in) :: etotal(trotter),rprimd(3,3),stressin(3,3,trotter)
535  real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
536  real(dp),intent(out) :: xred_next(3,natom,trotter)
537  real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter)
538 
539 !Local variables-------------------------------
540 !Options
541 !        Option for the Langevin algorithm correction
542  integer,parameter :: ilangevin=0
543 !        Tolerance for the SC cycle
544  real(dp),parameter :: tolerance=tol9
545 
546 !scalars
547  integer :: idum=-5
548  integer :: iimage,irestart,ndof,pitransform,prtstress,use_qtb,zeroforce
549  real(dp) :: dtion,eharm,eharm2,epot,friction,initemp,kt,kt_,rescale_temp
550  real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol
551 !arrays
552  real(dp) :: constraint_output(2),spring_prim(natom),stress_pimd(3,3,3),vel_cell(3,3)
553  real(dp),allocatable :: alea(:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:)
554  real(dp),allocatable :: inertmass(:),langev(:,:),mass(:,:),quantummass(:),spring(:,:)
555  real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:)
556 
557 ! *************************************************************************
558 
559  if (pimd_param%use_qtb==1.and.pimd_param%qtb_file_unit<=0) then
560    ABI_BUG('piqtb_force not open!')
561  end if
562 
563 !############# Initializations ###########################
564 
565  pitransform=pimd_param%pitransform
566 
567 !The following option forces the total of forces to be zero
568 !It prevents the translation of the center of mass
569 !If it is zero, no constraint on mass center is applied
570  zeroforce=1
571  if(pitransform==1) zeroforce=0
572  if(pitransform==2) zeroforce=0
573  if(pimd_param%constraint==1) zeroforce=0
574 
575 !Allocation of local arrays
576  ABI_MALLOC(xcart,(3,natom,trotter))
577  ABI_MALLOC(xcart_prev,(3,natom,trotter))
578  ABI_MALLOC(xcart_next,(3,natom,trotter))
579  ABI_MALLOC(forces_orig,(3,natom,trotter))
580  ABI_MALLOC(forces_pimd,(3,natom,trotter))
581  ABI_MALLOC(inertmass,(natom))
582  ABI_MALLOC(quantummass,(natom))
583 
584 !Fill in the local variables
585  use_qtb=pimd_param%use_qtb
586  ndof=3*natom*trotter
587  rescale_temp=one;if(zeroforce==1)rescale_temp=dble(ndof)/dble(ndof-3)
588  quantummass(1:natom)=pimd_param%amu   (pimd_param%typat(1:natom))*amu_emass
589  inertmass  (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
590  if(pitransform==1) inertmass=quantummass !compulsory for good definition of normal mode masses
591  if(pitransform==2) inertmass=quantummass !compulsory for good definition of staging masses
592  initemp=pimd_param%mdtemp(1)/rescale_temp;thermtemp=pimd_param%mdtemp(2)
593  friction=pimd_param%vis;dtion=pimd_param%dtion
594  kt=thermtemp*kb_HaK
595  forces_orig=forces
596 
597 !Masses and spring constants
598  select case(pitransform)
599  case(0)
600    ABI_MALLOC(mass,(natom,1))   ! This second dimension is needed
601    ABI_MALLOC(spring,(natom,1))
602    ABI_MALLOC(langev,(natom,1))
603  case(1,2)
604    ABI_MALLOC(mass,(natom,trotter))
605    ABI_MALLOC(spring,(natom,trotter))
606    ABI_MALLOC(langev,(natom,trotter))
607  end select
608  spring_prim(:)=quantummass(:)*dble(trotter)*kt*kt
609  call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,pitransform,trotter)
610 
611 !Initialize random forces
612  ABI_MALLOC(alea,(3,natom,trotter))
613  if (use_qtb==0) then
614    langev(:,:)=sqrt(two*friction*mass(:,:)*kt/dtion)
615  else
616    langev(:,:)=sqrt(two*friction*mass(:,:))
617  end if
618 
619 !Random number generator initialization
620  if(itimimage<=1) then
621    call pimd_langevin_random_init(pimd_param%irandom,idum)
622  end if
623 
624 !Compute cartesian coordinates
625  do iimage=1,trotter
626    call xred2xcart(natom,rprimd,xcart     (:,:,iimage),xred(:,:,iimage))
627    call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
628  end do
629 
630 !Determine if it is a restart or not
631 !If this is a calculation from scratch,generate random distribution of velocities
632  irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel)
633  if (irestart==0) then
634    call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
635  end if
636 
637 !Compute temperature at t
638  temperature1=pimd_temperature(mass,vel)*rescale_temp
639 
640 !################## Images evolution #####################
641 
642 !Generate random numbers
643  if (use_qtb==0) then
644    call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce)
645  else
646    call pimd_langevin_random_qtb(alea,langev,mass,natom,pimd_param%qtb_file_unit,trotter,zeroforce)
647  end if
648 
649 !Compute PIMD and Langevin contributions to forces
650  call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
651  call pimd_force_transform(forces,1,natom,pitransform,trotter) !compute staging forces
652  call pimd_forces(forces,natom,spring,pitransform,trotter,xcart)
653  call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,vel)
654  call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
655 & mass,natom,trotter,pimd_param%wtatcon,xcart)
656 
657 !Compute atomic positions at t+dt
658  if (itimimage<=1) then
659 
660 !  === 1st time step: single Taylor algorithm
661 !  Predict positions
662    call pimd_predict_taylor(dtion,forces_pimd,mass,natom,trotter,&
663 &   vel,xcart,xcart_next)
664 
665 !  Estimate the velocities at t+dt/2
666    vel=(xcart_next-xcart)/dtion
667 
668 !  Compute new temperature
669    temperature2=pimd_temperature(mass,vel)*rescale_temp
670 
671  else
672 
673 !  === Other time steps: Verlet algorithm + SC cycle
674 !  Predict positions
675    call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
676    call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
677 &   xcart,xcart_next,xcart_prev)
678 !  Self-consistent loop
679    temperature2=pimd_temperature(mass,vel)*rescale_temp
680    temp2_prev=temperature2; tol=one
681    do while (tol>tolerance)
682 !    Recompute a (better) estimation of the velocity at time step t
683      vel = (xcart_next - xcart_prev) / (two*dtion)
684      temperature2=pimd_temperature(mass,vel)*rescale_temp
685 !    Reestimate the force
686      call pimd_langevin_forces(alea,forces,forces_pimd,friction,&
687 &     langev,mass,natom,trotter,vel)
688      call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
689 &     mass,natom,trotter,pimd_param%wtatcon,xcart)
690 !    Compute new positions
691      call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
692 &     xcart,xcart_next,xcart_prev)
693 
694 !    Compute variation of temperature (to check convergence of SC loop)
695      tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev)
696      temp2_prev=temperature2
697 
698    end do ! End self-consistent loop
699 
700  end if ! itimimage==1
701 
702  call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
703  call pimd_coord_transform(xcart,-1,natom,pitransform,trotter)
704  call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
705 
706 !Compute contributions to energy
707  call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring_prim,trotter,xcart)
708 
709 !Compute stress tensor at t
710  if (use_qtb==0) then
711    call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
712  else
713    call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,temperature2,thermtemp,trotter,vel,volume,xcart)
714  end if
715  stress_pimd=-stress_pimd ! Translate pressure to stress
716 
717 !############# Final operations ############################
718 
719 !Print messages
720  vel_cell=zero;prtstress=1;if (prtvolimg>=2) prtstress=0
721  kt_=kt;if (use_qtb==1) kt_=temperature2*kb_HaK
722  call pimd_print(pimd_param%constraint,constraint_output,&
723 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
724 & itimimage,kt_,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
725 & stress_pimd,temperature1,temperature2,&
726 & pimd_param%traj_unit,trotter,vel,vel_cell,xcart,xred)
727 
728 !If possible, estimate the (transformed) velocities at t+dt
729  if (itimimage>1) then
730    call pimd_coord_transform(xcart_next,1,natom,pitransform,trotter)
731    call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
732    call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
733    vel = (three*xcart_next - four*xcart + xcart_prev)/(two * dtion)
734    call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
735    call pimd_coord_transform(xcart,-1,natom,pitransform,trotter)
736    call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
737  end if
738 
739 !Come back to reduced coordinates
740  do iimage=1,trotter
741    call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage))
742  end do
743 
744 !Free memory
745  ABI_FREE(xcart)
746  ABI_FREE(xcart_prev)
747  ABI_FREE(xcart_next)
748  ABI_FREE(forces_orig)
749  ABI_FREE(forces_pimd)
750  ABI_FREE(inertmass)
751  ABI_FREE(quantummass)
752  ABI_FREE(mass)
753  ABI_FREE(spring)
754  ABI_FREE(alea)
755  ABI_FREE(langev)
756 
757 end subroutine pimd_langevin_nvt