TABLE OF CONTENTS


ABINIT/m_pimd_langevin [ Modules ]

[ Top ] [ Modules ]

NAME

 m_pimd_langevin

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_langevin
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_pimd
33  use m_random_zbq
34 
35  use m_symtk,     only : matr3inv
36  use m_geometry,  only : xcart2xred, xred2xcart
37 
38 
39  implicit none
40 
41  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=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 + 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,trotter)=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]]

PARENTS

      predict_pimd

CHILDREN

      matr3inv,pimd_apply_constraint,pimd_energies,pimd_forces,pimd_initvel
      pimd_langevin_forces,pimd_langevin_random,pimd_langevin_random_bar
      pimd_langevin_random_init,pimd_mass_spring,pimd_print,pimd_stresses
      xcart2xred,xred2xcart

SOURCE

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

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

PARENTS

      predict_pimd

CHILDREN

      pimd_apply_constraint,pimd_coord_transform,pimd_energies
      pimd_force_transform,pimd_forces,pimd_initvel,pimd_langevin_forces
      pimd_langevin_random,pimd_langevin_random_init,pimd_langevin_random_qtb
      pimd_mass_spring,pimd_predict_taylor,pimd_predict_verlet,pimd_print
      pimd_stresses,xcart2xred,xred2xcart

SOURCE

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