Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Justin Lemkul]

On 5/12/17 2:52 PM, Dan Gil wrote:

I ran equilibrium simulations, and then a longer simulations with the
recommended sd integrator. Still, the mass_lambda dH/dL is nonzero and
sometimes it is the dominant contribution at certain lambdas. I don't think
it is a convergence issue because there aren't huge fluctuations in dH/dL
values seen in dhdl.xvg.

I also ran simulations where only the mass changed as a function of lambda.
Others (i.e. Electrostatic, Van der Waals, bonded, and restraints) were
kept constant. The mass still has an effect on the Helmholtz free energy
change. Is there possibly a mistake with how I've prepared the topology?



Can you provide the section of your .mdp with free energy options?  I don't see 
it in this chain of emails.  It has been suggested previously to leave 
mass-lambdas in one of the end states (A or B, doesn't matter) while doing a 
relative free energy calculation on the assumption that the contributions 
cancel.  Perhaps this is not the case, but I have very little else to contribute 
on the topic because I am not intimately familiar with those aspects of the 
code.  Maybe it's a constraint issue?  You're going from H to F and presumably 
the C-H bonds are to be constrained with the C-F should not be.  I doubt this is 
it (I have successfully done H->Cl transformations and the results match between 
GROMACS and CHARMM), but it's something to check.


-Justin


 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120 12.011
 opls_9620.24012.011
...
21   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
22   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
23   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998

... [ bonds ] ... [ angles ] ... etc...

On Thu, May 11, 2017 at 11:33 AM, Justin Lemkul  wrote:




On 5/11/17 9:49 AM, Dan Gil wrote:


Thank you Dr. Lemkul,

The simulations ran, but I have some a question about the results.

Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
the potential energy does not depend on mass, and my thinking was that the
mass contribution should be very small, if not zero.



It should.  Maybe it's a convergence issue?


-Justin

Best Regards,


Dan

On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul  wrote:




On 5/3/17 10:43 AM, Dan Gil wrote:

Hi Dr. Kausar,


If I want the molecule to have all of its nonbonded interactions with
the
rest of the system (solvent) at each lambda, don't I want to have vdwq
option ON for both couple-lambda0 and couple-lambda1?

What I am attempting to do is mutate molecule A to molecule B through
the
simulation.


Relative free energy calculations are a little bit of a strange beast;

you
do need to set the B-state (couple-lambda1) to "none" and the B-state
types
and charges are read from [atoms].


-Justin

Best Regards,



Dan

On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar <
tasneemkausa...@gmail.com



wrote:


In the last line of mdp file couple-lambda0 and couple-lambda1 have the


same vdwq variable. It indicates A and B state of system. So grompp is
complaining about the identical states.

You can also see the alchemistry.org for detail information about the
free
energy calculations.

On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:

Thank you Dr. Lemkul,



Following your advice, I was able to calculate the solvation free
energy

of


two molecules that I am interested in. The results, unfortunately, were


not


what I expected. I want to try an alternative method, that is, the

alchemical thermodynamic integration method.

If I understand this correctly, I should be able to calculate the free
energy difference between two systems A and B as I change the topology
of
the molecule from A to B in small steps.

I wasn't able to find a tutorial online, but I attempted to try the

method


and I obtained this error from grompp:


WARNING 1 [file grompp.mdp, line 43]:
  The lambda=0 and lambda=1 states for coupling are identical

I was hoping if you could help me understand what I am missing in my
topology or mdp file.

Topology:
 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180
 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120
 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120
 12.011
 opls_9620.24012.011
 4   opls_136  1   HEPTCH2  4 -0.120
 12.011
 opls_9620.24012.011
 5   opls_136  1   HEPTCH2  5 -0.120
 12.011
 

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Dan Gil]

I ran equilibrium simulations, and then a longer simulations with the
recommended sd integrator. Still, the mass_lambda dH/dL is nonzero and
sometimes it is the dominant contribution at certain lambdas. I don't think
it is a convergence issue because there aren't huge fluctuations in dH/dL
values seen in dhdl.xvg.

I also ran simulations where only the mass changed as a function of lambda.
Others (i.e. Electrostatic, Van der Waals, bonded, and restraints) were
kept constant. The mass still has an effect on the Helmholtz free energy
change. Is there possibly a mistake with how I've prepared the topology?

 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120 12.011
 opls_9620.24012.011
...
21   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
22   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
23   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998

... [ bonds ] ... [ angles ] ... etc...

On Thu, May 11, 2017 at 11:33 AM, Justin Lemkul  wrote:

>
>
> On 5/11/17 9:49 AM, Dan Gil wrote:
>
>> Thank you Dr. Lemkul,
>>
>> The simulations ran, but I have some a question about the results.
>>
>> Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
>> the potential energy does not depend on mass, and my thinking was that the
>> mass contribution should be very small, if not zero.
>>
>>
> It should.  Maybe it's a convergence issue?
>
>
> -Justin
>
> Best Regards,
>>
>> Dan
>>
>> On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul  wrote:
>>
>>
>>>
>>> On 5/3/17 10:43 AM, Dan Gil wrote:
>>>
>>> Hi Dr. Kausar,

 If I want the molecule to have all of its nonbonded interactions with
 the
 rest of the system (solvent) at each lambda, don't I want to have vdwq
 option ON for both couple-lambda0 and couple-lambda1?

 What I am attempting to do is mutate molecule A to molecule B through
 the
 simulation.


 Relative free energy calculations are a little bit of a strange beast;
>>> you
>>> do need to set the B-state (couple-lambda1) to "none" and the B-state
>>> types
>>> and charges are read from [atoms].
>>>
>>>
>>> -Justin
>>>
>>> Best Regards,
>>>

 Dan

 On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar <
 tasneemkausa...@gmail.com

>
> wrote:

 In the last line of mdp file couple-lambda0 and couple-lambda1 have the

> same vdwq variable. It indicates A and B state of system. So grompp is
> complaining about the identical states.
>
> You can also see the alchemistry.org for detail information about the
> free
> energy calculations.
>
> On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:
>
> Thank you Dr. Lemkul,
>
>>
>> Following your advice, I was able to calculate the solvation free
>> energy
>>
>> of
>
> two molecules that I am interested in. The results, unfortunately, were
>>
>> not
>
> what I expected. I want to try an alternative method, that is, the
>> alchemical thermodynamic integration method.
>>
>> If I understand this correctly, I should be able to calculate the free
>> energy difference between two systems A and B as I change the topology
>> of
>> the molecule from A to B in small steps.
>>
>> I wasn't able to find a tutorial online, but I attempted to try the
>>
>> method
>
> and I obtained this error from grompp:
>>
>> WARNING 1 [file grompp.mdp, line 43]:
>>   The lambda=0 and lambda=1 states for coupling are identical
>>
>> I was hoping if you could help me understand what I am missing in my
>> topology or mdp file.
>>
>> Topology:
>>  [ moleculetype ]
>> HEPT 3
>>
>>  [ atoms ]
>>  1   opls_135  1   HEPTCH3  1 -0.180
>>  12.011
>>  opls_9610.36012.011
>>  2   opls_136  1   HEPTCH2  2 -0.120
>>  12.011
>>  opls_9620.24012.011
>>  3   opls_136  1   HEPTCH2  3 -0.120
>>  12.011
>>  opls_9620.24012.011
>>  4   opls_136  1   HEPTCH2  4 -0.120
>>  12.011
>>  opls_9620.24012.011
>>  5   opls_136  1   HEPTCH2  5 -0.120
>>  12.011
>>  opls_9620.24012.011
>>  6   opls_136  1   HEPTCH2  6 -0.120
>>  12.011
>>  opls_9620.24012.011
>>  7   opls_135  1   HEPTCH3  7 -0.180
>>  12.011

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Justin Lemkul]

On 5/11/17 9:49 AM, Dan Gil wrote:

Thank you Dr. Lemkul,

The simulations ran, but I have some a question about the results.

Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
the potential energy does not depend on mass, and my thinking was that the
mass contribution should be very small, if not zero.



It should.  Maybe it's a convergence issue?

-Justin


Best Regards,

Dan

On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul  wrote:




On 5/3/17 10:43 AM, Dan Gil wrote:


Hi Dr. Kausar,

If I want the molecule to have all of its nonbonded interactions with the
rest of the system (solvent) at each lambda, don't I want to have vdwq
option ON for both couple-lambda0 and couple-lambda1?

What I am attempting to do is mutate molecule A to molecule B through the
simulation.



Relative free energy calculations are a little bit of a strange beast; you
do need to set the B-state (couple-lambda1) to "none" and the B-state types
and charges are read from [atoms].


-Justin

Best Regards,


Dan

On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar  wrote:

Thank you Dr. Lemkul,


Following your advice, I was able to calculate the solvation free energy


of


two molecules that I am interested in. The results, unfortunately, were


not


what I expected. I want to try an alternative method, that is, the
alchemical thermodynamic integration method.

If I understand this correctly, I should be able to calculate the free
energy difference between two systems A and B as I change the topology
of
the molecule from A to B in small steps.

I wasn't able to find a tutorial online, but I attempted to try the


method


and I obtained this error from grompp:

WARNING 1 [file grompp.mdp, line 43]:
  The lambda=0 and lambda=1 states for coupling are identical

I was hoping if you could help me understand what I am missing in my
topology or mdp file.

Topology:
 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120 12.011
 opls_9620.24012.011
 4   opls_136  1   HEPTCH2  4 -0.120 12.011
 opls_9620.24012.011
 5   opls_136  1   HEPTCH2  5 -0.120 12.011
 opls_9620.24012.011
 6   opls_136  1   HEPTCH2  6 -0.120 12.011
 opls_9620.24012.011
 7   opls_135  1   HEPTCH3  7 -0.180 12.011
 opls_9610.36012.011
 8   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
 9   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
10   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
11   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
12   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
13   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
14   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
15   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
16   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
17   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
18   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
19   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
20   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
21   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
22   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
23   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
... And so on with the bonds, pairs, angles, and dihedrals directive.

MDP File:

;Integration Method and Parameters
integrator   = md
nsteps   = 1
dt = 0.002
nstenergy= 1000
nstlog   = 5000

;Output Control
nstxout = 0
nstvout = 0

;Cutoff Schemes
cutoff-scheme= 

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Dan Gil]

Thank you Dr. Lemkul,

The simulations ran, but I have some a question about the results.

Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
the potential energy does not depend on mass, and my thinking was that the
mass contribution should be very small, if not zero.

Best Regards,

Dan

On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul  wrote:

>
>
> On 5/3/17 10:43 AM, Dan Gil wrote:
>
>> Hi Dr. Kausar,
>>
>> If I want the molecule to have all of its nonbonded interactions with the
>> rest of the system (solvent) at each lambda, don't I want to have vdwq
>> option ON for both couple-lambda0 and couple-lambda1?
>>
>> What I am attempting to do is mutate molecule A to molecule B through the
>> simulation.
>>
>>
> Relative free energy calculations are a little bit of a strange beast; you
> do need to set the B-state (couple-lambda1) to "none" and the B-state types
> and charges are read from [atoms].
>
>
> -Justin
>
> Best Regards,
>>
>> Dan
>>
>> On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar > >
>> wrote:
>>
>> In the last line of mdp file couple-lambda0 and couple-lambda1 have the
>>> same vdwq variable. It indicates A and B state of system. So grompp is
>>> complaining about the identical states.
>>>
>>> You can also see the alchemistry.org for detail information about the
>>> free
>>> energy calculations.
>>>
>>> On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:
>>>
>>> Thank you Dr. Lemkul,

 Following your advice, I was able to calculate the solvation free energy

>>> of
>>>
 two molecules that I am interested in. The results, unfortunately, were

>>> not
>>>
 what I expected. I want to try an alternative method, that is, the
 alchemical thermodynamic integration method.

 If I understand this correctly, I should be able to calculate the free
 energy difference between two systems A and B as I change the topology
 of
 the molecule from A to B in small steps.

 I wasn't able to find a tutorial online, but I attempted to try the

>>> method
>>>
 and I obtained this error from grompp:

 WARNING 1 [file grompp.mdp, line 43]:
   The lambda=0 and lambda=1 states for coupling are identical

 I was hoping if you could help me understand what I am missing in my
 topology or mdp file.

 Topology:
  [ moleculetype ]
 HEPT 3

  [ atoms ]
  1   opls_135  1   HEPTCH3  1 -0.180 12.011
  opls_9610.36012.011
  2   opls_136  1   HEPTCH2  2 -0.120 12.011
  opls_9620.24012.011
  3   opls_136  1   HEPTCH2  3 -0.120 12.011
  opls_9620.24012.011
  4   opls_136  1   HEPTCH2  4 -0.120 12.011
  opls_9620.24012.011
  5   opls_136  1   HEPTCH2  5 -0.120 12.011
  opls_9620.24012.011
  6   opls_136  1   HEPTCH2  6 -0.120 12.011
  opls_9620.24012.011
  7   opls_135  1   HEPTCH3  7 -0.180 12.011
  opls_9610.36012.011
  8   opls_140  1   HEPT  H  1  0.060  1.008
  opls_965   -0.12018.998
  9   opls_140  1   HEPT  H  1  0.060  1.008
  opls_965   -0.12018.998
 10   opls_140  1   HEPT  H  1  0.060  1.008
  opls_965   -0.12018.998
 11   opls_140  1   HEPT  H  2  0.060  1.008
  opls_965   -0.12018.998
 12   opls_140  1   HEPT  H  2  0.060  1.008
  opls_965   -0.12018.998
 13   opls_140  1   HEPT  H  3  0.060  1.008
  opls_965   -0.12018.998
 14   opls_140  1   HEPT  H  3  0.060  1.008
  opls_965   -0.12018.998
 15   opls_140  1   HEPT  H  4  0.060  1.008
  opls_965   -0.12018.998
 16   opls_140  1   HEPT  H  4  0.060  1.008
  opls_965   -0.12018.998
 17   opls_140  1   HEPT  H  5  0.060  1.008
  opls_965   -0.12018.998
 18   opls_140  1   HEPT  H  5  0.060  1.008
  opls_965   -0.12018.998
 19   opls_140  1   HEPT  H  6  0.060  1.008
  opls_965   -0.12018.998
 20   opls_140  1   HEPT  H  6  0.060  1.008
  opls_965   -0.12018.998
 21   opls_140  1   HEPT  H  7  0.060  1.008
  opls_965   -0.12018.998
 22   opls_140  1   HEPT  H  7  0.060  1.008
  opls_965   -0.12018.998
 23   opls_140  1   HEPT  H  7  0.060  1.008
 

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Justin Lemkul]

On 5/3/17 10:43 AM, Dan Gil wrote:

Hi Dr. Kausar,

If I want the molecule to have all of its nonbonded interactions with the
rest of the system (solvent) at each lambda, don't I want to have vdwq
option ON for both couple-lambda0 and couple-lambda1?

What I am attempting to do is mutate molecule A to molecule B through the
simulation.



Relative free energy calculations are a little bit of a strange beast; you do 
need to set the B-state (couple-lambda1) to "none" and the B-state types and 
charges are read from [atoms].


-Justin


Best Regards,

Dan

On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar 
wrote:


In the last line of mdp file couple-lambda0 and couple-lambda1 have the
same vdwq variable. It indicates A and B state of system. So grompp is
complaining about the identical states.

You can also see the alchemistry.org for detail information about the free
energy calculations.

On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:


Thank you Dr. Lemkul,

Following your advice, I was able to calculate the solvation free energy

of

two molecules that I am interested in. The results, unfortunately, were

not

what I expected. I want to try an alternative method, that is, the
alchemical thermodynamic integration method.

If I understand this correctly, I should be able to calculate the free
energy difference between two systems A and B as I change the topology of
the molecule from A to B in small steps.

I wasn't able to find a tutorial online, but I attempted to try the

method

and I obtained this error from grompp:

WARNING 1 [file grompp.mdp, line 43]:
  The lambda=0 and lambda=1 states for coupling are identical

I was hoping if you could help me understand what I am missing in my
topology or mdp file.

Topology:
 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120 12.011
 opls_9620.24012.011
 4   opls_136  1   HEPTCH2  4 -0.120 12.011
 opls_9620.24012.011
 5   opls_136  1   HEPTCH2  5 -0.120 12.011
 opls_9620.24012.011
 6   opls_136  1   HEPTCH2  6 -0.120 12.011
 opls_9620.24012.011
 7   opls_135  1   HEPTCH3  7 -0.180 12.011
 opls_9610.36012.011
 8   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
 9   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
10   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
11   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
12   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
13   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
14   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
15   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
16   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
17   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
18   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
19   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
20   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
21   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
22   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
23   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
... And so on with the bonds, pairs, angles, and dihedrals directive.

MDP File:

;Integration Method and Parameters
integrator   = md
nsteps   = 1
dt = 0.002
nstenergy= 1000
nstlog   = 5000

;Output Control
nstxout = 0
nstvout = 0

;Cutoff Schemes
cutoff-scheme= group
rlist= 1.0
vdw-type = cut-off
rvdw = 2.0

;Coulomb interactions
coulombtype  = pme
rcoulomb = 1.0
fourierspacing   = 0.4

;Constraints
constraints  = all-bonds

;Temperature coupling
tcoupl   = v-rescale
tc-grps  = system
tau-t= 0.1
ref-t= 300

;Pressure coupling
pcoupl = parrinello-rahman
ref-p = 1.01325

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Dan Gil]

Hi Dr. Kausar,

If I want the molecule to have all of its nonbonded interactions with the
rest of the system (solvent) at each lambda, don't I want to have vdwq
option ON for both couple-lambda0 and couple-lambda1?

What I am attempting to do is mutate molecule A to molecule B through the
simulation.

Best Regards,

Dan

On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar 
wrote:

> In the last line of mdp file couple-lambda0 and couple-lambda1 have the
> same vdwq variable. It indicates A and B state of system. So grompp is
> complaining about the identical states.
>
> You can also see the alchemistry.org for detail information about the free
> energy calculations.
>
> On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:
>
> > Thank you Dr. Lemkul,
> >
> > Following your advice, I was able to calculate the solvation free energy
> of
> > two molecules that I am interested in. The results, unfortunately, were
> not
> > what I expected. I want to try an alternative method, that is, the
> > alchemical thermodynamic integration method.
> >
> > If I understand this correctly, I should be able to calculate the free
> > energy difference between two systems A and B as I change the topology of
> > the molecule from A to B in small steps.
> >
> > I wasn't able to find a tutorial online, but I attempted to try the
> method
> > and I obtained this error from grompp:
> >
> > WARNING 1 [file grompp.mdp, line 43]:
> >   The lambda=0 and lambda=1 states for coupling are identical
> >
> > I was hoping if you could help me understand what I am missing in my
> > topology or mdp file.
> >
> > Topology:
> >  [ moleculetype ]
> > HEPT 3
> >
> >  [ atoms ]
> >  1   opls_135  1   HEPTCH3  1 -0.180 12.011
> >  opls_9610.36012.011
> >  2   opls_136  1   HEPTCH2  2 -0.120 12.011
> >  opls_9620.24012.011
> >  3   opls_136  1   HEPTCH2  3 -0.120 12.011
> >  opls_9620.24012.011
> >  4   opls_136  1   HEPTCH2  4 -0.120 12.011
> >  opls_9620.24012.011
> >  5   opls_136  1   HEPTCH2  5 -0.120 12.011
> >  opls_9620.24012.011
> >  6   opls_136  1   HEPTCH2  6 -0.120 12.011
> >  opls_9620.24012.011
> >  7   opls_135  1   HEPTCH3  7 -0.180 12.011
> >  opls_9610.36012.011
> >  8   opls_140  1   HEPT  H  1  0.060  1.008
> >  opls_965   -0.12018.998
> >  9   opls_140  1   HEPT  H  1  0.060  1.008
> >  opls_965   -0.12018.998
> > 10   opls_140  1   HEPT  H  1  0.060  1.008
> >  opls_965   -0.12018.998
> > 11   opls_140  1   HEPT  H  2  0.060  1.008
> >  opls_965   -0.12018.998
> > 12   opls_140  1   HEPT  H  2  0.060  1.008
> >  opls_965   -0.12018.998
> > 13   opls_140  1   HEPT  H  3  0.060  1.008
> >  opls_965   -0.12018.998
> > 14   opls_140  1   HEPT  H  3  0.060  1.008
> >  opls_965   -0.12018.998
> > 15   opls_140  1   HEPT  H  4  0.060  1.008
> >  opls_965   -0.12018.998
> > 16   opls_140  1   HEPT  H  4  0.060  1.008
> >  opls_965   -0.12018.998
> > 17   opls_140  1   HEPT  H  5  0.060  1.008
> >  opls_965   -0.12018.998
> > 18   opls_140  1   HEPT  H  5  0.060  1.008
> >  opls_965   -0.12018.998
> > 19   opls_140  1   HEPT  H  6  0.060  1.008
> >  opls_965   -0.12018.998
> > 20   opls_140  1   HEPT  H  6  0.060  1.008
> >  opls_965   -0.12018.998
> > 21   opls_140  1   HEPT  H  7  0.060  1.008
> >  opls_965   -0.12018.998
> > 22   opls_140  1   HEPT  H  7  0.060  1.008
> >  opls_965   -0.12018.998
> > 23   opls_140  1   HEPT  H  7  0.060  1.008
> >  opls_965   -0.12018.998
> > ... And so on with the bonds, pairs, angles, and dihedrals directive.
> >
> > MDP File:
> >
> > ;Integration Method and Parameters
> > integrator   = md
> > nsteps   = 1
> > dt = 0.002
> > nstenergy= 1000
> > nstlog   = 5000
> >
> > ;Output Control
> > nstxout = 0
> > nstvout = 0
> >
> > ;Cutoff Schemes
> > cutoff-scheme= group
> > rlist= 1.0
> > vdw-type = cut-off
> > rvdw = 2.0
> >
> > ;Coulomb interactions
> > coulombtype  = pme
> > rcoulomb = 1.0
> > fourierspacing   = 0.4
> >
> > ;Constraints
> > constraints  = all-bonds
> >
> > ;Temperature coupling
> > tcoupl   = 

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Tasneem Kausar]

In the last line of mdp file couple-lambda0 and couple-lambda1 have the
same vdwq variable. It indicates A and B state of system. So grompp is
complaining about the identical states.

You can also see the alchemistry.org for detail information about the free
energy calculations.

On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:

> Thank you Dr. Lemkul,
>
> Following your advice, I was able to calculate the solvation free energy of
> two molecules that I am interested in. The results, unfortunately, were not
> what I expected. I want to try an alternative method, that is, the
> alchemical thermodynamic integration method.
>
> If I understand this correctly, I should be able to calculate the free
> energy difference between two systems A and B as I change the topology of
> the molecule from A to B in small steps.
>
> I wasn't able to find a tutorial online, but I attempted to try the method
> and I obtained this error from grompp:
>
> WARNING 1 [file grompp.mdp, line 43]:
>   The lambda=0 and lambda=1 states for coupling are identical
>
> I was hoping if you could help me understand what I am missing in my
> topology or mdp file.
>
> Topology:
>  [ moleculetype ]
> HEPT 3
>
>  [ atoms ]
>  1   opls_135  1   HEPTCH3  1 -0.180 12.011
>  opls_9610.36012.011
>  2   opls_136  1   HEPTCH2  2 -0.120 12.011
>  opls_9620.24012.011
>  3   opls_136  1   HEPTCH2  3 -0.120 12.011
>  opls_9620.24012.011
>  4   opls_136  1   HEPTCH2  4 -0.120 12.011
>  opls_9620.24012.011
>  5   opls_136  1   HEPTCH2  5 -0.120 12.011
>  opls_9620.24012.011
>  6   opls_136  1   HEPTCH2  6 -0.120 12.011
>  opls_9620.24012.011
>  7   opls_135  1   HEPTCH3  7 -0.180 12.011
>  opls_9610.36012.011
>  8   opls_140  1   HEPT  H  1  0.060  1.008
>  opls_965   -0.12018.998
>  9   opls_140  1   HEPT  H  1  0.060  1.008
>  opls_965   -0.12018.998
> 10   opls_140  1   HEPT  H  1  0.060  1.008
>  opls_965   -0.12018.998
> 11   opls_140  1   HEPT  H  2  0.060  1.008
>  opls_965   -0.12018.998
> 12   opls_140  1   HEPT  H  2  0.060  1.008
>  opls_965   -0.12018.998
> 13   opls_140  1   HEPT  H  3  0.060  1.008
>  opls_965   -0.12018.998
> 14   opls_140  1   HEPT  H  3  0.060  1.008
>  opls_965   -0.12018.998
> 15   opls_140  1   HEPT  H  4  0.060  1.008
>  opls_965   -0.12018.998
> 16   opls_140  1   HEPT  H  4  0.060  1.008
>  opls_965   -0.12018.998
> 17   opls_140  1   HEPT  H  5  0.060  1.008
>  opls_965   -0.12018.998
> 18   opls_140  1   HEPT  H  5  0.060  1.008
>  opls_965   -0.12018.998
> 19   opls_140  1   HEPT  H  6  0.060  1.008
>  opls_965   -0.12018.998
> 20   opls_140  1   HEPT  H  6  0.060  1.008
>  opls_965   -0.12018.998
> 21   opls_140  1   HEPT  H  7  0.060  1.008
>  opls_965   -0.12018.998
> 22   opls_140  1   HEPT  H  7  0.060  1.008
>  opls_965   -0.12018.998
> 23   opls_140  1   HEPT  H  7  0.060  1.008
>  opls_965   -0.12018.998
> ... And so on with the bonds, pairs, angles, and dihedrals directive.
>
> MDP File:
>
> ;Integration Method and Parameters
> integrator   = md
> nsteps   = 1
> dt = 0.002
> nstenergy= 1000
> nstlog   = 5000
>
> ;Output Control
> nstxout = 0
> nstvout = 0
>
> ;Cutoff Schemes
> cutoff-scheme= group
> rlist= 1.0
> vdw-type = cut-off
> rvdw = 2.0
>
> ;Coulomb interactions
> coulombtype  = pme
> rcoulomb = 1.0
> fourierspacing   = 0.4
>
> ;Constraints
> constraints  = all-bonds
>
> ;Temperature coupling
> tcoupl   = v-rescale
> tc-grps  = system
> tau-t= 0.1
> ref-t= 300
>
> ;Pressure coupling
> pcoupl = parrinello-rahman
> ref-p = 1.01325
> compressibility = 4.5e-5
> tau-p = 5
>
> ;Free energy calculation
> free-energy  = yes
> init-lambda  = 0
> delta-lambda = 0.1
> couple-moltype   = HEPT
> couple-lambda0   = vdwq
> couple-lambda1   = vdwq
>
> On Mon, Apr 3, 2017 at 2:05 PM, Justin Lemkul  wrote:
>
> >
> >
> > On 4/3/17 2:01 PM, Dan Gil wrote:
> >
> >> Thank you Dr. Lemkul,
> >>
> >> I am trying to run 

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Dan Gil]

Thank you Dr. Lemkul,

Following your advice, I was able to calculate the solvation free energy of
two molecules that I am interested in. The results, unfortunately, were not
what I expected. I want to try an alternative method, that is, the
alchemical thermodynamic integration method.

If I understand this correctly, I should be able to calculate the free
energy difference between two systems A and B as I change the topology of
the molecule from A to B in small steps.

I wasn't able to find a tutorial online, but I attempted to try the method
and I obtained this error from grompp:

WARNING 1 [file grompp.mdp, line 43]:
  The lambda=0 and lambda=1 states for coupling are identical

I was hoping if you could help me understand what I am missing in my
topology or mdp file.

Topology:
 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120 12.011
 opls_9620.24012.011
 4   opls_136  1   HEPTCH2  4 -0.120 12.011
 opls_9620.24012.011
 5   opls_136  1   HEPTCH2  5 -0.120 12.011
 opls_9620.24012.011
 6   opls_136  1   HEPTCH2  6 -0.120 12.011
 opls_9620.24012.011
 7   opls_135  1   HEPTCH3  7 -0.180 12.011
 opls_9610.36012.011
 8   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
 9   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
10   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
11   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
12   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
13   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
14   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
15   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
16   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
17   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
18   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
19   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
20   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
21   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
22   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
23   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
... And so on with the bonds, pairs, angles, and dihedrals directive.

MDP File:

;Integration Method and Parameters
integrator   = md
nsteps   = 1
dt = 0.002
nstenergy= 1000
nstlog   = 5000

;Output Control
nstxout = 0
nstvout = 0

;Cutoff Schemes
cutoff-scheme= group
rlist= 1.0
vdw-type = cut-off
rvdw = 2.0

;Coulomb interactions
coulombtype  = pme
rcoulomb = 1.0
fourierspacing   = 0.4

;Constraints
constraints  = all-bonds

;Temperature coupling
tcoupl   = v-rescale
tc-grps  = system
tau-t= 0.1
ref-t= 300

;Pressure coupling
pcoupl = parrinello-rahman
ref-p = 1.01325
compressibility = 4.5e-5
tau-p = 5

;Free energy calculation
free-energy  = yes
init-lambda  = 0
delta-lambda = 0.1
couple-moltype   = HEPT
couple-lambda0   = vdwq
couple-lambda1   = vdwq

On Mon, Apr 3, 2017 at 2:05 PM, Justin Lemkul  wrote:

>
>
> On 4/3/17 2:01 PM, Dan Gil wrote:
>
>> Thank you Dr. Lemkul,
>>
>> I am trying to run grompp with the md integrator, but I am getting this
>> error:
>> "For proper sampling of the (nearly) decoupled state, stochastic dynamics
>> should be used"
>>
>> Should I ignore this warning with the -maxwarn option and try running it?
>> I
>> will see if I obtain comparable values for ethanol.
>>
>>
> People have already done such a comparison and that's why grompp is
> telling you this - it's better to use the Langevin integrator.  You'll get
> better sampling, particularly towards the end states.  Using -maxwarn tells
> grompp "you're trying to prevent me from making a mistake, but 

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Justin Lemkul]

On 4/3/17 2:01 PM, Dan Gil wrote:

Thank you Dr. Lemkul,

I am trying to run grompp with the md integrator, but I am getting this
error:
"For proper sampling of the (nearly) decoupled state, stochastic dynamics
should be used"

Should I ignore this warning with the -maxwarn option and try running it? I
will see if I obtain comparable values for ethanol.



People have already done such a comparison and that's why grompp is telling you 
this - it's better to use the Langevin integrator.  You'll get better sampling, 
particularly towards the end states.  Using -maxwarn tells grompp "you're trying 
to prevent me from making a mistake, but I want to do it anyway" :)  You'd 
better have a really, really good reason to try to override it.


-Justin


Best Regards,

Dan


On Mon, Mar 27, 2017 at 12:51 PM, Justin Lemkul  wrote:




On 3/26/17 9:40 PM, Dan Gil wrote:


Hi,

I am following Dr. Sander Pronk's and Dr. Justin Lemkul's tutorial on
calculating free energy of solvation. Is it possible and theoretically
sound to use the md integrator instead of the sd integrator for these
calculations?



Langevin dynamics gives better sampling so it is frequently used for free
energy calculations.  You may get comparable results with the leap-frog
integrator, but I haven never done a side-by-side comparison.

-Justin

I have already done a considerable amount of work using md integration, and

I want to make sure that the free energy values I calculate are consistent
with my previous work.

If using the md integrator is not sound, is there an alternative way of
calculating solvation energy that will be consistent?

Best Regards,

Dan



--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support
/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
send a mail to gmx-users-requ...@gromacs.org.



--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Dan Gil]

Thank you Dr. Lemkul,

I am trying to run grompp with the md integrator, but I am getting this
error:
"For proper sampling of the (nearly) decoupled state, stochastic dynamics
should be used"

Should I ignore this warning with the -maxwarn option and try running it? I
will see if I obtain comparable values for ethanol.

Best Regards,

Dan


On Mon, Mar 27, 2017 at 12:51 PM, Justin Lemkul  wrote:

>
>
> On 3/26/17 9:40 PM, Dan Gil wrote:
>
>> Hi,
>>
>> I am following Dr. Sander Pronk's and Dr. Justin Lemkul's tutorial on
>> calculating free energy of solvation. Is it possible and theoretically
>> sound to use the md integrator instead of the sd integrator for these
>> calculations?
>>
>>
> Langevin dynamics gives better sampling so it is frequently used for free
> energy calculations.  You may get comparable results with the leap-frog
> integrator, but I haven never done a side-by-side comparison.
>
> -Justin
>
> I have already done a considerable amount of work using md integration, and
>> I want to make sure that the free energy values I calculate are consistent
>> with my previous work.
>>
>> If using the md integrator is not sound, is there an alternative way of
>> calculating solvation energy that will be consistent?
>>
>> Best Regards,
>>
>> Dan
>>
>>
> --
> ==
>
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 629
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalem...@outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support
> /Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-requ...@gromacs.org.
>
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Using the md integrator for calculating free energy of solvation

[Justin Lemkul]

On 3/26/17 9:40 PM, Dan Gil wrote:

Hi,

I am following Dr. Sander Pronk's and Dr. Justin Lemkul's tutorial on
calculating free energy of solvation. Is it possible and theoretically
sound to use the md integrator instead of the sd integrator for these
calculations?



Langevin dynamics gives better sampling so it is frequently used for free energy 
calculations.  You may get comparable results with the leap-frog integrator, but 
I haven never done a side-by-side comparison.


-Justin


I have already done a considerable amount of work using md integration, and
I want to make sure that the free energy values I calculate are consistent
with my previous work.

If using the md integrator is not sound, is there an alternative way of
calculating solvation energy that will be consistent?

Best Regards,

Dan



--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

[gmx-users] Using the md integrator for calculating free energy of solvation

[Dan Gil]

Hi,

I am following Dr. Sander Pronk's and Dr. Justin Lemkul's tutorial on
calculating free energy of solvation. Is it possible and theoretically
sound to use the md integrator instead of the sd integrator for these
calculations?

I have already done a considerable amount of work using md integration, and
I want to make sure that the free energy values I calculate are consistent
with my previous work.

If using the md integrator is not sound, is there an alternative way of
calculating solvation energy that will be consistent?

Best Regards,

Dan
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.