Tip

All input and output files can be downloaded `here`

.

# Molecular Dynamics: Standard MD

Molecular dynamics (MD) is an important way of exploring dynamic pictures of chemical and biological phenomena. There are many MD methods in Qbics:

Standard MD: Calculating trajectories by integrating Newtonian equations of motion using forces from QM, MM, or QM/MM potentials;

Enhanced-sampling MD: Calculating trajectories by integrating Newtonian equations of motion using forces from QM, MM, or QM/MM potentials, and a

**biased potential**for realizing “random walking” on assigned collective variables (CVs). This is important for free energy calculations.Free-energy perturbations (FEPs). This is a special way to calculate free energy changes for two different configurations (there may be different numbers of atoms for the two configurations).

In this tutorial, we only describe standard MDs.

## The First Classic MD Example

We will use the protein system in Molecular Mechanics Energy to do a MD. Since the energy and force are calculated using CHARMM potential, this is a classic NVT MD.

```
charmm
parameters toppar_water_ions.str par_all36m_prot.prm
topology step2_solvator.psf
scaling14 1.0
rcutoff 12.0
rswitch 10.0
use_PBC
Lbox 64 64 64
electrostatic pme
PMEk 64 64 64
end
md
type NVT
dt 0.001
temp 300.
temp_nhc_length 5
temp_nhc_tau 0.05
num_steps 50000
output_freq 1000
end
mol
protein-opt.xyz
nd
ask
md charmm
nd
```

MD is controlled by `md ... end`

block. The meaning of options are shown below:

`type`

The type of MD. In this section, you can use`NVE`

(constant number of particles, volume, and energy, microcanonical ensemble) or`NVT`

(constant number of particles, volume, and temperature, canonical ensemble).`dt`

The time-step for MD in ps. For systems containing hydrogen atoms, it should be smaller than 0.001 ps (1 fs). If only heavy atoms exist, one may set a larger one like 0.002 ps (2 fs).`temp`

The temperature in Kelvin. For NVE simulations, this is the initial temperature. For NVT simulations, this is the target temperature that the system is controlled at. In this program, temperature is controlled using Nose-Hoover chain method.`temp_nhc_length`

The chain length in Nose-Hoover chain method. It should be an integer greater than 1.`temp_nhc_tau`

The time constant in Nose-Hoover chain method in ps. The larger this value is, the faster Nose-Hoover chain method will “pull back” the system to the target temperature.`num_steps`

The number of steps that the MD will be performed.`output_freq`

The frequency for MD result output.

Attention

In `protein.inp`

, we do not use the original coordinate file `step2_solvator.pdb`

since this structure is of high energy and MD will definitely crash. So, we use the optimized structure obtained in Geometry Optimization, i.e. `protein-opt.xyz`

. In classic MD, **one must optimize the structure before doing MD.**

Now we can run the MD simulation:

```
$ qbics protein.inp -n 4 > protein.out
```

After simulation, you can find 5 files:

`protein-md-traj.xyz`

The atomic coordinates of each output frame.`protein-md-vel.xyz`

The atomic velocities of each output frame.`protein-md-grad.xyz`

The atomic gradients of each output frame.`protein-md-log.txt`

Some physical quantities (like energy components, temperature) of each output frame.`protein-md-cv.txt`

It is used for free energy calculations.

## An Ab Initio MD Example

You can do an ab initio MD (AIMD) using the input in Density Functional Theory Calculations. For example, the following is an NVT-AIMD at B3LYP-D3BJ/def2-SVP level of theory:

```
basis
def2-svp
end
scf
charge 0
spin2p1 1
type r
end
grimmedisp
type bj
end
md
type NVT
dt 0.0005
temp 300.
temp_nhc_length 5
temp_nhc_tau 0.05
num_steps 100
output_freq 1
end
mol
C -1.72885348 -0.27751594 -0.54883607
C -0.54335149 -0.47251267 0.35111257
C -1.54903120 0.59933911 0.65616163
H -1.57461433 0.15313436 -1.55374643
H -2.53379418 -1.03326873 -0.54728357
H -0.50081907 -1.36760843 0.99612354
H -2.22562112 0.47045579 1.51913545
H -1.26635100 1.65697053 0.51257787
C 0.85387745 -0.06614383 -0.15308545
O 1.09143365 0.21649289 -1.34006650
N 1.82735682 -0.03716326 0.80612988
H 2.75631019 0.21299272 0.57405353
H 1.62607490 -0.26347660 1.74464638
end
task
md b3lyp
end
```

For AIMD, the time step should be **smaller than 1 fs**, in this case, `dt 0.0005`

which is 0.5 fs.

## An QM/MM MD Example

For both examples in Quantum Mechanics/Molecular Mechanics (QM/MM) Energy, just change `energy`

to `md`

and add `md ... end`

block to configure MD options. For example, for ubiquitin, the following input does an NVT-MD at PHO-B3LYP-D3BJ/def2-SVP/CHARMM level of theory:

```
charmm
parameters par_all36m_prot.prm
topology step1_pdbreader.psf
scaling14 1.0
rcutoff 12.0
rswitch 13.0
end
basis
def2-svp
end
scf
charge 0
spin2p1 1
type R
end
grimmedisp
type bj
end
qmmm
qm_region 779 792-815
end
md
type NVT
dt 0.0005
temp 10.
temp_nhc_length 5
temp_nhc_tau 0.05
num_steps 20
output_freq 1
end
mol
protein-qmmm-opt.xyz
end
task
md b3lyp/charmm
end
```

Note that, as done for the protein system in classic MD, we have used the optimized structure `protein-qmmm-opt.xyz`

as the initial one for MD.