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.

protein.inp

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:

aimd.inp

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:

protein_qmmm.inp

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.