[gmx-users] soft-core and coulomb transformation

[David Mobley]

Hi, Bharat, and David vdS.

This indeed is looking suspiciously like a real problem. FWIW, I did
reproduce your original problem using my own input files and PME with
gromacs 3.3.1.

It's still a bit perplexing, in that you seem to be getting basically
the "right" answers for the forward transformation (at least, by
comparing the PME forward results with the cutoff forward results) but
the reverse transformation is crazily wrong. I'm having a hard time
figuring out what sort of bug could cause this to happen.

The only thing I can come up with so far that might cause this is a
bug in the handling of B state 1-4 interactions with PME. Since
methane has no B state 1-4 electrostatic interactions, the
ethane->methane transformation would be OK both with and without PME.
But methane->ethane would go awry with PME because ethane still has
1-4 electrostatic interactions.

If I remember correctly, though, you reported that the vacuum
calculations work fine? That wouldn't be consistent with my
explanation above, and perhaps would suggest that the problem is
interactions with water, somehow?

Just speculation at this point.

I'm doing some more testing on this end (trying to come up with a more
minimal example). I'm checking to see whether discharging phenol in
water gives the same results as charging phenol in water.

David


On 11/11/07, bharat v. adkar <bharat at sscu.iisc.ernet.in> wrote:
>
> Hi David,
>    I am keeping the subject same so that it is helpful to track in
> future.... :)
>
>
> On Sat, 10 Nov 2007, David Mobley wrote:
> > Dear Bharat,
> >
> > OK, I went ahead and ran with your topologies using my own run
> > scripts. My runs haven't finished yet,  but looking at just
> > <dV/dlambda> values from equilibration,  I seem to be seeing roughly
> > the same trends you do. In particular:
> > (1) The forward simulation always gives dv/dlambda values that are
> > fairly close to -2 at each lambda
> > (2) The reverse simulation gives a range of <dv/dlambda> values at
> > different lambda values, and switches from positive to negative as
> > lambda increases.
> >
> > The integral of <dV/dlambda> in the forward case is NOT the negative
> > of the integral of <dV/dlambda> for the reverse case, which is what it
> > should be.
> >
> > Obviously, something is wrong here, and it doesn't seem to be your run
> > input files. I don't see anything obviously wrong with your topology
> > files, either -- as far as I can tell your charge groups are OK.
> >
> > My suggestion is to try to strip this down to a simpler case where it
> > will become more clear what the problem is. Maybe ethane->methane
> > rather than with the capping on the ends.  (Did you say you have the
> > same problem with ethane->methane? It would be easier to troubleshoot
> > those topologies).
>
> As I mentioned earlier, ethane to methane transformation also shows the
> same behaviour. I am almost sure that the topologies and other input
> files are correct. I am pasting below the output of the
> coulomb-transformation for both, ethane to methane (forward) and methane
> to ethane (reverse). Topology is the same as i mentioned few mails
> back...
>
> Forward charge transformation
> lmbd  dg/dl    error
> 0.00 -8.286304 0.0108319
> 0.05 -8.32708  0.0114018
> 0.10 -8.328157 0.0095819
> 0.15 -8.344236 0.00893341
> 0.25 -8.353463 0.00984712
> 0.40 -8.368163 0.0104302
> 0.50 -8.397035 0.00969758
> 0.60 -8.415603 0.00908416
> 0.75 -8.448185 0.00941215
> 0.85 -8.484247 0.00922024
> 0.90 -8.481465 0.00885981
> 0.95 -8.495802 0.00961897
> 1.00 -8.501283 0.00957299
>
> Reverse charge transformation
> 0.00 13.770658 0.074592
> 0.05 13.331217 0.0778554
> 0.10 13.075527 0.0722786
> 0.15 12.483814 0.0750718
> 0.25 11.789282 0.0674654
> 0.40 10.730386 0.0875609
> 0.50 10.008684 0.0807128
> 0.60  9.266026 0.0736956
> 0.75  7.999017 0.0738001
> 0.85  7.388243 0.0832792
> 0.90  6.991334 0.0828869
> 0.95  6.595921 0.0796503
> 1.00  6.180366 0.0893983
>
> As can be clearly seen, both transformations are not equivalent. (These
> are output from the 1 ns simulations.)
>
>
> >
> > If you get to the point where you're convinced it has nothing to do
> > with your topology, you could submit a bugzilla -- but it would be a
> > good idea to be more sure that there is no problem with your topology,
> > first. Also it would be helpful to have more evidence about where the
> > bug might be, if there is one (hence the importance of stripping this
> > down to the "minimum" topology necessary to reproduce the problem).
> >
> > You might also, for example, just try turning off the charges on
> > methane in water. You could also try turning on the charges in methane
> > in water. These two should be equivalent (except for the sign) of
> > course.
> >
> > Keep me posted on what you find out.
>
> to further see what happens, i used cut-off instead of pme and that seems
> to be giving okay results. i tried running only 20 ps runs as in
> equilibration, and there the values look as per expectations. then tried
> with Ace-Ala-Nac <-> Ace-Gly-Nac also. both seems okay. below are the
> values i obtained:
>
> Ethane <-> Methane
> forward
> 0.00 -8.277302 0.0747561
> 0.05 -8.346413 0.105338
> 0.10 -8.311737 0.0775508
> 0.15 -8.326182 0.0749244
> 0.25 -8.383022 0.0741275
> 0.40 -8.309543 0.199901
> 0.50 -8.388351 0.0807475
> 0.60 -8.334778 0.0769525
> 0.75 -8.444722 0.0556352
> 0.85 -8.541355 0.0529021
> 0.90 -8.562161 0.079475
> 0.95 -8.588182 0.0679287
> 1.00 -8.535591 0.0505447
>
> reverse
> 0.00 8.647486 0.0848346
> 0.05 8.571583 0.0678888
> 0.10 8.410289 0.04688
> 0.15 8.336588 0.085721
> 0.25 8.536206 0.0457025
> 0.40 8.452587 0.0730119
> 0.50 8.397571 0.0540216
> 0.60 8.379067 0.0775835
> 0.75 8.345908 0.0889559
> 0.85 8.408099 0.082048
> 0.90 8.284268 0.0632315
> 0.95 8.1428   0.123793
> 1.00 8.303987 0.0505205
>
>
> and for Ace-Ala-Nac <-> Ace-Gly-Nac
> forward
> 0.00 -2.113457 0.160573
> 0.05 -2.056746 0.0790143
> 0.10 -2.619748 0.242052
> 0.15 -2.470358 0.247365
> 0.25 -2.121426 0.623961
> 0.35 -2.210045 0.0567349
> 0.45 -2.173477 0.0483762
> 0.50 -2.156053 0.167631
> 0.55 -2.156565 0.0948302
> 0.65 -2.194102 0.0810602
> 0.75 -2.302156 0.306185
> 0.85 -2.28119  0.0545011
> 0.90 -1.936869 0.147152
> 0.95 -2.275215 0.190621
> 1.00 -2.221756 0.124669
>
> reverse
> 0.00 2.30627  0.169452
> 0.05 2.102584 0.181941
> 0.10 2.531891 0.27292
> 0.15 2.491334 0.200524
> 0.25 2.020721 0.0577124
> 0.35 2.42911  0.133624
> 0.45 2.113785 0.24776
> 0.50 2.312891 0.056794
> 0.55 2.147666 0.0652882
> 0.65 2.259201 0.0969329
> 0.75 2.341642 0.0767809
> 0.85 2.381157 0.19062
> 0.90 2.293333 0.238674
> 0.95 2.197692 0.175085
> 1.00 1.921228 0.0853432
>
>
> in both the cases, the values do not change much, but they appear to be
> satisfying.
>
> So everythings approximately zeroes down to the suspision that PME
> implementation for TI is buggy.
>
> Developers, please have a look at and confirm.
>
>
> bharat
>
>
>
>
>
> >
> > David
> >
> >
> > Hope that helps. Keep me posted on what you find.
> >
> > On Nov 7, 2007 10:42 PM, bharat v. adkar <bharat at sscu.iisc.ernet.in> wrote:
> >>
> >> Hi David,
> >>
> >> This is regarding our discussion on gromacs mailing list on "soft-core
> >> and coulomb transformation". Sorry for delay in reply. I am  attaching a
> >> TGZ file which contains a script to run the simulations and a MDP folder
> >> with all required files.
> >>
> >> The system is Ace-Ala-Nac, i.e., alanine capped at both ends.
> >>
> >> The MDP files are made to run the FORWARD simulations, i.e., from
> >> alanine->glycine. To run the REVERSE simulations, in all MDP files,
> >> "define = -DFE_coul_rev" has to be defined.
> >>
> >> I am sorry for a lot of botherance. but i am not getting the same values,
> >> within error limit, as i mentioned in the mailing list. I will send my set
> >> of values to you most probably by tomorrow.
> >>
> >> Thanks.
> >>
> >> bharat
> >>
> >> --
> >> This message has been scanned for viruses and
> >> dangerous content by MailScanner, and is
> >> believed to be clean.
> >>
> >>
> >
> >
>
> --
> This message has been scanned for viruses and
> dangerous content by MailScanner, and is
> believed to be clean.
>
>