4. Running molecular dynamics

This functionality of the cp2kmdpy package can be used to run molecular dynamics simulations.

class cp2kmdpy.md.MD(molecules=None, functional=None, box=None, cutoff=None, scf_tolerance=None, basis_set=[None], basis_set_filename=None, potential_filename=None, fixed_list=None, periodicity=None, simulation_time=None, time_step=None, ensemble=None, project_name=None, temperature=None, pressure=None, n_molecules=None, thermostat=None, traj_type=None, traj_freq=None, seed=None, input_filename=None, output_filename=None, initial_coordinate_filename=None, use_atom_name_as_symbol=True)[source]

Base class for running molecule optimization.

Parameters:
  • molecules – Molecule(s) in the simulation box
  • functional (string) – DFT XC functional to be used
  • box (mBuild box) – Box in which the molecules should be placed
  • cutoff (float, optional) – Plane wave cutoff (Ry) for DFT calculation
  • scf_tolerance (float, optional) – Tolerance for each SCF cycle
  • basis_set (dictionary with key being the atomic symbol and value being the basis set, both strings) – Basis set for each atomic kind
  • basis_set_filename (string, optional, defaults to BASIS_MOLOPT) – Filename for the basis set
  • potential_filename (string, optional, defaults to GTH_POTENTIALS) – Filename for the pseudopotential to be used
  • fixed_list (string, optional, for example: if atom # 1 to 10 shall be fixed => fixed_list='1..10') – list of elements to be frozen
  • periodicity – Periodicity of the box
  • simulation_time (unyt quantity, optional) – Total time of the simulation
  • time_step (unyt qunatity, example: time_step=0.1*u.fs # when time step is 0.1 femtosecond) – Time step for integrating the trajectory
  • ensemble (string, optional) – MD simulation ensemble
  • project_name (string, optional) – Name of the project
  • temperature (unyt quantity) – Simulation temperature
  • pressure (unyt quantity) – Simulation pressure
  • n_molecules (list of number of molecules of each kind) – Number of molecules in the box of each kind
  • thermostat (string, optional) – type of thermostat
  • traj_type (string, optional) – Output trajectory format
  • traj_freq (integer, optional) – Trjectory snapshot capture frequency
  • seed (integer, optional) – Random number seed for MD simulation
  • input_filename (string, optional) – Name that should be given to the input file
  • output_filename (string, optional) – Name that should be given to the output file
  • initial_coordinate_filename (string, optional) – Structure filename if user is providing their own initial structure for md
  • use_atom_name_as_symbol (boolean, optional) – If you want the atom name to be same as its symbol. Useful if you have want to use different basis sets for a single atom type in different environment. If you are setting it to false, make sure that the particles in your mbuild molecule have particle.element and particle.name attributes set (particle.element will be the symbol and particle.name will be the name used).

Instantiating MD first requires the specification of the molecule(s), functional, box , the basis set, and the number of molecules of each type. It can be done as follows:

import cp2kmdpy
md_instance = cp2kmdpy.md.MD(
    molecules=[molecule],
    functional=functional,
    box=box,
    basis_set=basis_set,
    n_molecules=n_molecules
)

4.1. MD class

4.1.1. molecules

molecules is the list of molecules present in the simulation box. Each element of the list must belong to the mbuild.Compound class.

For example, for a mixture of ethane and propane, molecules can be defined as follows:

import mbuild

ethane = mbuild.load("CC", smiles=True)
propane = mbuild.load("CCC", smiles=True)
molecules= [ethane, propane]

4.1.2. functional

functional is the name of the DFT XC functional to be used for the MD simulation. functional should be a string.

functional = 'PBE'

4.1.3. box

box is the simulation box which will be populated with the molecules. box should be an instance of the mbuild.Box class.

import mbuild
box = mbuild.Box(lengths=[1.0, 1.0, 1.0], angles=[90., 90., 90.])

4.1.4. cutoff

While using the GPW method in CP2K, the plane wave cutoff needs to be defined. cutoff is the plane wave cutoff in Ry for optimization simulation. It should be a float.

cutoff=600

4.1.5. scf_tolerance

scf_tolerance controls the tolerance of each SCF cycle in the optimization. Should be a float.

scf_tolerance=1e-7

4.1.6. basis_set

With basis_set, the basis set to be used for each element present in the simulation can be defined. It is a dictionary with keys being the element symbol and values being the basis set for that particular element. For example, if the molecule contains carbon and oxygen atoms and DZVP-MOLOPT-SR-GTH basis set is desired, then the basis_set can be defined as:

basis_set={'C':'DZVP-MOLOPT-SR-GTH','O':'DZVP-MOLOPT-GTH'}

4.1.7. basis_set_filename

The name of the file to look into for the basis set specified. Should be a string. Default value is BASIS_MOLOPT.

4.1.8. potential_filename

The name of the file to look into for pseudopotential information. Should be a string. Default value is GTH_POTENTIALS.

4.1.9. fixed_list

During molecular dynamics simulations, sometimes constraining the movement of some atoms is desired. This can be achieved by specifying fixed_atoms. For example, if atom number 1 to 80 needs to be fixed, it can be specified as follows

fixed_list='1..80'

4.1.10. periodicity

This attribute controls the periodicity of the box. Default value is 'XYZ' Consider a system which is periodic in only X and Y direction. It can be specified as follows:

periodicity='XY'

4.1.11. simulation_time

Time for which simulation should be run. It should be a unyt_quantity. For example if the simulation should be run for 1 ps , simulation_time can be set as:

import unyt as u

simulation_time=1*u.ps

The default value of simulation_time is 10 ps.

4.1.12. time_step

Time step for the molecula dynamics simulation. Should be a unyt_quantity.

import unyt as u

time_step=1*u.fs

4.1.13. ensemble

Ensemble for the molecular dynamics simulation. Must be a string of capital letters. Deafult value is 'NVT'

4.1.14. project_name

Name of the project. Should be string. Default value is 'sample_project'

4.1.15. temperature

Temperature of the ensemble, if required. Must be specified as a unyt_quantity.

import unyt as u

temperature=298*u.K

4.1.16. pressure

Pressure of the ensemble, if required. Must be specified as a unyt_quantity.

import unyt as u

pressure=1*u.bar

4.1.17. n_molecules

List containing the number of molecules of each kind. Each elemnet of the list must be a positive integer.

4.1.18. thermostat

Type of thermostat to be used. Must be a string. Default value is 'NOSE'.

4.1.19. traj_type

Output trajectory format. Must be a string. Default is 'PDB'.

4.1.20. traj_freq

Frequency of storing simulation trajectory. Must be positive integer.

4.1.21. seed

Seed for randomly initializing velocities. Must be a positive integer.

4.1.22. input_filename

Name you want to give to the input file that is generated using cp2kmdpy. Default name will be depend on the project name.

4.1.23. output_filename

Name you want to give to the output file. Default name will be depend on the project name.

4.1.24. initial_coordinate_filename

In case structure generation from mBuild is not wanted/required, initial coordinates can be specified through this attribute.

4.1.25. use_atom_name_as_symbol

If you want to make the atom name different from its symbol, set this flag to False. In that case, you will have to set the particle.name and particle.element attributes of the particles in the mbuild molecule.

4.2. Running molecular dynamics

  1. Copy the setter.py file from cp2kmdpy/setter.py and change any settings if needed
  2. Instantiate a cp2kmdpy.md.MD class and set all the attributes
  3. Run the md_initialization method on that instance
  4. Generate input files using setter.md_files function
  5. Run using the run_md method from cp2kmdpy.runners

4.3. Examples

4.3.1. Dinitrogen NpT MD

# ### Loading modules
import numpy as np
import unyt as u
import mbuild as mb
from cp2kmdpy.molecule_optimization import Molecule_optimization # for single molecule optimization
from cp2kmdpy.md import MD # for running MD
from cp2kmdpy import runners
import setter


# Defining the molecule we want to simulate


class N2(mb.Compound):
    def __init__(self):
        super(N2, self).__init__()

        nitrogen1= mb.Particle(pos=[0.0, 0.0, 0.0], name='N')
        nitrogen2= mb.Particle(pos=[0.2, 0.0, 0.0], name='N')
        self.add([nitrogen2,nitrogen1])
        self.add_bond((nitrogen2,nitrogen1))

#creating an instance of the MD class, also defining the parametrs of our md simulation
molecule_list=[N2()]
box=mb.box.Box(lengths=[1,1,1])
q=MD(molecules=molecule_list,box=box,cutoff=200,functional='PBE',basis_set={'N':'DZVP-MOLOPT-GTH'},periodicity='XYZ',n_molecules=[2],traj_type='PDB',seed=1,project_name='N2_npt')


#Setting temperature and ensemble, could have also been set while creating an instance of the MD class

q.temperature=273.15*u.K
q.ensemble='NPT_I'
q.simulation_time=10*u.fs
q.pressure=1*u.bar

#Initializing q

q.md_initialization()


#generating input files
setter.md_files(q)


#running md
runners.run_md(q.input_filename,q.output_filename )