Hi Anthony,

if I was to amend the Amber force field by adding new interaction potentials I would not generate the missing parameters as described in the above paper. I would go the "standard" Amber way, computing the missing interaction potentials on a small model compound. 

I give You a short example. Suppose You want to simulate a compound containing complexed Magnesium in your compound an you miss the following interaction potentials (atom types from the Amber gaff force field):

• c3-n4-mg
• n4-mg-n4
• n4-mg-os
• n4-mg-ca
• mg-ca-ca
• mg-os-c3
• ca-mg-os

Then You construct the smallest structure containing all these interactions (Mg / P1, fig01.png). The atom masses and Van-der-Waals radii you take from another force field (e.g. parm99.dat). Then  you optimize the geometry using gaussian, as you are planning,  with:

# HF 6-31G* Test opt

Afterwards you scan the interaction potentials (bonding distances and angles) around the equilibrium bond lengths with Gaussian.

# HF 6-31G* Test scan

.... stuff about structure ...

angle_c3-n4-mg    104.0     scan 20     1.0

Find attached an output file of this gaussian scan job. I varried the angles with 20 steps a 1 degree around the eq. angle and bond lengths with 40 steps a 0.01 angstroms around the eq. bond distance. You will get a harmonic interaction potentail (see fig02.png) when you plot the energies. Fit this plot to a harmonic potential

V = 0.5 k (R - R0) + V0

and you get a harmonic force constant k. As you get Hartree energies from gaussian you have to convert to your foce fields units. For Amber i think you need (but please recheck yourself):

Hartree / angstroms2 = 627.6 kcal / (mol * angstroms2)

Hartree / (1°)2 = 2.06 * 106 kcal / (mol * rad2)

This way you can get the missing parameters more in accord to the amber force field definition. The way they do it in the paper above I would only derive parameters to use them exclusively for your compound and not mix it with another force field. 

Good luck and regards

Sebastian

Hi Anthony,

it has been some time, that I have ussed Gaussian and I have no hands on the software to date, but as I remember the force constants to any vibration are given in the standard output of a FREQ calculation. Somewhere close to the reported normal modes and the reduced masses and .... .

Or are you asking of how to assign force constants to individual bonds in your molecule, for example to derive parameters for to use a classical force field? You might be interested in the publication I attached. But if you only need to parameterize parts of a froce  field I would go for a simple computation on a small compound containing your structure of interest. 

Regards Sebastian

Article Determination of precise harmonic force constants for alanin...

"Or are you asking of how to assign force constants to individual bonds in your molecule, for example to derive parameters for to use a classical force field? You might be interested in the publication I attached. "

I've went through the attached publication. It was very interesting, I will see if I can get my hands on a copy of the software to generate the force constants from the internal coordinate hessian. 

So I am thinking. I want to include a new compound into an Amber force field (but using gromacs). I understand that Amber was parameterised using HF/6-31G* (I think). Should I optimise my new structure using this ab initio model chemistry, extract the force constants using software such as FUERZA, and then plug the values in to the Amber force field?

Thanks for the help. 

Hi Anthony,

if I was to amend the Amber force field by adding new interaction potentials I would not generate the missing parameters as described in the above paper. I would go the "standard" Amber way, computing the missing interaction potentials on a small model compound. 

I give You a short example. Suppose You want to simulate a compound containing complexed Magnesium in your compound an you miss the following interaction potentials (atom types from the Amber gaff force field):

• c3-n4-mg
• n4-mg-n4
• n4-mg-os
• n4-mg-ca
• mg-ca-ca
• mg-os-c3
• ca-mg-os

Then You construct the smallest structure containing all these interactions (Mg / P1, fig01.png). The atom masses and Van-der-Waals radii you take from another force field (e.g. parm99.dat). Then  you optimize the geometry using gaussian, as you are planning,  with:

# HF 6-31G* Test opt

Afterwards you scan the interaction potentials (bonding distances and angles) around the equilibrium bond lengths with Gaussian.

# HF 6-31G* Test scan

.... stuff about structure ...

angle_c3-n4-mg    104.0     scan 20     1.0

Find attached an output file of this gaussian scan job. I varried the angles with 20 steps a 1 degree around the eq. angle and bond lengths with 40 steps a 0.01 angstroms around the eq. bond distance. You will get a harmonic interaction potentail (see fig02.png) when you plot the energies. Fit this plot to a harmonic potential

V = 0.5 k (R - R0) + V0

and you get a harmonic force constant k. As you get Hartree energies from gaussian you have to convert to your foce fields units. For Amber i think you need (but please recheck yourself):

Hartree / angstroms2 = 627.6 kcal / (mol * angstroms2)

Hartree / (1°)2 = 2.06 * 106 kcal / (mol * rad2)

This way you can get the missing parameters more in accord to the amber force field definition. The way they do it in the paper above I would only derive parameters to use them exclusively for your compound and not mix it with another force field. 

Good luck and regards

Sebastian

Hi Sebastian,

That's a very comprehensive set of instruction. Thank you very much for that. It was clear and should definitely help. I think I only have one remaining question. The compound I am trying to structure in my MD simulation is a derivative of glucose. Amber does not include it (neither do any of the other force fields). Given that it only includes N, C, O and H, is it possible to 'stitch' the force field together using existing parameters (pulling apart amino acids as an example) or is it a 100% correct way to totally reparamterise a new compound from scratch?

Hi Anthony,

somehow I missed Your last question. Please excuse my delayed answers. I would choose a force-filed which is reported to give good results for systems similar to yours in the literature and I would target the main (most difficult) properties which are important to answer your research question. 

Then I would go and parametrize the only missing interactions in the force field! Maybe I was not clear in my answer above. Take the molecule in figure one in my previous post. I parametrized this molecule using the gaff force field from Amber which is intended for "small"  molecules to be used together with the amber force-field for the protein etc. I took all parameters from this (gaff) force field and only any missing parameter (above all Mg containing interactions) I parameterized as described above. The Electrostatics (=Partial Charges) for this compound should then be generated in accord with the force-field you use. I attached a link to the Amber tutorials, you will find a lot of information there (even if you do not use the Amber software). The AmberTools package (see link) is free of charge and you might be able to use some of these tools to prepare your stuff. I do not work any more with MD and have no idea about gromacs, but this software might provide similar preparation tools.

Fot the Amber/gaff force-filed the Partial Charges are derived via a RESP fit to the Electrostatic potential, e.g from a gaussian simulation:

# HF 6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) opt

The fit is done with Antechamber from the AmberTools. If you want it less complicated You could just take the Mulliken-Charges from the Gaussian simulation. Look at the Amber tutorial, you will find a lot of informations!

The I would run some test simulations and check the behaviour of the new compound. For a glucose derivative check the ring structure/dynamics. If then Your simulations show strange behavious, .... Then I would think  about reparametrizing everything from scratch for this compound. But I would really try to avoid a mixing of force-filelds in a simulation (exept those intended to be mixed e.g. Amber99SB and gaff).

Best Regards

Sebastian

http://ambermd.org/tutorials/

http://ambermd.org/#AmberTools

Mr. Nash,

They are directly obtained as values, using "freq" or "Freq=Raman" keywords. An example:

115 116 117

?A ?A ?A

Frequencies -- 5546.4817 5548.4316 5560.1790

Red. masses -- 1.0879 1.0879 1.0509

Frc consts -- 19.7188 19.7319 19.1423

IR Inten -- 12.0701 9.6398 352.3205

Raman Activ -- 14.4955 14.3278 129.1547

Depolar -- 0.0587 0.0552 0.5899

Atom AN X Y Z X Y Z X Y Z

1 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

2 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

3 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

4 6 0.00 -0.07 0.05 0.00 0.00 0.00 0.00 0.00 0.00

5 6 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00

6 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

7 7 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

..........................................................................................................................................................................................................................

29 7 0.00 0.00 0.00 0.00 0.00 0.00 -0.06 0.00 0.00

30 7 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

31 8 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

118 119 120

?A ?A ?A

Frequencies -- 5584.0741 5586.0175 5608.1146

Red. masses -- 1.0874 1.0871 1.0467

Frc consts -- 19.9767 19.9865 19.3954

IR Inten -- 14.3847 4.0979 76.5424

Raman Activ -- 15.8794 8.2938 34.7180

Depolar -- 0.7145 0.4025 0.0791

Atom AN X Y Z X Y Z X Y Z

1 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

2 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

3 6 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

.................................................................................................................................................................................................................................

40 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

41 1 0.00 0.02 -0.04 0.00 0.00 0.00 0.00 0.00 0.00

42 1 -0.35 -0.16 0.09 0.00 0.00 0.00 0.00 0.00 0.00

43 1 0.08 -0.76 0.51 0.00 -0.01 0.00 0.00 0.00 0.00

44 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

124 125 126

?A ?A ?A

Frequencies -- 7195.8766 7877.3162 7882.2722

Red. masses -- 1.0866 1.0756 1.0757

Frc consts -- 33.1509 39.3253 39.3769

IR Inten -- 1.9118 71.4026 75.2443

Raman Activ -- 20.2606 8.5695 7.8491

Depolar -- 0.2007 0.3908 0.2089

Atom AN X Y Z X Y Z X Y Z