Discussion:
[gmx-users] Strange dgdl-value together with lincs
Maik Goette
2008-02-14 17:10:59 UTC
Permalink
Hi

While trying to find, whats wrong with some of my free energy
calculation, I stumbled upon a strange thing, I don't understand.

I'm simulating the interconversion of ethane to methanol in vacuum.
I'm splitting the process into 3 steps:
1. charge of perturbed atoms -> 0
2. morph of LJ-parameters and bonds
3. charge of perturbed atoms -> B-state values

Each system is morphed with 1 fs timestep and within 50 steps.
The total energies at the "borders" of the simulations (qqoff/vdw,
vdw/qqon, qqon/vdw and vdw/qqoff) perfectly match and show no jump; the
following simulation always got energies and velocities from the one before.

Now, when running the LJ morph A->B with lincs, the dgdl looks like this:
Loading Image...
For B->A:
Loading Image...
Doing both without lincs yields:
Loading Image...
For B->A
Loading Image...

I find this spike in the first step extremely strange. Additionally it
only occurs, when using shake or lincs as constraint solver. I guess, in
a longer process it will contribute to the integral in a neglectable
manner, but maybe theres some bug. The latest CVS shows the same behaviour.

Any comments on that?

Regards
--
Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Berk Hess
2008-02-14 17:18:19 UTC
Permalink
I assume your are using a united atom model.
What are the bonded parameters for the methanol H
in the ethane topology?
Does it have the normal methanol bond and angle potential?
If these parameters are zero, I would indeed expect a spike.

Berk.
Date: Thu, 14 Feb 2008 18:10:59 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: [gmx-users] Strange dgdl-value together with lincs
Hi
While trying to find, whats wrong with some of my free energy
calculation, I stumbled upon a strange thing, I don't understand.
I'm simulating the interconversion of ethane to methanol in vacuum.
1. charge of perturbed atoms -> 0
2. morph of LJ-parameters and bonds
3. charge of perturbed atoms -> B-state values
Each system is morphed with 1 fs timestep and within 50 steps.
The total energies at the "borders" of the simulations (qqoff/vdw,
vdw/qqon, qqon/vdw and vdw/qqoff) perfectly match and show no jump; the
following simulation always got energies and velocities from the one before.
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_lincs_01.png
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_lincs_10.png
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_nolincs_01.png
For B->A
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_nolincs_10.png
I find this spike in the first step extremely strange. Additionally it
only occurs, when using shake or lincs as constraint solver. I guess, in
a longer process it will contribute to the integral in a neglectable
manner, but maybe theres some bug. The latest CVS shows the same behaviour.
Any comments on that?
Regards
--
Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
Maik Goette
2008-02-15 08:27:24 UTC
Permalink
Berk,

here comes the topology for this run (vdw morph):
Its OPLSAA and the hydrogens have bond terms. This is the standard
OPLSAA FF for ethane and methanol.
Thanks.

[ atoms ]
; nr type resnr residue atom cgnr charge mass
typeB chargeB massB
1 opls_135 1 E2M CB 1 0.0 12.011
opls_157 0.0 12.011
2 opls_140 1 E2M HB1 1 0.0 1.008
opls_156 0.0 1.008
3 opls_140 1 E2M HB2 1 0.0 1.008
opls_156 0.0 1.008
4 opls_140 1 E2M HB3 1 0.0 1.008
opls_156 0.0 1.008
5 opls_135 1 E2M CA 2 0.0 12.011
opls_154 0.0 15.9994
6 opls_140 1 E2M HA1 2 0.0 1.008
opls_155 0.0 1.008
7 opls_140 1 E2M HA2 2 0.0 1.008
DUM 0.0 1.008
8 opls_140 1 E2M HA3 2 0.0 1.008
DUM 0.0 1.008

[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1 0.10900 284512.0 0.10900 284512.0
1 3 1 0.10900 284512.0 0.10900 284512.0
1 4 1 0.10900 284512.0 0.10900 284512.0
1 5 1 0.15290 224262.4 0.14100 267776.0
5 6 1 0.10900 284512.0 0.09450 462750.4
5 7 1 0.10900 284512.0 0.10900 284512.0
5 8 1 0.10900 284512.0 0.10900 284512.0

[ pairs ]
; ai aj funct c0 c1 c2 c3
2 6 1
2 7 1
2 8 1
3 6 1
3 7 1
3 8 1
4 6 1
4 7 1
4 8 1

[ angles ]
; ai aj ak funct c0 c1 c2
c3
2 1 3 1 107.800 276.144 107.800 276.144
2 1 4 1 107.800 276.144 107.800 276.144
2 1 5 1 110.700 313.800 109.500 292.880
3 1 4 1 107.800 276.144 107.800 276.144
3 1 5 1 110.700 313.800 109.500 292.880
4 1 5 1 110.700 313.800 109.500 292.880
1 5 6 1 110.700 313.800 108.500 460.240
1 5 7 1 110.700 313.800 110.700 313.800
1 5 8 1 110.700 313.800 110.700 313.800
6 5 7 1 107.800 276.144 107.800 276.144
6 5 8 1 107.800 276.144 107.800 276.144
7 5 8 1 107.800 276.144 107.800 276.144

[ dihedrals ]
; ai aj ak al funct c0 c1 c2
c3 c4 c5
2 1 5 6 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.94140 2.82420 0.00000
-3.76560 0.00000 0.00000
2 1 5 7 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000
2 1 5 8 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000
3 1 5 6 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.94140 2.82420 0.00000
-3.76560 0.00000 0.00000
3 1 5 7 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000
3 1 5 8 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000
4 1 5 6 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.94140 2.82420 0.00000
-3.76560 0.00000 0.00000
4 1 5 7 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000
4 1 5 8 3 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000 0.62760 1.88280 0.00000
-2.51040 0.00000 0.00000

Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Post by Berk Hess
I assume your are using a united atom model.
What are the bonded parameters for the methanol H
in the ethane topology?
Does it have the normal methanol bond and angle potential?
If these parameters are zero, I would indeed expect a spike.
Berk.
Date: Thu, 14 Feb 2008 18:10:59 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: [gmx-users] Strange dgdl-value together with lincs
Hi
While trying to find, whats wrong with some of my free energy
calculation, I stumbled upon a strange thing, I don't understand.
I'm simulating the interconversion of ethane to methanol in vacuum.
1. charge of perturbed atoms -> 0
2. morph of LJ-parameters and bonds
3. charge of perturbed atoms -> B-state values
Each system is morphed with 1 fs timestep and within 50 steps.
The total energies at the "borders" of the simulations (qqoff/vdw,
vdw/qqon, qqon/vdw and vdw/qqoff) perfectly match and show no jump; the
following simulation always got energies and velocities from the one before.
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_lincs_01.png
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_lincs_10.png
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_nolincs_01.png
For B->A
http://www.mpibpc.mpg.de/groups/grubmueller/start/people/mgoette/images/dgdl_nolincs_10.png
I find this spike in the first step extremely strange. Additionally it
only occurs, when using shake or lincs as constraint solver. I guess, in
a longer process it will contribute to the integral in a neglectable
manner, but maybe theres some bug. The latest CVS shows the same behaviour.
Any comments on that?
Regards
--
Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
.
Berk Hess
2008-02-19 10:20:14 UTC
Permalink
----------------------------------------
Date: Fri, 15 Feb 2008 09:27:24 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Berk,
Its OPLSAA and the hydrogens have bond terms. This is the standard
OPLSAA FF for ethane and methanol.
Thanks.
OK, I found the problem.

You have a delta_lambda of 0.02, which is so high that the kinetic energy
of the constrained degrees of freedom is no longer negligible.
I am not sure if constraints with kinetic energy are formally correct.

But anyhow, you had unconstrained_start=yes, which is probably incorrect,
since I assume your equilibration run had delta_lambda=0 and therefore
in the starting configuration the constraint velocity is zero.
The spike in dG/dl therefore corresponds to building up this kinetic energy.

Now even if you would set unconstrained_start=no, this would not help,
since Gromacs would not put velocities into the constraints at startup.
I have managed to fix this in the code (not committed).
But when delta_lamda is so large, the is another issue.
The coordinates are constrained after the update with lambda of the current
step, but the constraints in the updated configuration should actually obey
the lengths for the next step.
I have to think if things can still be somewhat correct with such a large delta_lambda,
before I starting correcting these two issues.

Berk.

_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
Berk Hess
2008-02-19 12:52:27 UTC
Permalink
----------------------------------------
From: gmx3 at hotmail.com
To: gmx-users at gromacs.org
Subject: RE: [gmx-users] Strange dgdl-value together with lincs
Date: Tue, 19 Feb 2008 11:20:14 +0100
----------------------------------------
Date: Fri, 15 Feb 2008 09:27:24 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Berk,
Its OPLSAA and the hydrogens have bond terms. This is the standard
OPLSAA FF for ethane and methanol.
Thanks.
OK, I found the problem.
You have a delta_lambda of 0.02, which is so high that the kinetic energy
of the constrained degrees of freedom is no longer negligible.
I am not sure if constraints with kinetic energy are formally correct.
But anyhow, you had unconstrained_start=yes, which is probably incorrect,
since I assume your equilibration run had delta_lambda=0 and therefore
in the starting configuration the constraint velocity is zero.
The spike in dG/dl therefore corresponds to building up this kinetic energy.
Now even if you would set unconstrained_start=no, this would not help,
since Gromacs would not put velocities into the constraints at startup.
I have managed to fix this in the code (not committed).
But when delta_lamda is so large, the is another issue.
The coordinates are constrained after the update with lambda of the current
step, but the constraints in the updated configuration should actually obey
the lengths for the next step.
I have to think if things can still be somewhat correct with such a large delta_lambda,
before I starting correcting these two issues.
Berk.
I have now properly fixed this issue in the CVS head branch.
The constraint lengths are now exactly correct at every step
and with unconstrained_start=no the spike in dG/dl due to the velocity is gone.
Note that the error would normally be extreme small, since delta_lambda
is normally extremely small.

Berk.

_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
Maik Goette
2008-02-19 15:07:50 UTC
Permalink
Thanks, Berk

I'm also quite sure, that this spike doesn't explain the deviations I
see in adding up single DG contributions. Maybe I just bring it up once
more:

Assume, you take the interconversion of ethane to methanol in solvent in
one step.
You sample hardcore at 11 evenly spaced lambda values for, lets say, 2
ns each.
You get the dG/dl mean from every run and integrate via simpson to get a
total DG.

Now, you do the same, but split your process into 3. Switching charges
to zero, morph the vdw and bondeds, switching charges back on. All
hardcore to keep comparability.
You take 3x11 lambda snapshots and sample for 666ps each.
You do the integration 3 times like above.
The sum of these 3 DGs should yield the same DG as the run above, I suppose.

Unluckily, both numbers deviate by roughly 10 kJ. I have the feeling,
that hydrogen bonding could play a role, but intuitively, I'd say it
shouldn't differ in the free energy of both.

Thanks again for the help

Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Post by Berk Hess
----------------------------------------
From: gmx3 at hotmail.com
To: gmx-users at gromacs.org
Subject: RE: [gmx-users] Strange dgdl-value together with lincs
Date: Tue, 19 Feb 2008 11:20:14 +0100
----------------------------------------
Date: Fri, 15 Feb 2008 09:27:24 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Berk,
Its OPLSAA and the hydrogens have bond terms. This is the standard
OPLSAA FF for ethane and methanol.
Thanks.
OK, I found the problem.
You have a delta_lambda of 0.02, which is so high that the kinetic energy
of the constrained degrees of freedom is no longer negligible.
I am not sure if constraints with kinetic energy are formally correct.
But anyhow, you had unconstrained_start=yes, which is probably incorrect,
since I assume your equilibration run had delta_lambda=0 and therefore
in the starting configuration the constraint velocity is zero.
The spike in dG/dl therefore corresponds to building up this kinetic energy.
Now even if you would set unconstrained_start=no, this would not help,
since Gromacs would not put velocities into the constraints at startup.
I have managed to fix this in the code (not committed).
But when delta_lamda is so large, the is another issue.
The coordinates are constrained after the update with lambda of the current
step, but the constraints in the updated configuration should actually obey
the lengths for the next step.
I have to think if things can still be somewhat correct with such a large delta_lambda,
before I starting correcting these two issues.
Berk.
I have now properly fixed this issue in the CVS head branch.
The constraint lengths are now exactly correct at every step
and with unconstrained_start=no the spike in dG/dl due to the velocity is gone.
Note that the error would normally be extreme small, since delta_lambda
is normally extremely small.
Berk.
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
.
Berk Hess
2008-02-19 15:21:16 UTC
Permalink
There could be many problems at different points.
One can't say much without having looked in detail at all the data.
Methanol will strongly hydrogen bond to water. So you need to be very careful
that you sample the process of (dis)appearing this hydrogen bond carefully,
both in the 1 and 3 step process. But the "path" could be quite different
in the 1 and 3 step process. So you have to choose your alpha, sigma
and lambda points carefully.

Maybe you could post links to dG/dlambda plots for the 1+3 parts?

Berk.
Date: Tue, 19 Feb 2008 16:07:50 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Thanks, Berk
I'm also quite sure, that this spike doesn't explain the deviations I
see in adding up single DG contributions. Maybe I just bring it up once
Assume, you take the interconversion of ethane to methanol in solvent in
one step.
You sample hardcore at 11 evenly spaced lambda values for, lets say, 2
ns each.
You get the dG/dl mean from every run and integrate via simpson to get a
total DG.
Now, you do the same, but split your process into 3. Switching charges
to zero, morph the vdw and bondeds, switching charges back on. All
hardcore to keep comparability.
You take 3x11 lambda snapshots and sample for 666ps each.
You do the integration 3 times like above.
The sum of these 3 DGs should yield the same DG as the run above, I suppose.
Unluckily, both numbers deviate by roughly 10 kJ. I have the feeling,
that hydrogen bonding could play a role, but intuitively, I'd say it
shouldn't differ in the free energy of both.
Thanks again for the help
Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Post by Berk Hess
----------------------------------------
From: gmx3 at hotmail.com
To: gmx-users at gromacs.org
Subject: RE: [gmx-users] Strange dgdl-value together with lincs
Date: Tue, 19 Feb 2008 11:20:14 +0100
----------------------------------------
Date: Fri, 15 Feb 2008 09:27:24 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Berk,
Its OPLSAA and the hydrogens have bond terms. This is the standard
OPLSAA FF for ethane and methanol.
Thanks.
OK, I found the problem.
You have a delta_lambda of 0.02, which is so high that the kinetic energy
of the constrained degrees of freedom is no longer negligible.
I am not sure if constraints with kinetic energy are formally correct.
But anyhow, you had unconstrained_start=yes, which is probably incorrect,
since I assume your equilibration run had delta_lambda=0 and therefore
in the starting configuration the constraint velocity is zero.
The spike in dG/dl therefore corresponds to building up this kinetic energy.
Now even if you would set unconstrained_start=no, this would not help,
since Gromacs would not put velocities into the constraints at startup.
I have managed to fix this in the code (not committed).
But when delta_lamda is so large, the is another issue.
The coordinates are constrained after the update with lambda of the current
step, but the constraints in the updated configuration should actually obey
the lengths for the next step.
I have to think if things can still be somewhat correct with such a large delta_lambda,
before I starting correcting these two issues.
Berk.
I have now properly fixed this issue in the CVS head branch.
The constraint lengths are now exactly correct at every step
and with unconstrained_start=no the spike in dG/dl due to the velocity is gone.
Note that the error would normally be extreme small, since delta_lambda
is normally extremely small.
Berk.
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
.
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
David Mobley
2008-02-20 16:38:53 UTC
Permalink
Maik,
Post by Maik Goette
Assume, you take the interconversion of ethane to methanol in solvent in
one step.
You sample hardcore at 11 evenly spaced lambda values for, lets say, 2
ns each.
You get the dG/dl mean from every run and integrate via simpson to get a
total DG.
Hold on, if you are doing this with hardcore, there is a singularity
in dG/dlambda. It won't be numerically integrable so whatever you
compute will be in error. See my recent JCP paper that I referred you
to before.

I think if you are using hardcore, you *have* to use slow growth and
the Jarzynski or Crooks expressions, or do some sort of polynomial fit
to the singularity in dG/dlambda and integrate the polynomial.
Integrating dG/dlambda directly will simply fail since you can't
numerically integrate a singularity.

If you're using softcore, you have more options, and simple numerical
integration can work.

For what it's worth, too, a delta_lambda of 0.02 is insanely fast
"slow growth".

David
Maik Goette
2008-02-21 10:34:58 UTC
Permalink
David, Berk,
Post by David Mobley
Hold on, if you are doing this with hardcore, there is a singularity
in dG/dlambda. It won't be numerically integrable so whatever you
compute will be in error. See my recent JCP paper that I referred you
to before.
thats what I thought.
With the latest patched code, I get really good results for the
three-step process and, for the hardcore LJ a deviation, which really
comes from the singularities, I guess.
Now, the 10ns hardcore slow growth and the 3-step discrete lambda
sampling yield perfectly the same number and this is fine.
I will now continue with computing the non-equilibrium values for the
switching process, but I think, those numbers will be ok, too.
Post by David Mobley
For what it's worth, too, a delta_lambda of 0.02 is insanely fast
"slow growth".
agreed ;)
This was just for bug tracking issues. The step I will use for
production will be slightly smaller :)

Thanks for all your help.

Btw, cause I already uploaded the pics, here they are:
The 1-step hardcore process:
Loading Image...
3-step QQ off:
Loading Image...
3-step QQ on:
Loading Image...
3-step LJ and bonded softcore:
Loading Image...
3-step LJ and bonded hardcore:
Loading Image...

This actually proofs again the problem of singularities while using
hardcore LJ morphs. Anyway, the softcore part looks fine and the sum of
the single contributions fits quite well.

As I found out recently, the GROMACS version on the cluster I used, was
without the november bug fix, so I can't tell which fix helped with the
problem. Sorry for that.

Regards

Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Post by David Mobley
Maik,
Post by Maik Goette
Assume, you take the interconversion of ethane to methanol in solvent in
one step.
You sample hardcore at 11 evenly spaced lambda values for, lets say, 2
ns each.
You get the dG/dl mean from every run and integrate via simpson to get a
total DG.
Hold on, if you are doing this with hardcore, there is a singularity
in dG/dlambda. It won't be numerically integrable so whatever you
compute will be in error. See my recent JCP paper that I referred you
to before.
I think if you are using hardcore, you *have* to use slow growth and
the Jarzynski or Crooks expressions, or do some sort of polynomial fit
to the singularity in dG/dlambda and integrate the polynomial.
Integrating dG/dlambda directly will simply fail since you can't
numerically integrate a singularity.
If you're using softcore, you have more options, and simple numerical
integration can work.
For what it's worth, too, a delta_lambda of 0.02 is insanely fast
"slow growth".
David
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
.
van Bemmelen
2008-02-20 07:14:09 UTC
Permalink
Maik,

This is just a hunch, but since you're turning charges on (instead of
off) in the third step, maybe this is also a case of bug 175
(http://bugzilla.gromacs.org/show_bug.cgi?id=175). That is, if you're
using PME for electrostatics.

This bug has been fixed for the CVS version. Another "fix" would be to
do all steps in the reverse direction.

Greetz,
Jeroen
Post by Maik Goette
Unluckily, both numbers deviate by roughly 10 kJ. I have the
feeling, that hydrogen bonding could play a role, but
intuitively, I'd say it shouldn't differ in the free energy of both.
Thanks again for the help
Maik Goette, Dipl. Biol.
Berk Hess
2008-02-20 08:36:47 UTC
Permalink
If I recall correctly, Maik also tried the CVS version.

This bug only has an effect on the PME mesh part.
The effect of this bug is much smaller than 10 kJ/mol.

Berk.
Subject: RE: [gmx-users] Strange dgdl-value together with lincs
Date: Wed, 20 Feb 2008 08:14:09 +0100
From: J.J.M.vanBemmelen at student.TUDelft.NL
To: gmx-users at gromacs.org
Maik,
This is just a hunch, but since you're turning charges on (instead of
off) in the third step, maybe this is also a case of bug 175
(http://bugzilla.gromacs.org/show_bug.cgi?id=175). That is, if you're
using PME for electrostatics.
This bug has been fixed for the CVS version. Another "fix" would be to
do all steps in the reverse direction.
Greetz,
Jeroen
Post by Maik Goette
Unluckily, both numbers deviate by roughly 10 kJ. I have the
feeling, that hydrogen bonding could play a role, but
intuitively, I'd say it shouldn't differ in the free energy of both.
Thanks again for the help
Maik Goette, Dipl. Biol.
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
Maik Goette
2008-02-20 08:38:55 UTC
Permalink
Jeroen,

I'm aware of that bug and using the 3.3.2_pre, with the bugfix included.
Btw, I don't see how to sample "backwards" ;)
Probably you mean in the case of slow growth, no?
In that case, the problem occurs in both directions.

Berk, there are quite a few files (44 in that case).
I will sample the system a bit longer (66ns for each case in total) and
make a gz-download.
Anyway, to me they look normal.

Regards

Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Post by van Bemmelen
Maik,
This is just a hunch, but since you're turning charges on (instead of
off) in the third step, maybe this is also a case of bug 175
(http://bugzilla.gromacs.org/show_bug.cgi?id=175). That is, if you're
using PME for electrostatics.
This bug has been fixed for the CVS version. Another "fix" would be to
do all steps in the reverse direction.
Greetz,
Jeroen
Post by Maik Goette
Unluckily, both numbers deviate by roughly 10 kJ. I have the
feeling, that hydrogen bonding could play a role, but
intuitively, I'd say it shouldn't differ in the free energy of both.
Thanks again for the help
Maik Goette, Dipl. Biol.
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
.
Berk Hess
2008-02-20 08:45:58 UTC
Permalink
No, I meant just the 1+3 plots of as a function of lambda.

Berk.
Date: Wed, 20 Feb 2008 09:38:55 +0100
From: mgoette at mpi-bpc.mpg.de
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Jeroen,
I'm aware of that bug and using the 3.3.2_pre, with the bugfix included.
Btw, I don't see how to sample "backwards" ;)
Probably you mean in the case of slow growth, no?
In that case, the problem occurs in both directions.
Berk, there are quite a few files (44 in that case).
I will sample the system a bit longer (66ns for each case in total) and
make a gz-download.
Anyway, to me they look normal.
Regards
Maik Goette, Dipl. Biol.
Max Planck Institute for Biophysical Chemistry
Theoretical & computational biophysics department
Am Fassberg 11
37077 Goettingen
Germany
Tel. : ++49 551 201 2310
Fax : ++49 551 201 2302
Email : mgoette[at]mpi-bpc.mpg.de
mgoette2[at]gwdg.de
WWW : http://www.mpibpc.gwdg.de/groups/grubmueller/
Post by van Bemmelen
Maik,
This is just a hunch, but since you're turning charges on (instead of
off) in the third step, maybe this is also a case of bug 175
(http://bugzilla.gromacs.org/show_bug.cgi?id=175). That is, if you're
using PME for electrostatics.
This bug has been fixed for the CVS version. Another "fix" would be to
do all steps in the reverse direction.
Greetz,
Jeroen
Post by Maik Goette
Unluckily, both numbers deviate by roughly 10 kJ. I have the
feeling, that hydrogen bonding could play a role, but
intuitively, I'd say it shouldn't differ in the free energy of both.
Thanks again for the help
Maik Goette, Dipl. Biol.
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
.
_______________________________________________
gmx-users mailing list gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
van Bemmelen
2008-02-20 12:07:14 UTC
Permalink
Hi Maik,
From: Maik Goette <mgoette at mpi-bpc.mpg.de>
Subject: Re: [gmx-users] Strange dgdl-value together with lincs
Jeroen,
I'm aware of that bug and using the 3.3.2_pre, with the bugfix
included.
OK. No problem there then. As I said, it was just a guess.
Btw, I don't see how to sample "backwards" ;) Probably you
mean in the case of slow growth, no?
In that case, the problem occurs in both directions.
Oh, LOL, sorry if I was unclear. With "reverse direction" I just meant
making sure charges are always turned off (instead of on) in your
topology by going from the A- to the B-state, in whichever alchemical
free energy scheme you intend to use.

But all of this doesn't matter anyway, since you were already aware of
the bug. :)

Regards,
Jeroen
Loading...