I want to generate the dihedral constants for a force field.
Doing a potential energy scan of by rotating around the 2-3 bond using an QM software we can get the variation of the dihedral energy as a function of the dihedral angle, something like what is shown in the figure. This profile can be fit using a suitable function, for example an expansion in cosines.
But when we have to implement this force field in a MD software like LAMMPS (which is what i intend to use), it requires that the dihedral constants be specified for all possible bonds associated with a dihedral. For example in the cartoon, the dihedral constants need to be specified for 1-2-3-4 and also for 1-2-3-5. My question is how do I split the constants obtained from fitting the ab-inito dihedral scan.
Though I have come across a few, I would like to be pointed to some exhaustive review papers on force field parametrization.
dear Salman,
The main issues are three:
- the QM scan should be relaxed, i.e. optimized at each angle, and so it should be the MM scan. You want to stay close to the minimum of the potential energy surface at each angle, and measure there your potential.
- the total MM torsional potential is not given only by the explicit dihedral parameters, but there's also an intramolecular mostly non-bonded contribution given by LJ and electrostatic interactions (but also bond stretching, ...), let's call it U_NB(phi). To measure it, you perform a MM scan with the explicit dihedral parameters all set to zero.
- as you point out, the QM scan gives the sum of the terms present in the force field for a given torsion. For instance the QM scan about the 2-3 bond will give you the potential of 1234+1235, accounting also that they could be phase shifted of an angle gamma. So at the end, what you want to get is
U_QM (phi) = U_1234 (phi) +U_1234 (phi+gamma)+U_NB(phi)
There are no unique solutions in the most general case; one can be of setting the two potentials to be equal, and then modify only one of the two to eventually match U_QM.
You might have a look at the attached tutorial - is it focus on NAMD but most of the concepts should apply for any MD code.
Regards
L
Salman,
Your calculated dihedral potential from QM should include both potentials of 1-2-3-4 and 1-2-3-5. So, the easiest way is to consider all the estimated potential for dihedral of 1-2-3-4 and set the K value equal to zero for dihedral 1-2-3-5. This would gave you the right relative potentials because for your specific molecules 1-2-3-4 an 1-2-3-5 are always associated.
Thanks
For example if i have a molecule like in this figure. There are 4 dihedrals in this
1-2-3-4
1-2-3-5
3-5-6-7
3-5-6-8
Given that I have all other parameters for the force field ( non-bonded terms, bond coefficients and angle coefficients), how do I get constants for all dihedrals?
Do I do a QM scan along a particular dihedral letting the rest of the molecule relax, then do a similar scan using my incomplete force field and fit the difference?
Then repeat the same procedure for the rest of the dihedrals, but including the already fit dihedral coefficients?
dear Salman,
The main issues are three:
- the QM scan should be relaxed, i.e. optimized at each angle, and so it should be the MM scan. You want to stay close to the minimum of the potential energy surface at each angle, and measure there your potential.
- the total MM torsional potential is not given only by the explicit dihedral parameters, but there's also an intramolecular mostly non-bonded contribution given by LJ and electrostatic interactions (but also bond stretching, ...), let's call it U_NB(phi). To measure it, you perform a MM scan with the explicit dihedral parameters all set to zero.
- as you point out, the QM scan gives the sum of the terms present in the force field for a given torsion. For instance the QM scan about the 2-3 bond will give you the potential of 1234+1235, accounting also that they could be phase shifted of an angle gamma. So at the end, what you want to get is
U_QM (phi) = U_1234 (phi) +U_1234 (phi+gamma)+U_NB(phi)
There are no unique solutions in the most general case; one can be of setting the two potentials to be equal, and then modify only one of the two to eventually match U_QM.
You might have a look at the attached tutorial - is it focus on NAMD but most of the concepts should apply for any MD code.
Regards
L
Salman, your conclusion was what I meant. The problem with the method I suggested is that is only limited to the molecule type under study.
Thank you
the answers and the pdf you shared was helpful. I was wondering if it is also correct to do a rigid scan instead of a relaxed scan using QM and then also calculated U_NB(phi) using a rigid scan and then fit the difference to both the dihedral contributions?
it makes sense to compare rigid with rigid, however since the rigid scan is overestimating all the energies (excluded the minimum), you will end up with wrong barriers and relative energies of conformers. How much "wrong" depends on the molecule.
So, using a rigid scan is OK only if you are interested in reproducing the minimum energy geometry, and not if your target is to obtain the correct population of conformers and interconversion rates.
Thanks a lot, I implemented your suggestions and was able to get quite a decent fit. I have anohter question though, what if my molecule has more than one dihedrals (mine has two).
My molecule looks something like the one in the figure and the way I did it was, I fixed 3-5-6-7 and 3-5-7-8 and scanned along 1-2-3-4 to obtain the dihedral constants for the 2-3 dihedral.
Then used the constants for the 2-3 dihedral and scanned along 3-5-6-7 letting everything relax to obtain dihedral constants for the 5-6 dihedral.
In that case which two (or more) dihedrals 1 & 2 are close (their energy profiles depend one on each other), you have to repeat the procedure until it converges, so:
a) measure the MM potential for 1, U_MM1
b) correct the FF parameters for 1 for mininising U_MM1-U_QM1
c) measure the FF potential for 2, U_MM2
d) correct the FF parameters for 2 for minimising U_MM2-U_QM2
e) if { |U_MM1-U_QM1|
- Badges
- Science topic
- Similar topics
- Chemistry
- Theoretical Chemistry
- Computational Chemistry
- Molecular Dynamics