[gmx-users] need help with free energy of solvation calculati on

[Moore, Jonathan (J)]

David M., 

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.

Thanks,
Jonathan

____________________________
Jonathan Moore, Ph.D.
Research and Engineering Sciences - New Products
Core R&D
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


-----Original Message-----
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


Jonathan,

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:

Jonathan,

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
>calculation before.

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.

David
_______________________________________________
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.