.. tip::
    All input and output files can be downloaded :download:`here <_static/md/md.zip>`.

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 :doc:`mmenergy` to do a MD. Since the energy and force are calculated using CHARMM potential, this is a classic NVT MD.

.. code-block:: 
    :caption: 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 :doc:`geoopt`, i.e. ``protein-opt.xyz``. In classic MD, **one must optimize the structure before doing MD.**

Now we can run the MD simulation:

.. code-block:: bash

    $ 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 :doc:`qmenergy`. For example, the following is an NVT-AIMD at B3LYP-D3BJ/def2-SVP level of theory:

.. code-block:: 
    :caption: 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 :doc:`qmmmenergy`, 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:

.. code-block:: 
    :caption: 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.