The phase angle for dihedrals is almost always 0 or 180 in most force fields. In concert with the multiplicity and force constant/amplitude, you will get the shape of the periodic curve.

Hi Justin, 

Thanks for the answer. After looking through Amber's parm99.dat file, I figured it will be one or the other. I'm still yet to figure out which do I pick for my specific arrangement of unknown bond connections. Next on the list is deciding equilibrium bond lengths and bond angles. 

An initial guess by analogy is reasonable, but dihedrals have very poor transferability.

For angles and dihedrals, QM optimizations and comparison with any available crystal data would be suitable.

I do believe I took the approach you mentioned ("For angles and dihedrals, QM optimizations and comparison with any available crystal data would be suitable."). 

Would you be so helpful to just clarify something for me?

Using an angle between X-Y-Z atoms, I assume the angle in any forcefield is the "equilibrium angle", i.e., any difference between the actual angle and the equilibirum angle will incur a potential energy penalty. Would the resulting X-Y-Z angle of the ab initio optimised structure be the "equilibrium angle" for this particular unknown entry in the MD forcefield?

I thought of taking the measured bond angle and bond length from my ab initio calculation and using them in the force field. For force constants, I would look for an existing force constant in the force field that satisfies both (a) identical atom elements (not types) and in the same order, and (b) the angle/distance is the same. After all, surely the force constant is a relationship between bonded element types and their distance/angle.

Thanks

Yes, the equilibrium angle is either the angle in the optimized geometry, or in an available crystal structure, or some compromise between the two. The force constant should be based on the vibrational frequency (easily done in QM). Analogy is probably OK, but if you really want to go about it rigorously, it's not hard to do.

Thanks Justin. These 'simple' things are sometimes as clear as mud until someone with more experience explains it. 

For the angles, I have used the measure QM angles for my angle equilibrium values in for the unknown parameters in my force field. I've yet to discovered how the vibrational frequency output in Gaussian translate into a force constant. That's next on the list. At the moment I am assuming a angle/distance with atom type relationship and using an analogy. 

Thanks

http://www.gaussian.com/g_tech/g_ur/k_freq.htm

Thanks Justin. I've read that several times, but funnily enough, it wasn't until tonight that I finally noticed the "frc consts" entry above the labeled vibrational mode!  

Hi Justin, I would appreciate your thoughts on this one. The force constants (mDyn/A) returned by Gaussian are of the vibrational normal modes (VNM). From the looks of it, there is not a one-to-one relationship between VNM and a particular atomic motion e.g., stretch, twisting or torsional twist. A VNM may contain multiple atomic motions. Perhaps I am wrong, but from the looks of it I am going to need to perform a long list of potential energy surface scan in Gaussian and then fit the potential energy to a particular model to derive force constants. 

Does this sound about right?

Unless you're using a program like CHARMM that can read the Gaussian frequency output and assign each of the actual bonded vibrational states, it will can be difficult to map everything exactly.  CHARMM's MOLVIB does this, I don't know what else might have this capability.  For bonds and angles, a straightforward method is 3-point scans (in which you fix the value of the target bond length or angle at some value that deviates from the equilibrium value, e.g. +/- 5 or 10 degrees for angles) and then optimize the rest, calculate the energy, and compare with your force field.  Of course, this is all (by necessity) iterative as you refine other parameters, since the energy will change. Welcome to the fun parts of force field parametrization :)

I am very new (well, a year in but not exactly dedicated to it) to force field optimisation. My PhD was on rational-designed TM helical proteins and free energy calculations. Back then it was a case of being plagued by the long correlation times of membrane systems. Now it's something else. 

Anyway, I think I have found a reasonably good way of doing it, but as you said, it will be difficult. I have found a method by Seminario et al., FREQ my structure which was returned from RED. I then need to open MATLAB and start figuring out how I can isolate a 3x3 hessian matrix of the atom of interest from the whole system Hessian. Then, diagonalise, get the eignen vectors and then derive the force constant from the eigenvalue (plus other maths stuff). But who knows! This might be completely wrong. 

Time to carry on digging. Thanks :) 

http://alma.karlov.mff.cuni.cz/bio/99_Studenti/00_Dalsi/ParamFit/2013_ParamFit_AmberTools13.pdf