@@ 30,10 +30,395 @@ Module **PIM_qtb** generates trajectories based on several classical and semic
These trajectories can be used to sample initial conditions for intramolecular vibrationalenergy redistribution (IVR) dynamics.
Applications of the module
__________________________
This module has been used to calculate time correlation function of an anharmonic oscillator parameterized on a hydroxyle bond OH and to compute pair correlation functions of a 13 atoms' neon cluster described by a LennardJones potential.
Description of the module
=========================
The module implements various methods based on Langevin dynamics to
sample initial conditions for IVR or to directly exploit the generated
trajectories. The methods implemented are: classical Langevin dynamics,
Quantum Thermal Bath (QTB) and two variants of adaptive QTB (adQTBr and
adQTBf).
Classical Langevin dynamics

Classical Langevin dynamics is described by a stochastic differential
equation :
.. math:: m \dot v = \nabla U m \gamma v + F(t)
:label: eqLGV
where :math:`\nabla U` is the vector forces exerted on the sets of
atoms, :math:`m` the mass of atoms, :math:`\gamma` the damping term
determining the strength of a viscous force and :math:`F(t)` a
stochastic noise. The stochastic noise is gaussian and deltacorrelated.
To ensure the classical fluctuationdissipation theorem, one might write
The parameter :math:`\gamma_{f,j}` are then modified to satisfy
:math:`\Delta_{FDT} = 0` (eq. :eq:`eqDFDT`). In this method, one should take care of checking results convergence by decreasing the :math:`\alpha` parameter.
Input file
==========
To run PaPIM using one of the Langevin methods, one must set the
parameter *sampling\_type* in the *sampling* section to one of the
following values:
 classical\_langevin
 qtb
 adqtbr
 adqtbf
 In this case the parameters *n\_equilibration\_steps* and
*n\_mc\_steps* are ignored and the section *langevin* is read.

 The section *langevin* must specify the following parameters:
 *dt* : time step of the Langevin dynamics (REAL)
 *lgv\_nsteps* : number of Langevin steps between each IVR sample
(INTEGER)
 *lgv\_nsteps\_therm* : number of thermalization steps (INTEGER)
 *integrator* : integration method (two splitting methods are
currently implemented: BAOAB, ABOBA (see reference
[Lei]_ )) (STRING,
default=“ABOBA”)
 *damping* : base damping coefficient for production runs
(:math:`\gamma` in eq. :eq:eqLGV) (REAL)
 *damping\_therm* : base damping coefficient for thermalization
(:math:`\gamma` in eq. :eq:eqLGV) (REAL)
 *qtb\_frequency\_cutoff* : cutoff frequency for the QTB kernel (REAL)
 *adqtb\_agammas* : (Only for adqtbr and adqtbf) adaptation speed
coefficient for adQTB (:math:`A_\gamma` in eq. :eq:`eqgammadapt`)(REAL)
 *adqtb\_alpha* : (Only for adqtbf) Width of the lorentzian used to
represent the dissipative kernel :math:`\gamma_f(\omega)`
(:math:`\alpha` in eq. :eq:`eqlorentzgenlgv`) (REAL)
 *write\_spectra* : write average random force autocorrelation
function ff, velocity autocorrelation function vv and velocity random
force crosscorrelation function vf spectra (LOGICAL, default=.FALSE.)
 *write\_trajectories* : write Langevin trajectories in x,y,z,px,py,pz
format (LOGICAL, default=.FALSE.)
Remark: all physical quantities are specified in Hartree atomic units.
Output files
============
The Langevin module is plugged to the IVR subroutines and thus can
output the same correlation functions as the classical MC sampling.
Additionally, it can write the Langevin trajectories and spectra
obtained directly from them.
Langevin trajectories

If the parameter *write\_trajectories* of the *langevin* section of the
input file is set to TRUE, Langevin trajectories are saved. Trajectory
files follow the following format:
::
num_of_atoms
At_symbol(1) X Y Z Px Py Pz
At_symbol(2) X Y Z Px Py Pz
.
.
At_symbol(n) X Y Z Px Py Pz
num_of_atoms
At_symbol(1) X Y Z Px Py Pz
At_symbol(2) X Y Z Px Py Pz
.
.
At_symbol(n) X Y Z Px Py Pz
.
.
.
This corresponds to an extended XYZ format with information on momenta.
It is readable by visualization software such as VMD to display the
trajectories.
The module outputs multiple trajectory files depending on the number of
independent trajectories (blocks) and the number of MPI processes. The
naming follows the rules:
 ``xp.traj.xyz`` for 1 block and 1 process
 ``xp_proci.traj.xyz`` for 1 block and multiple processes
 ``xp_proci_blockj.traj.xyz`` for multiple blocks and processes
QTB analysis files

In addition to the trajectories, several files can be edited during the
simulations. They are useful to carefully check the convergence of the
adaptive QTB, notably by calculating :math:`\Delta_{FDT}(\omega)` (eq. :eq:`eqDFDT`).
 ``ff_vv_vf_spectra.out`` spectra of random force and velocity
autocorrelation and random force velocity crosscorrelation functions
Parameters for this potential are specified in an external text file.
The file name is given in the input file using the parameter
*lennard\_jones\_parameters* in section *system*. The parameters to
specify are:
 *epsil* : depth of the potential well :math:`\epsilon` (in Kelvin)
(eq. :eq:`eqLJ_pot`)
 *sigma* : distance for which the potential cancels :math:`\sigma` (in
Å) (eq. :eq:`eqLJ_pot`)
 *r\_cont* : minimum distance for which a confining potential
:math:`r_{cont}` defined in eq. :eq:`eqLJ_cont` is applied (in Å)
The QTB and both adaptive methods were tested on a Ne13 cluster in order
to reproduce results from reference [Man]_.
The LennardJones parameters which have been used are
:math:`epsil=34.9`, :math:`sigma=2.78` and :math:`r\_cont=10.` 5 runs of
8000 steps with 16000 initial time steps are used with all four methods
(Langevin, QTB, adQTBr,adQTBf). Damping term is set to 5.0e5 atomic
units and adaptive coefficients :math:`A_\gamma` and :math:`\alpha` for
adQTBf to 5.0e6 atomic units. Pair correlation function is then
computed from the trajectories output with a Python script
``compute_g2r.py``. Results are shown in figure :numref:`fig_Ne13g2r` and are in
agreement with the ones of ref. [Man]_.
.. _fig_Ne13g2r:
.. figure:: Ne13_g2r.png
Pair correlation function of Ne\ :math:`_{13}` cluster obtained with
Langevin, QTB, adQTBr and adQTBf implemented with Langevin module
in PaPIM. Reference curve calculated with Path Integral Molecular
Dynamics (PIMD)
In this particular case, adaptive QTB leads to significantly better
results than both classical Langevin and QTB when comparing them to the
reference results obtained with PIMD (Path Integral Molecular Dynamics).
Implementation
==============
Langevin module is built with the fewest modifications possible in the
main and previous code of PaPIM. The main program of the sampler is in
the file ``langevin.f90``. It is structured in the same fashion as the
existing samplers (``PIM.f90`` and ``ClassMC.f90``) and only provides
the subroutine *langevin\_sampling* to the main program.
Source files

The Langevin module is divided in multiple files:
 ``langevin.f90``: contains the Langevin sampler and links the main
code with the other files of the module
 ``langevin_integrator.f90``: subroutines to integrate Langevin
equations
 ``langevin_analysis.f90``: spectral analysis tools for Langevin and
(ad)QTB trajectories
 ``qtb_random.f90``: generation of QTB colored noise and adaptation
subroutines for adQTB
Other modifications

Some other routines have been modified during the implementation of
Langevin module.
 ``PaPIM.f90``: main code ; add calls to Langevin module
 ``GlobType.f90``: add declarations for Langevin
 ``ReadFiles.f90``: read input files
Outlook
=======
The next step in the implementation is to add the Wigner Langevin
dynamics to sample the initial conditions for the IVR.
Furthermore, one should note that the current implementation is
functional with the preCP2K version of PaPIM. Thus, some work must be
done to ensure that everything is compatible with the last version of
the code.
Compiling
...
...
@@ 105,6 +490,11 @@ __________
.. [Man] E. Mangaud, S. Huppert, T. Pl'e, P. Depondt, S. Bonella, F. Finocchi, Quantum thermal bath with enforced fluctuationdissipation theorem for reliable simulations of nuclear quantum effects, Journal of Chemical Theory and Computation, Submitted (2018).
.. [Bri] F. Brieuc, Y. Bronstein, H. Dammak, P. Depondt, F. Finocchi, M. Hayoun, Zeropoint energy leakage in quantum thermal bath molecular dynamics simulations, J. Chem. Th. Comput. 12 (2016) 5688–5697.
.. [Hern] J. Hern'andezRojas, F. Calvo, E. G. Noya, Applicability of Quantum Thermal Baths to Complex ManyBody Systems with Various Degrees of Anharmonicity, Journal of Chemical Theory and Computation 11 (2015) 861–870.
.. [Lei] B. Leimkuhler, C. Matthews, Rational Construction of Stochastic Numerical Methods for Molecular Sampling, Applied Mathematics Research eXpress (2012).