[gmx-users] need help with free energy of solvation calculati on
Moore, Jonathan (J)
JMoore2 at dow.com
Tue Jan 10 23:56:42 CET 2006
Thanks for the help and suggestions...
>I am not sure why your A and B states would be the same,
>though; your toplogy file (mostly) looks OK (see below).
>If you can send the full topology and coordinate files
>and an mdp file as attachments I can take a more detailed look.
It looks like my B states was not being defined because of the way I was defining my atoms in my .itp where I defined the molecule. I was using the form:
[ atoms ]
; nr type resnr residu atom cgnr charge typeB chargeB
4 CH1 1 H000 C3 2 0.232 CH1 0.180
with the masses defined in the [ atomtypes ] section of the ff45a4nb.itp file instead. Defining the masses in the [ atomtypes ] section worked OK for a non-FE simulation (i.e., no B-state to be defined), but not for FE. My B state definition didn't begin to work until I explicitly defined the A- and B-state masses in the [ atoms ] section like this:
[ atoms ]
; nr type resnr residu atom cgnr charge mass typeB chargeB massB
4 CH1 1 H000 C3 2 0.232 13.0190 CH1 0.180 13.0190
Also, I saw somewhere in the manual that you need to define B state parameters for angles, torsions, etc. that involve perturbed atoms even if that particular interaction doesn't change.
Now I think my A and B states are being defined properly. I am getting several warnings of this type from grommp:
"No default Proper Dih. types for perturbed atoms, using normal values"
This warning occurs each time my dihedral type gd_17 appears, for example:
[ dihedrals ]
; ai aj ak al funct type
13 3 4 7 1 gd_17 gd_17
Atoms 4 and 7 are being perturbed, so I understand that it is expecting a B-state dihedral type. However, as you can see above, I have defined the the B-state dihedral type, so I don't know why it is complaining. Anyone have any ideas? Looking in gmxdump, I think the parameters that are being assigned for those proper torsions are correct, so I guess I'll just ignore the warnings for now using the -maxwarn option and try the FE calculation again.
>Also (and here's where my ignorance of the paper comes in),
>are you actually intending to change some of the atoms from
>hydrogens to carbons, and appear new sites? Your topology
>doesn't look like it's appearing any new sites right now,
>I don't think (don't you need dummy atoms initially for that?
The polar hydrogens that I'm changing are explicitly represented (charges but no vdw interactions), so I don't need dummy atoms for them.
Jonathan Moore, Ph.D.
Research and Engineering Sciences - New Products
The Dow Chemical Company
1702 Building, Office 4E
Midland, MI 48674 USA
Phone: (989) 636-9765
Fax: (989) 636-4019
E Mail: jmoore2 at dow.com
From: gmx-users-bounces at gromacs.org [mailto:gmx-users-bounces at gromacs.org] On Behalf Of David Mobley
Sent: Tuesday, January 10, 2006 10:31 AM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] need help with free energy of solvation calculation
I responded to this four days ago but my e-mail bounced because it was too long, so let me try again. Here's what I wrote before:
I can't speak too much to the details of the system you're working on, but there are a couple things I would try checking. I haven't read the paper, though...
>The reported value in Table 3 of that paper is 393 kJ per mol. I tried
>to calculate the value with Gromacs 3.3rc3 using thermodynamic
>integration with soft cores (since the hydrogens that change to methyls
>have no sigma or
>epsilon) and 20 lambda values between 0 and 1 (50 ps of equilibration
>followed by 150 ps of production at each lambda value). I have reasonable
>confidence that I have the force field implemented correctly because I
>calculated exactly the same values for the various energy terms as one of
>the paper's authors did for a configuration of an isolated cellobiose
>molecule provided by the author. However, the free energy difference that I
>calculated was the wrong sign and an order of magnitude too large, so I must
>have done something very wrong. One thing that has me perplexed is that if
>I do a gmxdump of the .tpr file, it looks like the A and B states are
>exactly the same for things like atom type, charge, and mass...suggesting
>that I'm not specifying the B state successfully. But if the A and B states
>are exactly the same, then why am I getting anything at all for dgdl? I'm
>relatively new to GROMACS and have never done a coupling-parameter-type
One thing to be aware of is that the minute you switch on soft-core, the potential energy becomes a function of lambda. Hence you can actually do no perturbation and still get nonzero dgdl at most of your lambda values. It's only the integral that's zero. Odd, I know, but that's how it works. I am not sure why your A and B states would be the same, though; your toplogy file (mostly) looks OK (see below). If you can send the full topology and coordinate files and an mdp file as attachments I can take a more detailed look.
Also, I would do this in two steps: First change the charges, then change the vdW. Then you can use linear (not soft core) scaling for the electrostatics part and only change the vdW with soft core. You can do the electrostatics part with relatively few lambda values this way, as it's often pretty smooth. Maybe this isn't necessary here since you shouldn't have any charge sites with zero vdW radii, though (which is usually the reason you want to do this in two steps) so maybe this isn't a problem.
Also (and here's where my ignorance of the paper comes in), are you actually intending to change some of the atoms from hydrogens to carbons, and appear new sites? Your topology doesn't look like it's appearing any new sites right now, I don't think (don't you need dummy atoms initially for that? I don't have any exp with those sorts of calculations). And if you are changing atom types which involve a change in mass you'll probably have a dEkin/dlambda term which is nonzero also. Since the free energy difference really requires integrating dH/dlambda you probably need to take that into account too (or maybe you're already doing so?).
I can probably help you troubleshoot a little more if you'll send topologies, etc. Otherwise, I would try starting with something which has been set up in GROMACS before and make sure you can get the right answer. I can provide a topology for something simple like disappearing methane in water (and a right answer) if that would be helpful.
gmx-users mailing list
gmx-users at gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users