The <code>sigma_nu</code> keyword is an '''optional''' keyword in the <code>[[input]]</code> file, unsed in combination with <code>smooth = .true.</code>. It takes 1 float argument, in units of <math>ps^{-1}</math>. It controls the locality of the smoothing applied to the DoS when the LOWESS filter is enabled (<code>smooth = .true.</code>).

=== Summary ===

{| class="wikitable"
|+ Summary for <code>sigma_nu</code> keyword
! style="text-align:left;"| Required/optional
! style="text-align:left;"| Type
! style="text-align:left;"| Accepted values
! style="text-align:left;"| Default
! style="text-align:left;"| See also
|-
| Optional
| Float
| Any positive float
| 1/[[tau]]
| [[Input file]], [[smooth]]
|}

=== Example ===

smooth = .true.
sigma_nu = 2.

File:Smooth 01.png

2021-06-24T07:40:57Z

Miguel Caro:

Smooth

2021-06-24T07:40:46Z

Miguel Caro: Created page with "The <code>smooth</code> keyword is an '''optional''' keyword in the <code>input</code> file. It takes 1 logical argument, <code>.true.</code> or <code>.false.</code>. When..."

The <code>smooth</code> keyword is an '''optional''' keyword in the <code>[[input]]</code> file. It takes 1 logical argument, <code>.true.</code> or <code>.false.</code>. When used, it applies a [https://en.wikipedia.org/wiki/Local_regression LOWESS filter] to smooth out the DoS. This can help stabilize entropy calculations with insufficient sampling due to the noise in the DoS. The locality of the smoothing is optionally tuned with the [[sigma_nu]] keyword. See an example below:

[[File:Smooth_01.png|center]]

=== Summary ===

{| class="wikitable"
|+ Summary for <code>temperature</code> keyword
! style="text-align:left;"| Required/optional
! style="text-align:left;"| Type
! style="text-align:left;"| Accepted values
! style="text-align:left;"| Default
! style="text-align:left;"| See also
|-
| Optional
| Logical
| .false.|.true.
| .false.
| [[Input file]], [[sigma_nu]]
|}

=== Example ===

smooth = .true.

Estimate velocities

2021-06-24T07:28:18Z

Miguel Caro: Created page with "The <code>estimate_velocities</code> keyword is an '''optional''' keyword in the <code>input</code> file. It takes 1 logical argument, <code>.true.</code> or <code>.false...."

The <code>estimate_velocities</code> keyword is an '''optional''' keyword in the <code>[[input]]</code> file. It takes 1 logical argument, <code>.true.</code> or <code>.false.</code>. It is used to enable a quadratic interpolation of the velocities from the positions, useful when the velocities in the trajectory file are not reliable (e.g., positions and velocities are not synchronous). See the discussion in [[Tutorial 1: standard molar entropy of liquid water#Velocity estimation vs. velocities from trajectory|tutorial 1]] for more details.

=== Summary ===

{| class="wikitable"
|+ Summary for <code>temperature</code> keyword
! style="text-align:left;"| Required/optional
! style="text-align:left;"| Type
! style="text-align:left;"| Accepted values
! style="text-align:left;"| Default
! style="text-align:left;"| See also
|-
| Optional
| Logical
| .false.|.true.
| .false.
| [[Input file]]
|}

=== Example ===

estimate_velocities = .false.

File:Tutorial01 04.png

2021-06-24T07:21:58Z

Miguel Caro:

Tutorial 1: standard molar entropy of liquid water

2021-06-24T07:21:43Z

Miguel Caro:

''Note added on 24/06/2021: This tutorial has been updated to accommodate defaults in newer Gromacs versions and includes more information on the issues to keep in mind when using the velocity estimation option. We have also updated the Python code to accommodate Python3 syntax.''

The first tutorial is the simplest calculation you can carry out with DoSPT: estimating the standard molar entropy of liquid water.

=== Generate the trajectory ===

'''Note:''' ''If you have already generated a trajectory file, you can skip this step. Also note that Gromacs (or else) topology files are '''not''' required to run a DoSPT calculation: you only need a trajectory file generated with your MD package of choice (*.xyz and *.gro are currently supported formats) and DoSPT input files.''

DoSPT cannot generate a molecular dynamics trajectory, therefore we need to generate our MD trajectory using some other package. At the moment DoSPT has native support for .xyz and .gro trajectory files. Since [http://www.gromacs.org Gromacs] is an easy-to-use and freely available package, we will be generating our trajectory with Gromacs. We will assume that you are already familiar with MD calculations.

[[File:100_water_molecules.png|thumb|200px|Initial configuration with 100 TIP4P water molecules]]

To avoid wasting time on equilibrating the system, you can download this pre-equilibrated system made up of [http://dospt.org/tutorials/tutorial01/water.gro 100 water molecules]. The model used is TIP4P, which contains one dummy atom per water molecule. Therefore each water molecule is assigned a total of 4 atoms (three real and one dummy). The file name is "water.gro".

Now we generate a trajectory from this initial configuration. We are not going to worry about the quality of the dynamics for this tutorial, since we only need an example trajectory for post-processing with DoSPT. The Gromacs files needed are "topol.top":

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/tip4p.itp"

[ system ]
Water

[ molecules ]
SOL 100

And an "md.mdp" file with some basic options:

integrator = md
;integrator = md-vv ; use this to avoid issues with asynchronous velocities and positions
dt = 0.001
nsteps = 20000
constraints = h-bonds
constraint_algorithm = lincs
lincs_order = 4
nstxout = 4
nstvout = 4
cutoff-scheme = group
ns_type = grid
nstlist = 5
rlist = 0.7
coulombtype = PME
rcoulomb = 0.7
vdw-type = cut-off
rvdw = 0.7
tc-grps = SOL
Tcoupl = nose-hoover
tau_t = 0.2
ref_t = 298.15

Note that we are running dynamics for a total of 20ps and we are saving position and velocity information every 4 steps, which is every 4fs in this case given our 1fs time step. It is also very important to note that if we are using the "md" integrator, which is a "leap-frog" integrator, velocities and positions are offset by half a time step, and cannot be used to accurately estimate angular momenta. We will compare below between "md" and "md-vv", which uses the Velocity Verlet integrator and allows us to directly use velocities from the simulation.

We generate our topology files and run the dynamics:

gmx grompp -f md.mdp -c water.gro -p topol.top -o md.tpr

gmx mdrun -v -deffnm md

After a few seconds the simulation will be finished. We now need to extract the required trajectory information to a file called "traj.gro":

gmx trjconv -f md.trr -s md.tpr -o traj.gro -dt 0.004 -ndec 5 -vel

Note that we are extracting velocity and position information every 0.004ps = 4fs, and we are requesting 5 decimal places for positions (the default 0.001 nm precision of Gromacs is not enough for the 2PT analysis). If we ran with the "md" integrator, we should not use the velocities because of the mentioned issue with leap-frog, so we could as well have left the "-vel" option out.

The generated trajectory file is quite large, 132MB in this case. You'll have to take memory limitations into account when extracting trajectory info.

=== Prepare the DoSPT input files ===

Now that we have generated a trajectory DoSPT comes into play. We need 5 files to run this simulation, which should all be located within the same directory:

* traj.gro (trajectory file from the MD simulation)
* input (all the DoSPT input options)
* groups (defines the groups of atoms)
* supergroups (all groups of the same type)
* masses (defines atomic masses)

The following input file will be enough for our purposes:

# Number of points in the trajectory (5001)
points = 5001
# Trajectory period in ps (20)
tau = 20.
# Size of the box in nm
cell = 1.472 1.472 1.472
# Temperature in K
temperature = 298.15
# Trajectory info
format = gro
# Estimate velocities from positions
estimate_velocities = .true.

Note that we had to enable velocity estimation because of the problem with leap-frog integration. The cell size can be taken from the last line in the "water.gro" file. The units in this input file are strictly given in the units specified in the comments. Comments always start with a "#" symbol.

The "masses" file assigns a mass in a.m.u. to the atomic labels. For TIP4P, atomic labels are "OW", "HW1", "HW2" and "MW", for oxygen, the two hydrogens, and the dummy atom, respectively. Our "masses" file looks like this:

OW 16.00
HW1 1.01
HW2 1.01
MW 0.00

We have assigned zero mass to the dummy atom. You will need to adapt the mass definitions depending on the specific topology of your molecule when dummy atoms are included.

Now we need to generate the "groups" file. The mandatory structure is as follows:

natoms ngroups
n1 s1
n1_1 n1_2 n1_3 ...
n2 s2
n2_1 n2_2 n2_3 ...
...

The first line contains the total number of atoms "natoms" in the system followed by the number of groups "ngroups". This is then followed by two lines for each group, one line with the number of atoms in the group ("n1", "n2", etc.) and the symmetry number of the group ("s1", "s2", etc.) and another line with the labels of all the atoms belonging to that group. In our present example, TIP4P water molecules have 4 atoms, and their symmetry number is 2; therefore our "groups" file will look like this:

400 100
4 2
1 2 3 4
4 2
5 6 7 8
...

The easiest way to generate this file is with a Python script "setup_groups.py", e.g.:

print(400, 100)

for i in range(0,100):
print(4,2)
print(4*i+1, 4*i+2, 4*i+3, 4*i+4)

Then run:

python3 setup_groups.py > groups

Finally, we need a "supergroups" file. Since this is a monocomponent system, all the groups, from 1 to 100, belong to the same supergroup. The "supergroups" file can be constructed like this:

1-100

=== Run DoSPT and analyze the results ===

At this point we have everything we need to run a DoSPT calculation. Simply type "DoSPT" and look at the program's progress being printed on the screen. At the end of the execution several files will have been created. The entropy information is contained in the "entropy" file. Type "tail -2 entropy" to get this info printed on the screen:

# Supergroup; Entropies in eV/K [trans, rot, vib]; Total entropy in eV/K; Total entropy in J/K
1 5.46004868E-02 1.05552205E-02 1.56976588E-06 6.51572771E-02 6.28682368E+03

The last number tells you that the total entropy for this system is 6286.2 J/K, if this system constituted one molecule of the substance under study. Since we have 100 molecules, this means that the calculated molar entropy is 62.862 J/K/mol. Compare this to the experimental value of 69.9 J/K/mol. Note that the discrepancy with experiment is mostly due to the water model, force field used and simulation box size (finite-size effects).

The density of states for supergroups is written to file "dos_sg". You can visualize it using for instance "gnuplot":

set xlabel "Frequency (ps^{-1})"
set ylabel "DoS (ps)"
plot "dos_sg" u 1:2 t "Trans" w l, "" u 1:3 t "Rot" w l, "" u 1:4 t "Vib" w l

[[File:tutorial01_01.png|center]]

As expected there are no vibrational contributions to the DoS since TIP4P is a rigid model. Since there are no high frequency vibrational contributions, we can plot a close-up:

[[File:tutorial01_02.png|center]]

Since the integral of the DoS should equal the total number of degrees of freedom (DoF), we can use the "smooth cumulative" option of gnuplot to get a quick estimate of this number (the distance between data points is 0.05 ps^-1):

set key bottom
set ylabel "Integrated no. of DoF"
plot "dos_sg" u 1:($2*0.05) smooth cumulative t "Trans DoF" w l, "" u 1:($3*0.05) smooth cumulative t "Rot DoF" w l

[[File:tutorial01_03.png|center]]

It's easy to see that the computed number of DoFs is less than expected: 300 translational and 300 vibrational - adjusted to the specifics of your MD code, which might be subtracting the linear and/or angular velocities of the center of mass of your system (i.e., minus 3 degrees of freedom for each). The problem of underestimation of the number of degrees of freedom is associated to equilibration, total trajectory time and system size, and severely affects the rate of convergence of calculated entropy results. Check out the [[DoS renormalization]] article to learn how DoSPT can correct the leading error and significantly speed up convergence of results with respect to the number of sampled configurations.

=== Velocity estimation vs. velocities from trajectory ===

As mentioned above, with the leap-frog integrator the velocities from Gromacs are useless, because they are offset by half a time step from the positions. This means that the calculated angular momenta would be wrong. '''DoSPT''' allows you to use the [[estimate_velocities]] keyword to discard the velocities given in the trajectory file and use a quadratic interpolation to estimate them from the positions. The problem with this is that you need very fine sampling to get accurate velocities with this method, leading to large trajectory files. Instead, the preferred method is to use a Velocity Verlet integrator ("md-vv" in Gromacs) to get synchronous velocities and positions out of the trajectory file. This is more accurate and requires less fine sampling. Below are examples with ("md") and without ("md-vv") velocity estimation for different sampling frequencies:

[[File:tutorial01_04.png|center]]

The issue is clear. Particularly, high-frequency modes are suppressed by the velocity smoothing, which should be used with caution.

=== Homework ===

Repeat this tutorial with other water models to get an idea of the impact of the force field on the results. Also, repeat the calculation for different starting configurations/number of molecules (with the same force field) to get an idea of the statistical variation of entropy results and the impact of system size.

File:Tutorial01 03.png

2021-06-24T06:59:19Z

Miguel Caro: Miguel Caro uploaded a new version of File:Tutorial01 03.png

File:Tutorial01 02.png

2021-06-24T06:58:34Z

Miguel Caro: Miguel Caro uploaded a new version of File:Tutorial01 02.png

2020-04-20T16:14:22Z

Miguel Caro:

This is a list of publications that highlight the capabilities and/or make use of the DoSPT code. Papers by DoSPT authors are in '''bold'''.

== 2020 ==

S. Mukherjee, B. Bagchi; 
''Entropic Origin of the Attenuated Width of the Ice–Water Interface''; 
[https://pubs.acs.org/doi/abs/10.1021/acs.jpcc.0c02030 J. Phys. Chem. C 124, 7334 (2020)].

== 2018 ==

S.R. Varanasi, Y. Subramanian, S.K. Bhatia; 
''High Interfacial Barriers at Narrow Carbon Nanotube--Water Interfaces''; 
[https://pubs.acs.org/doi/abs/10.1021/acs.langmuir.8b00616 Langmuir 34, 8099 (2018)].

== 2017 ==

'''Miguel A. Caro, Olga Lopez-Acevedo, and Tomi Laurila'''; 
''Redox Potentials from Ab Initio Molecular Dynamics and Explicit Entropy Calculations: Application to Transition Metals in Aqueous Solution''; 
[http://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00314 J. Chem. Theory Comput. 13, 3432 (2017)].

== 2016 ==

'''Miguel A. Caro, Tomi Laurila, and Olga Lopez-Acevedo'''; 
''Accurate schemes for calculation of thermodynamic properties of liquid mixtures from molecular dynamics simulations''; 
[http://aip.scitation.org/doi/pdf/10.1063/1.4973001 J. Chem. Phys. 145, 244504 (2016)].

2019-05-29T13:11:47Z

Miguel Caro:

A <code>traj.xyz</code> [[trajectory file]] should adhere to the following convention:

Number_of_atoms
Name_or_comment
Atom_name_1 posx posy posz (velx vely velz)
Atom_name_2 posx posy posz (velx vely velz)
...
Atom_name_nat posx posy posz (velx vely velz)

Where the units for positions are Angstrom and the units for velocities (which are optional if you use the <code>[[estimate_velocities]] = .true.</code> option) are Angstrom/fs.

For a 100-molecule TIP4P water system (including dummy atoms) this looks something like this (without velocities):

400
Water
OW 7.50 8.43 5.11
HW1 6.76 8.30 4.53
HW2 7.96 9.19 4.76
MW 7.47 8.51 4.99
OW 5.95 4.22 0.07
HW1 6.69 4.82 0.12
HW2 6.18 3.63 -0.65
MW 6.08 4.22 -0.02
.
.
.
OW 2.70 8.23 4.53
HW1 2.47 7.54 3.91
HW2 1.92 8.35 5.07
MW 2.57 8.16 4.52