Molecular Mechanics and Molecular Dynamics

Notes are in inverse chronological order: newest first

# Fit MM3 dihedral potentials

In the previous note, I wrote how derive or improve force field potentials. I used the test case of improving the force field description for the interring torsion in bithiophene molecules. The step 4 of that procedure was the fitting of the energy different between QC and FF (without the torsional potentials to improve) PES. Those data must be fitted with the function used in the force field to describe that potential. In case of MM3 force field (as implemented in TINKER), this function is (in gnuplot format):

f(x)=a+0.5*(b*(1+cos(x-e))+c*(1+cos(2*x-f))+d*(1+cos(3*x-g)))

The fit is done via the shift parameter a (not used in the force field), the force field parameters b, c and d and the cosine phases e, f and g. Since the phase is constrained to be between zero and 360°, the general formula to get the proper value for the phase $\phi$ is:

$\phi = (\frac{abs(\phi)}{360}-int(\frac{abs(\phi)}{360}))*360$

Potentials and phases can be used in the new dihedral potential.

# Derive force field parameters from Quantum Chemical calculations: The Bithiophene test case

The results and the accuracy of MM and MD calculations strongly depend on how good the force field is to describe the properties of the molecular systems under investigation. The choice of a proper force field is, therefore, crucial. Since force field are parameterized on a large number of a certain kind of molecules (e.g., proteins, biological system, water, organic molecules, etc.) to be as general as possible, some parameters can be not so accurate to reproduce a particular property. Or, simply, can happen some parameters are missing. In these cases, quantum chemical calculations can be used to derive missing force field parameters or improve them.

In this scrapnote I will show how a way to derive (or improve) the torsional parameters for the rotation of the two thiophene rings around the central bond in the bithiophene molecules. Figure 1 shows a bithiophene molecule (2T) in its syn and anti conformations.

Fig. 1. Syn (or s-gauche-cis) and anti (or s-gauche-trans) conformations of bithiophene molecule. In the syn conformation, a numerical index has been assigned to the atoms involved in the interring torsion.

Referring to the syn conformation in figure 1, the torsion we want to improve is that described in the forcefield as the sum of the following four dihedral potentials: 1-2-5-4, 1-2-5-6, 3-2-5-4 and 3-2-5-6. The dihedrals 1-2-5-6 and 3-2-5-4 are described by the same set of parameters, since they involve the same kind of atoms.

Experimentally, it is possible to measure the torsion for the syn and anti conformations of 2T in gas phase as 37º ± 5º and 148º ± 3º, respectively. Also the relative populations of the two minima is obtained experimentally (see Figure 2-Top): since the population of the anti conformation is the biggest, this conformation is the absolute minima. The comparison between different QC potential energy surfaces (PESs) for the torsion calculated at B3LYP and MP2 theory levels with different basis sets is shown in Figure 2-Bottom.

Fig. 2. Experimental positions and populations of the two minima for the interring torsional potential of 2T in vacuo (top) and systematic comparison of PESs calculated using B3LYP and MP2 methods with different basis sets (bottom).

The position and the relative energy of the minima, are not the only important properties of the PES: the location of the transition state between the minima and the high of the energy barrier are also crucial aspect to consider to properly describe this torsional potential. As shown in Figure 3, many force fields lack in the description of the syn local minima and, mostly, in the high of the rotational barrier for the transition state at 90º.

Fig. 3. Comparison between the B3LYP 6-31++g** PES and different force fields: mm2 (in different variants), cvff_aug and burchart-dreiding. None of them is able to well reproduce the key properties of the interring torsional potential in 2T molecules in vacuo.

To improve the torsion using the QC PES is quite straightforward and the simple scheme in Figure 4 can be used to illustrate the procedure:

Fig. 4. Schematic procedure to derive or improve force field parameters from quantum chemical calculations. The number of the step indicated in orange.

The first step consists in calculate the PES with QC calculations using a good method and basis set. A systematic comparison can be done as shown in Figure 2 and Figure 4 (step 1). Once the target QC PES is calculated, is time to move to the step 2, which consists in obtain the PES for the torsion at force field level. For this step, we need to suppress all the potential terms which involve the torsion we want to describe. In this case, the dihedral potentials for the dihedral angles 1-2-5-4, 1-2-5-6 and 3-2-5-4 and 3-2-5-6 are set to zero. At this point, the difference between the two PESs (step 3) can be considered as the new torsional potential to use in the force field. Indeed, in this energy difference all the differences from the QC and FF calculations are included. The torsional potential obtained in step 3, need to be fitted with the functional form used in the force field to describe dihedrals1 (VFF(x) in step 4). The fitting will provide the numerical values for the parameters to put in the force field2. The last two steps consist in redo a torsional scan with the complete force field and compare the result with the QC one3: Figure 4 shown now a complete agreement between the QC and FF PES for the interring torsion of 2T in vacuo.

NOTE: the images and the plots have illustrative purpose only: you should not use them for your research.

# Thermostats and Barostats: how to use them?

This is my personal scheme about “when” and “how” apply thermostats and barostats in molecular dynamics:

Algorithms:

• Berendsen: good for equilibration, less for NVT simulations, since does not sample properly the NVT ensemble.
• Nose'-Hoover ; Andersen: sample correctly the NVT ensemble. To use in simulations to collect statistic.
• Nose'-Hoover + Parrinello-Raman: for NPT simulations

Coupling constants:

• equilibration: tau_t <= 0.01ps ; tau_p=0.5-1.0ps
• md runs: tau_t=0.1-0.5ps; tau_p=5-10ps