[gmx-users] fitting dihedral angle function parameters using gromacs

[chris.neale at utoronto.ca]

I have never done this myself. Nevertheless, I'm going to take a shot  
since you haven't got an answer yet.

Dihedral parameters can be thought of as consistently developed fudge  
factors that are used to fine-tune a parameter set. The word  
consistent here implies that there is a method to developing them that  
is specific to each forcefield and you should follow that method.  
Therefore everything that I say here is only useful in terms of  
helping you to get an idea of what these 3 steps might look like and  
not necessarily what they should look like in your particular case.

I call dihedrals fudge factors since they don't imply anything of any  
physical relevance in and of themselves when you have explicit  
hydrogens. For example, the emphasis that they may place on staggered  
conformations of sp3 carbons might be able to come out of the NB LJ  
parameters directly if those parameters weren't representing some type  
of convoluted average over different interaction types (which they  
likely are). For united atom ff's the dihedrals likely do represent  
something of physical relevance... the hydrogens.

Step 2 is easy. Simply set the dihedral parameters to zero and run  
some MD simulation. I suspect that you can do this in vacuum, but  
again that depends on your forcefield. The ab initio calculations are  
done in vacuum so I don't see why this step would need any solvent.

Step 3 would start with calculating the probability of dihedral angle  
value occurrence as a function of dihedral angle value by  
histogramming your MD sampling. You can then turn this into a free  
energy by E=-KbT*ln(P). You then look back at your results from step 1  
and figure out the dihedral function that should be applied on top of  
the free energy profile that you get in step 2 in order to (we hope)  
give you the conformational preference that you desire out of your MD  
simulations. This of course would require that you did ab initio  
calculations for multiple dihedral rotamers so that you know the  
relative energies (and therefore probabilites of occurrence) of each.

Note that different force fields put emphasis on different things. For  
example, I believe that the CHARMM ff puts as it's number one priority  
the ability to reproduce actual experimental values. Therefore in this  
case, you would be also searching the literature for rotameric  
preferences and figuring out the applied dihedral potential that,  
applied over top of your explicit solvent simulation, would give you  
the correct distribution of dihedral angles.

As for your mdp file for step 2, I'd think that you want it to be  
exactly as you will eventually use for simulation with respect to the  
cutoffs.

Let me note once again that this is just what comes to my mind. I  
haven't done it and ff development is just about the most difficult  
topic in our field.

Chris.

-- original message --

Dear Gromacs developer:

I am asking for your help about fitting a dihedral angle parameters using
gromacs md simulation. One article (link:
http://www.mpip-mainz.mpg.de/~andrienk/journal_club/opls.pdf) says that
fitting a dihedral angle potential parameters needs three steps: 1)
perform ab initio calculation; 2) do md simulation by setting the dihedral
angle parameters to zero; 3) fit all the parameters by comparing the md
energies and ab initio energies.

I am now working on a dihedral angle fitting and I have done step 1. But I
don't know how to do step 2 and 3. For example, if I want to fit a
CT-CT-CT-CT dihedral angle for butane, how do I prepare for my *.mdp file,
and after md simulation, how to fit the md energies to the ab initio
energies?
Would you please help me?

Thank you.