The Basics of Molecular Simulations: Part-2
In Part-1, I talked about the fundamental idea behind molecular simulations and its applications in the field of material science and drug discovery in brief. I explained how ensembles are defined to mimic the experimental conditions to conduct molecular simulations. Furthermore, I described Grand Canonical Monte Carlo (GCMC) and Molecular Dynamics (MD) simulation techniques with used case scenarios. For those, who didn’t check the Part-1, please read it before going further.
In Part-2, I will talk about the system definition and input forcefield parameters needed to carry out these simulations. I have used molecular simulations for materials, so I will talk about it from that perspective. I have very limited knowledge about biological systems but still, I will add a note on whatever little I know of.
System Definition
Crystallographic Information File (CIF)
To initiate molecular simulations, we need to provide the crystal structure of the molecular system as an input. This information is stored in a text file called CIF that contains information about the crystal system, cell lengths, cell angles, cell volume and coordinates of the atoms constituting the system, etc. These details are obtained from the experiments. Here is the link for some sample files.
The Cambridge Structural Database (CSD) is a database that has over one million structures.
Note: For biological systems, all the known protein structures are stored in the Protein Data Bank (PDB). When we deal with proteins, we need PDB files instead of CIF files.
Force Field Parameters
As I explained in Part-1, potential energy is calculated between the atoms to run molecular simulations and to calculate it, certain parameters are required which are called forcefield parameters.
The interactions between atoms can be classified into two categories, bonded and non-bonded. Van der Waals and electrostatic interactions are non-bonded. Bonded interactions (bond, angle bending, and dihedrals/torsional)are the one which involves movements of atoms around bonds.
Non-bonded parameters
1. Atomic partial charges (Coulombic or electrostatic interactions):
To calculate electrostatic interactions, charges need to be assigned to the atoms. The most adopted approaches to calculate charges are based on density functional theory (DFT) calculations such as Mulliken population analysis, density derived electrostatic and chemical (DDEC) analysis and electrostatic potential analysis. They are highly accurate and reliable but they require a lot of computational time. These methods are good to use when we are working with a smaller system but pose a significant problem when we are dealing with a large system or many systems all at a time.
E.g., my research was to screen and assess the adsorption properties of many Metal-organic frameworks (MOFs). So, I used another method called charge equilibration (Qeq). It is not very accurate but requires comparatively very low computational power and time. So, there is a tradeoff involved between accuracy and time. My advice is, if the system is not very sensitive to atomic charges or you need to get faster results, then the Qeq method is a better choice. More information about it can be found in the paper below.
Note: For biological systems, atomic charges for all the protein structures are already known so no need to worry about deriving charges for proteins.
2. Van der Waals parameters:
Van der Waals interactions are also termed as Lennard-jones (LJ) potential. ε and σ are LJ parameters. Generic or classical force fields such as the Universal force field (UFF) and Dreiding are the most extensively utilized force fields for molecular simulations. L-J parameters reported in these forcefields have been validated rigorously. All the elements have standard parameters that can be used for any system consisting of these elements.
Note: CHARMM forcefield is widely used for biological systems.
Bonded parameters
As I have already described earlier that every system is composed of atoms connected. These atoms are linked to each other like a mass-spring system (Mass represents the atoms and the spring represents the bonds between them). Bonded parameters represent the spring and equilibrium constants. Bonded interactions consist of the bond, angle bending and torsional potentials.
Bonded parameters are also obtained from generic force fields, such as CVFF, Universal force field (UFF), or Dreiding force field but not all the parameters are available. So the unavailable parameters are derived or tuned using an empirical fitting approach or quantum mechanical (QM) calculations or experimental data. As it is obvious that bonded parameters depend on the molecular geometry and will be different for each geometry, unlike non-bonded parameters which are defined for individual atoms, independent of the configuration.
Recently, a new forcefield known as Reactive force field (ReaxFF) is becoming popular because it can be applied for a variety of systems having nonidentical molecular geometry.
Summary
References:
In the next part, I will cover the software and other tools which I used to perform molecular simulations.
Thank you for reading. I think, as a member of the science community it is our responsibility to pass on the information which we acquire, to the fellow researchers. This would save a lot of time for them which would be very valuable for their research. If you want to be a part of this initiative, please let me know. Connect with me here. If you liked this article, you may also like the article below.
