2c. The missing parameters were:
CG2O3 CG2O5 CG331 35.00 116.00 ! adapted from "ACO, acetone adm 11/08"
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.
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
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?
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).
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.
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.
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?