Page 1 of 6 1 2 3 4 5 6 >
Topic Options
#23873 - 04/02/10 07:30 PM generating parameters for pyruvate
RD2015 Offline
Forum Member

Registered: 01/29/10
Posts: 85
Hi,

This is my first time attempting the protocol outlined on http://dogmans.umaryland.edu/~kenno/#CGenFF.

I am working on generating parameters for pyruvate. And I want to check to see if my approach is correct.

This is what I have done:
*Note: I wasn't able to follow everything in the tutorial and I grossly simplified many steps.

1. I created the rtf file by analogy.
2. I ran generate_tutorial.inp
3. I put the resultant gjf file into Gaussian 09.
4. I looked at the Hessian written in the internal coordinates and extracted the diagonal entries as my force constants.
5. I scaled the values into kcal/(mol Angstroms)
6. For the dihedral terms I divided by multiplicity squared, and I also used the minimized coordinates to get delta (by multiplying Xmin by multiplicity and subtracting 180 degrees).
7. I entered only the missing parameters into the rtf.
8. I then ran gaussian again with scrf=(pcm,solvent=water) to simulate the presence of water.
9. I used these Mulliken charges in my rtf.
10. I ran a minimization (without waters) and got a distorted structure.
11. I adjusted the charges so that the methyl group had .09 on each hydrogen and .21 on the carbon and reran the minimization (without waters).
12. There was less distortion but the structure is not perfect.


So at this point I have a few questions:
1. Will this approach work?
2. Can I just keep tweaking charges?
3. What should my goal be? Should I just try to match the minimized gaussian structure as closely as possible using the following criteria: "typically 0.03 for the bonds and 3 for the angles".

Thanks

Top
#23888 - 04/05/10 02:53 PM Re: generating parameters for pyruvate [Re: RD2015]
Kenno Offline
Forum Member

Registered: 12/19/05
Posts: 1535
Loc: Baltimore, MD
Originally Posted By: RD2015
I am working on generating parameters for pyruvate. And I want to check to see if my approach is correct.
Thank you for your interest. This looks like a worthwile endeavor. Just for your reference, here are the suggested atom types for pyruvate (it's important to start from good atom types).
Code:
ATOM C1     CG2O3
ATOM O11    OG2D2
ATOM O12    OG2D2
ATOM C2     CG2O5
ATOM O2     OG2D3
ATOM C3     CG331
ATOM H31    HGA3
ATOM H32    HGA3
ATOM H33    HGA3
Also, don't forget to define the two improper dihedrals!

Originally Posted By: RD2015

This is what I have done:
*Note: I wasn't able to follow everything in the tutorial and I grossly simplified many steps.
Can you please elaborate on which steps you were not able to follow and why? I'd like to improve the tutorial but we get precious little feedback.



Originally Posted By: RD2015
4. I looked at the Hessian written in the internal coordinates and extracted the diagonal entries as my force constants.
This generally doesn't work well because the diagonal elements of the Hessian in most cases cannot be trivially correlated with the parameters.

Originally Posted By: RD2015
6. For the dihedral terms I divided by multiplicity squared,
confused What did you divide? By what multiplicity? Why? I don't really understand what you're doing but it's unlikely to be correct. Please perform a potential energy scan around the C-C-C-O dihedral and fit a single 2-fold term with phase 180 (or 0) to it. This is the only correct way to fit this dihedral.

Originally Posted By: RD2015
and I also used the minimized coordinates to get delta (by multiplying Xmin by multiplicity and subtracting 180 degrees).
Danger, Will Robinson! How about reading the FAQ?

Originally Posted By: RD2015
8. I then ran gaussian again with scrf=(pcm,solvent=water) to simulate the presence of water.
9. I used these Mulliken charges in my rtf.
This is wrong on so many levels.
  • Using an implicit solvent to generate target data for parameterizing a force field that is predominantly meant to work with explicit solvent is inconsistent. It is like using a condensed phase force field to generate target data for parameterizing a DFT method...
  • Mulliken is already not such a great method for calculating charges (basis set dependence,...), and combining it with PCM takes the last bit of relevance out of it.
  • Straight Mulliken charges are generally not useful in force field design; force fields that use QM-based charges use an electrostatic fitting method such as RESP.
  • CHARMM is not one of those force fields. Again, read the fine FAQ.


Originally Posted By: RD2015
10. I ran a minimization (without waters) and got a distorted structure.
11. I adjusted the charges so that the methyl group had .09 on each hydrogen and .21 on the carbon and reran the minimization (without waters).
12. There was less distortion but the structure is not perfect.
Given everything you did wrong, it is hardly surprising you get a "distorted structure". Also, you probably forgot to define improper dihedrals on the carbonyl groups.

Originally Posted By: RD2015
So at this point I have a few questions:
1. Will this approach work?
2. Can I just keep tweaking charges?
3. What should my goal be? Should I just try to match -snip-
No, no and no. We didn't make up our parameterization procedure overnight on a napkin in the bar after drinking a bit to much. It is the result of decades of experience and extensive validation by many research groups. Please re-read the CGenFF paper and also look up some of the references (especially the MacKerell and the Jorgensen ones, if you don't believe me). Then, start all over again and try to follow the tutorial as closely as you can. I'll understand if you can't figure out how to use molvib on your own, and in the specific case of pyruvate, you may get away with it, but you cannot optimize the charges without water interactions, and there's no way around doing a potential energy scan for the central dihedral. You have access to Gaussian 09; all the calculations should be light enough to run on any desktop PC produced in the last 6 years.

Edit: oh yeah, make sure you use the downloadable version of the tutorial, else you won't have all the files.
http://dogmans.umaryland.edu/~kenno/cgenff/download.html

Top
#23894 - 04/05/10 06:37 PM Re: generating parameters for pyruvate [Re: Kenno]
RD2015 Offline
Forum Member

Registered: 01/29/10
Posts: 85
Thank you for your advice. I'll start over and get back to you with any new questions.

Top
#23905 - 04/06/10 08:00 PM Re: generating parameters for pyruvate [Re: RD2015]
RD2015 Offline
Forum Member

Registered: 01/29/10
Posts: 85
Hi,

I restarted using the atom types you suggested. Everything works fine but I am now stuck at the "Target data" portion of the tutorial.

Here's what I have done so far:
1. Redid atom types as you suggested above.
2a. Changed zeed to 1 C1 1 C2 1 C3
2b. Ran generate tutorial twice to find out what parameters were missing
2c. The missing parameters were:

ANGLES
CG2O3 CG2O5 CG331 35.00 116.00 ! adapted from "ACO, acetone adm 11/08"

DIHEDRALS
CG2O3 CG2O5 CG331 HGA3 0.1000 3 0.00 !adapted from "BTON, butanone
OG2D2 CG2O3 CG2O5 CG331 0.3000 2 180.00 !adapted from "BIPHENYL ANALOGS unmodified,peml"

3. Ran gaussian on pyr_opt_freq_mp2.gjf

4. Assigned .09 charge to each hydrogen. And rescaled Mulliken charges:
i) x = (Sum of Mulliken charges for hydrogens) - (# of hydrogens * .09)
ii) x/(remaining # of atoms)=y
iii) (mulliken charge on remaining atom) +y= Charge inputed into pyruvate.str file

Now I have a few questions:

1. In the tutorial you prescribe using Merz-Kollman charges. I can't find them in the output gaussian log file. Is it okay to use Mulliken charges?

2. I am really confused about what to do next. Am I supposed to generate different monohydrate pyruvates and then minimize in gaussian to obtain an energy? (If so how many do I need? 3?)

3. Assuming that's the right idea, how do I generate these files that are analogous to "prld_water_hf_H1_0.gjf"?

4. I was also looking at the next script "water_constr_tutorial.inp" and I noticed that it requires a file called @residue_mp2.crd. I realize this comes from "pyr_opt_freq_mp2" and is somewhere under "Input Orientation" in the log file. What is the standard way to generate the crd from the log file?

Thanks again

EDIT: I forgot to ask about the improper terms.

I am not sure I understand when the improper terms need to be included. Are the two impropers I need to include: C2 attached to C1 C3 O2, and C1 attached to O11, O12, C2?


Edited by RD2015 (04/06/10 09:00 PM)
Edit Reason: forgot impropers

Top
#23913 - 04/07/10 01:13 PM Re: generating parameters for pyruvate [Re: RD2015]
Kenno Offline
Forum Member

Registered: 12/19/05
Posts: 1535
Loc: Baltimore, MD
Originally Posted By: RD2015
-snip-
2c. The missing parameters were:
Code:
ANGLES
CG2O3  CG2O5  CG331    35.00    116.00               ! adapted from "ACO, acetone adm 11/08"

DIHEDRALS
CG2O3  CG2O5  CG331  HGA3       0.1000  3     0.00   !adapted from "BTON, butanone
OG2D2  CG2O3  CG2O5  CG331      0.3000  2   180.00   !adapted from "BIPHENYL ANALOGS unmodified,peml"
These look like good initial guesses to me.

Originally Posted By: RD2015
1. In the tutorial you prescribe using Merz-Kollman charges. I can't find them in the output gaussian log file. Is it okay to use Mulliken charges?
Nope, still not OK. Read what I wrote about Mulliken charges in my previous post. We haven't changed our mind overnight.

I don't know about Gaussian 09, but Gaussian 03 produces an output section titled "Electrostatic Properties Using The MP2 Density" which contains (among other things) the Merz-Kollman charges, under "Charges from ESP fit" (if you used the input from the tutorial). Be careful, it does this 2 times, once for the input geometry and a second time for the minimized geometry. It's the second one you want to look at!

Also note that we use the Merz-Kollman charges only as an initial guess. If you're not going to generate water interactions and optimize the charges based on those, then forget about QM-based charges altogether and just use charges by analogy to existing compounds (but then again, the analogy will not be very good). Again, see this FAQ.

Originally Posted By: RD2015
4. Assigned .09 charge to each hydrogen. And rescaled Mulliken charges:
i) x = (Sum of Mulliken charges for hydrogens) - (# of hydrogens * .09)
ii) x/(remaining # of atoms)=y
iii) (mulliken charge on remaining atom) +y= Charge inputed into pyruvate.str file
Holy contrivance, batman! I have to repeat that you can't just make up your own scheme for assigning charges. Also, what is wrong with the + and - operators? Scaling charges is not a good idea because you're also scaling the dipole moment and all polar interactions the molecule would make in condensed phase. Here's what you should do instead (using the Merz-Kollman charges):
  • Sum the 4 charges on the CH3 group together.
  • Subtract (3 x .09 = .27) to get the charge on the carbon.
  • Assign .09 to each hydrogen.
In other words: the charge on the whole CH3 group should remain the same. Easy, isn't it?

Originally Posted By: RD2015
2. I am really confused about what to do next. Am I supposed to generate different monohydrate pyruvates and then minimize in gaussian to obtain an energy? (If so how many do I need? 3?)
The more water interactions you have, the less your system is underdefined and the easier it gets to generate relevant charges. Look at figure 2 in the CGenFF paper: we put waters wherever we saw possibility for a hydrogen-bond-like interaction. However, for pyruvate, it's slightly more complicated because you have carbonyl groups. For carbonyl groups, we usually generate two water interactions per oxygen: one where the oxygen donates a hydrogen bond in the lone pair direction and one where the oxygen donates a hydrogen bond in the direction of the C=O bond (ie. linear). Out of the two lone pairs each oxygen has, you should choose the one that suffers the least contamination by other parts of the molecule. For example, on the C2=O carbonyl, you should choose the lone pair that points away from the carboxylate. In total, you would have about 8 water interactions (2 per oxygen + 2 for the methyl hydrogens, assuming that the molecule is planar so that the 3rd methyl hydrogen is equivalent to one of the others for symmetry reasons).

Originally Posted By: RD2015
3. Assuming that's the right idea, how do I generate these files that are analogous to "prld_water_hf_H1_0.gjf"?
Since you're the first one asking this question on the forums and it's not a totally trivial thing to do, I'm willing to generate the input files for you so that other people can use them as an example. So if you can't figure it out on yourself, just use the forum's "File Manager" feature to attach your toppar file and your gaussian output to your next post, and I'll take care of it when I find a little hole in my schedule.

Originally Posted By: RD2015
4. I was also looking at the next script "water_constr_tutorial.inp" and I noticed that it requires a file called @residue_mp2.crd. I realize this comes from "pyr_opt_freq_mp2" and is somewhere under "Input Orientation" in the log file. What is the standard way to generate the crd from the log file?
The standard way is to take an existing crd file for the same molecule and copy-pasting the cartesian coordinates from the gaussian log file into it. Beware: the fields in the crd format are fixed-width. I have a script that does this but I may have forgotten to include it in the "downloadable" tutorial.

Note that because your molecule has a net charge, you should not apply the 1.16 scaling factors to your interaction energies, and also, you cannot compare the QM and CGenFF dipole moments, as the dipole moment is ill-defined for a charged entity. We do have a procedure to align the QM and MM geometry to the origin and to each other in a way that yields a reasonable comparison for most molecules, but this is not documented in the paper nor in the tutorial. I'm busy right now, but if you're interested, I'll try to post the CHARMM script one of the coming days.

Originally Posted By: RD2015
I am not sure I understand when the improper terms need to be included. Are the two impropers I need to include: C2 attached to C1 C3 O2, and C1 attached to O11, O12, C2?
That's correct.

Top
#23940 - 04/08/10 04:05 PM Re: generating parameters for pyruvate [Re: Kenno]
RD2015 Offline
Forum Member

Registered: 01/29/10
Posts: 85
Here are my files. Thanks for your help.


Attachments
params-RD2015.tar.gz (338 downloads)


Top
#24042 - 04/16/10 05:55 PM Re: generating parameters for pyruvate [Re: RD2015]
Kenno Offline
Forum Member

Registered: 12/19/05
Posts: 1535
Loc: Baltimore, MD
Sorry for my late reply, it slipped a bit off my radar. Attached is your toppar stream file, after the following corrections:
  • completed your IC table
  • impropers must be defined both in the topology and parameter section
  • we're going to abolish wildcard impropers in a future CGenFF release so I changed the wildcard atoms into explicit ones
  • when editing machine-readable ASCII files, please use a font in which you can differentiate between O and 0... corrected several of these mistakes.
Also attached are the Gaussian input files for the water interactions. Note that there's some redundancy: you'll have to make a choice between O2_0 and O2_180 based on which one suffers the least contamination. This may become obvious once you have the interaction distances and energies. Same for O12_0 and O12_180 . Get back to me when you have a list of HF water interaction distances and energies, then I'll show you the next step (reproducing them with CHARMM).

Oh yeah, there's also a file named pyru_center.gjf . You can run it if you want to have dipole moments that are (somewhat) comparable with the CHARMM ones.


Attachments
pyru_water_gauss.tgz (366 downloads)


Top
#24050 - 04/18/10 11:35 PM Re: generating parameters for pyruvate [Re: Kenno]
RD2015 Offline
Forum Member

Registered: 01/29/10
Posts: 85
Thanks again for the help.

I had a problem running g09 on pyru_water_hf_H31.gjf. The run terminates via linke 9999 error. I tried increasing resources and also put in the calcall option. And it still crashes.

I've attached the files.


Attachments
pyru_water.tar.gz (374 downloads)


Top
#24055 - 04/19/10 01:44 PM Re: generating parameters for pyruvate [Re: RD2015]
Kenno Offline
Forum Member

Registered: 12/19/05
Posts: 1535
Loc: Baltimore, MD
What happened is so obvious I'm not going to tell you. Carefully look at your output. Shotgun debugging almost never works!

Top
#24069 - 04/21/10 01:02 PM Re: generating parameters for pyruvate [Re: Kenno]
Kenno Offline
Forum Member

Registered: 12/19/05
Posts: 1535
Loc: Baltimore, MD
This doesn't bode well for your ability to bring this parameter optimization to a good end, but here's a hint. Your output contains the following:
Code:
 Optimization stopped.
    -- Number of steps exceeded,  NStep=  20
This is indeed the most common cause of Gaussian's "Error termination request processed by link 9999." Your geometry optimization ran out of cycles. Why did it run out of cycles? Look at the geometry. Or, in your case, since you're optimizing only one degree of freedom, you could just look at the value of that degree of freedom in your Gaussian output. What do you see? Why do you think this would happen?


Edited by Kenno (04/21/10 03:09 PM)
Edit Reason: spelling error

Top
Page 1 of 6 1 2 3 4 5 6 >

Moderator:  alex, chmgr, John Legato