[gmx-users] The sign of deuterium order parameter calculated by g_order

[Jason Zhu]

Dear All,

I am studying the effects of cholesterols in POPC lipid bilayer. The
Gromacs is the version 4.6.5 in my simulations.

By using the g_order, I could output the deuterium order parameter of POPC
lipids for a long-term equilibration.

"g_order_mpi -s msd.tpr -f msd-mol.xtc -n sn1.ndx -d z -od deuter_sn1.xvg
-o order_sn1.xvg"

The results in the file of "deuter_sn1.xvg" are given as below:

@title "Deuterium order parameters"
@xaxis  label "Atom"
@yaxis  label "Scd"
@TYPE xy
   1   0.204493
   2   0.223424
   3   0.222852
   4   0.228435
   5   0.230187
   6   0.229785
   7   0.223191
   8   0.217353
   9   0.204112
  10   0.193993
  11   0.176565
  12   0.162058
  130.13807
  14   0.112887

As in the experimental papers about the deuterium order parameter measured
by NMR, the Scd is negative and round -0.2.

I understand the sign of the deuterium order parameter is not measurable in
conventional NMR. I am wondering the values calculated and output by
g_order in Gromacs is the absolute value (|Scd|) or the opposite value
(-Scd) by default.

Sorry to bother you with this simple question. But I didn't find any clear
answer from mail list or manual. It would be appreciated if you could give
your insights or comments.

Best,
Wenpeng Zhu
-- 
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] Periodic Molecule's Free Energy Calculation Error

[Jason Zhu]

Dear Alex,

Many thanks for your replies and helping.

I totally agree with your comments about infinite molecules and periodic
boundary condition.

Here in our simulations, the hBN sheet covers the box dimensions (6nm*6nm)
without free edges. By using PBC and "periodic-molecules = yes", we could
calculate the surface energy of hBN sheet without effects of edges.

The parameters of force field we are using are borrowed from the following
papers. We introduced them into the Gromos force field. We are having a new
publication using this modified force field in which all the parameters are
shown explicitly. I could send it to you when it is published online. The
functions and parameters are fitted by the setting of "nrexcl=3" in these
papers. We couldn't make any changes for this. But we could try to use
larger value of "nrexcl" to fit the force field of hBN and include
short-range non-bonded interactions into bonded interactions. Thank you for
your suggestions.

1. Hilder, T. A. et al. Validity of current force fields for simulations on
boron nitride nanotubes. IET Micro & Nano Letters 5, 150-156,
doi:10.1049/mnl.2009.0112 (2010).

2. Kamath, G. & Baker, G. A. Are ionic liquids suitable media for boron
nitride exfoliation and dispersion? Insight via molecular dynamics. RSC
Advances 3, 8197-8202, doi:10.1039/c3ra40488a (2013).

3. Wu, J., Wang, B., Wei, Y., Yang, R. & Dresselhaus, M. Mechanics and
Mechanically Tunable Band Gap in Single-Layer Hexagonal Boron-Nitride.
Materials Research Letters 1, 200-206, doi:10.1080/21663831.2013.824516
(2013).

Best,
Jason


Message: 5
Date: Wed, 12 Jul 2017 19:12:26 -0600
From: Alex 
To: Discussion list for GROMACS users 
Subject: Re: [gmx-users] Periodic Molecule's Free Energy Calculation
Error
Message-ID: <4b927ef2-9a56-6655-8053-d284cf88b...@gmail.com>
Content-Type: text/plain; charset=utf-8; format=flowed


>
> I wonder if the "couple-intramol = yes" is a must. Does it have any
> influence on the output results if we turn off the intra-molecular
> non-bonded interactions of a large infinite molecule?
>
The answer to your question has nothing to do with Gromacs, but with
understanding the difference between crystals and biomolecules (for
which Gromacs was designed).
Also (unrelated), it is a common misconception to believe that PBC makes
something infinite -- the effective size of your system is entirely
determined by the supercell size (proof: consider the ripples in hBN and
determine the lowest wavelength of the ripple that can propagate -- it
is commensurate with the box size). In an infinite system, you can have
an immensely long wave (though not infinite, as shown by Landau a while
back). PBC does not make anything infinite, it is a mathematical way of
avoiding surfaces.
>
> There is no universal force field for HBN, so I am using a modified
> gromos54a7_atb force field, i.e., manually adding the parameters for
> boron and nitrogen to the bonded & nonbonded .itp files.
Oh, I know that there is no force fields for these structures. ;) My
question was about which Gromacs ff you were using to insert your
parameters, and, most importantly, where those parameters came from.

> The parameters are obtained from literature.
>
What literature? All bio-style ff adaptations of solid-state potentials
(e.g. Tersoff-Brenner for hBN) I am aware of make it very clear that
"intramolecular" interactions between atoms sharing up to a fairly
distant covalently bound neighbor are limited to bonds and angles. This
comes from the math involved in developing potentials for crystals.
There was a recent question regarding this very problem here, which was
solved by setting a larger nrexcl value. In your case, you solved it
with turning off intramolecular coupling. In fact, if you set your
nrexcl to something like 4 or 5, you may not even need to turn off the
coupling. But then again, I don't know where the parameters came from.

Alex
-- 
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] (no subject)

[Jason Zhu]

Dear Alex,

Many thanks for your replies and helping.

I totally agree with your comments about infinite molecules and periodic
boundary condition.

Here in our simulations, the hBN sheet covers the box dimensions (6nm*6nm)
without free edges. By using PBC and "periodic-molecules = yes", we could
calculate the surface energy of hBN sheet without effects of edges.

The parameters of force field we are using are borrowed from the following
papers. We introduced them into the Gromos force field. We are having a new
publication using this modified force field in which all the parameters are
shown explicitly. I could send it to you when it is published online. The
functions and parameters are fitted by the setting of "nrexcl=3" in these
papers. We couldn't make any changes for this. But we could try to use
larger value of "nrexcl" to fit the force field of hBN and include
short-range non-bonded interactions into bonded interactions. Thank you for
your suggestions.

1. Hilder, T. A. et al. Validity of current force fields for simulations on
boron nitride nanotubes. IET Micro & Nano Letters 5, 150-156,
doi:10.1049/mnl.2009.0112 (2010).

2. Kamath, G. & Baker, G. A. Are ionic liquids suitable media for boron
nitride exfoliation and dispersion? Insight via molecular dynamics. RSC
Advances 3, 8197-8202, doi:10.1039/c3ra40488a (2013).

3. Wu, J., Wang, B., Wei, Y., Yang, R. & Dresselhaus, M. Mechanics and
Mechanically Tunable Band Gap in Single-Layer Hexagonal Boron-Nitride.
Materials Research Letters 1, 200-206, doi:10.1080/21663831.2013.824516
(2013).

Best,
Jason


Message: 5
Date: Wed, 12 Jul 2017 19:12:26 -0600
From: Alex 
To: Discussion list for GROMACS users 
Subject: Re: [gmx-users] Periodic Molecule's Free Energy Calculation
Error
Message-ID: <4b927ef2-9a56-6655-8053-d284cf88b...@gmail.com>
Content-Type: text/plain; charset=utf-8; format=flowed


>
> I wonder if the "couple-intramol = yes" is a must. Does it have any
> influence on the output results if we turn off the intra-molecular
> non-bonded interactions of a large infinite molecule?
>
The answer to your question has nothing to do with Gromacs, but with
understanding the difference between crystals and biomolecules (for
which Gromacs was designed).
Also (unrelated), it is a common misconception to believe that PBC makes
something infinite -- the effective size of your system is entirely
determined by the supercell size (proof: consider the ripples in hBN and
determine the lowest wavelength of the ripple that can propagate -- it
is commensurate with the box size). In an infinite system, you can have
an immensely long wave (though not infinite, as shown by Landau a while
back). PBC does not make anything infinite, it is a mathematical way of
avoiding surfaces.
>
> There is no universal force field for HBN, so I am using a modified
> gromos54a7_atb force field, i.e., manually adding the parameters for
> boron and nitrogen to the bonded & nonbonded .itp files.
Oh, I know that there is no force fields for these structures. ;) My
question was about which Gromacs ff you were using to insert your
parameters, and, most importantly, where those parameters came from.

> The parameters are obtained from literature.
>
What literature? All bio-style ff adaptations of solid-state potentials
(e.g. Tersoff-Brenner for hBN) I am aware of make it very clear that
"intramolecular" interactions between atoms sharing up to a fairly
distant covalently bound neighbor are limited to bonds and angles. This
comes from the math involved in developing potentials for crystals.
There was a recent question regarding this very problem here, which was
solved by setting a larger nrexcl value. In your case, you solved it
with turning off intramolecular coupling. In fact, if you set your
nrexcl to something like 4 or 5, you may not even need to turn off the
coupling. But then again, I don't know where the parameters came from.

Alex
-- 
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] Periodic Molecule's Free Energy Calculation Error

[Jason Zhu]

Dear Justin,

Many thanks for your reply and explanation.

Now I understand the difference between "couple-intramol = no"
and "couple-intramol = yes" and their consequences.

According to the manual, it is recommended to use "couple-intramol = yes"
for relatively large molecules. Why is that?

I don't know why "couple-intramol = no" doesn't work on large molecules or
infinite molecules with periodic boundary condition.

Is it because the code of "couple-intramol = no" doesn't consider molecules
longer than box size or for any consideration.

Looking forward to hearing from you. Thank you again.

Best,
Jason

On 7/12/17 9:12 PM, Alex wrote:
>
>>
>> I wonder if the "couple-intramol = yes" is a must. Does it have any
influence
>> on the output results if we turn off the intra-molecular non-bonded
>> interactions of a large infinite molecule?
>>
> The answer to your question has nothing to do with Gromacs, but with
> understanding the difference between crystals and biomolecules (for which
> Gromacs was designed).

There is a functional difference between coupling intramolecular
interactions
and not, which will affect the computed free energy.  With couple-intramol
= no,
you get the correct vacuum state of the solute molecule; with
couple-intramol =
yes, you get perturbed intramolecular interactions.  This has important
consequences for what one is trying to compute.

-Justin

> Also (unrelated), it is a common misconception to believe that PBC makes
> something infinite -- the effective size of your system is entirely
determined
> by the supercell size (proof: consider the ripples in hBN and determine
the
> lowest wavelength of the ripple that can propagate -- it is commensurate
with
> the box size). In an infinite system, you can have an immensely long wave
> (though not infinite, as shown by Landau a while back). PBC does not make
> anything infinite, it is a mathematical way of avoiding surfaces.
>>
>> There is no universal force field for HBN, so I am using a modified
>> gromos54a7_atb force field, i.e., manually adding the parameters for
boron and
>> nitrogen to the bonded & nonbonded .itp files.
> Oh, I know that there is no force fields for these structures. ;) My
question
> was about which Gromacs ff you were using to insert your parameters, and,
most
> importantly, where those parameters came from.
>
>> The parameters are obtained from literature.
>>
> What literature? All bio-style ff adaptations of solid-state potentials
(e.g.
> Tersoff-Brenner for hBN) I am aware of make it very clear that
"intramolecular"
> interactions between atoms sharing up to a fairly distant covalently bound
> neighbor are limited to bonds and angles. This comes from the math
involved in
> developing potentials for crystals. There was a recent question regarding
this
> very problem here, which was solved by setting a larger nrexcl value. In
your
> case, you solved it with turning off intramolecular coupling. In fact, if
you
> set your nrexcl to something like 4 or 5, you may not even need to turn
off the
> coupling. But then again, I don't know where the parameters came from.
>
> Alex

--
==

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] Periodic Molecule's Free Energy Calculation Error

[Jason Zhu]

Hi Alex,

Thank you for your reply. Here below is my response to your comments.


hBN is hardly a common subject of simulation for Gromacs folks, but let's
see...

1. When you run the simulation in vacuum, do you get the same error?
Does a 300K vacuum simulation result in reasonable sheet behavior? What
about NVT?


The answer is yes. I ran more tests today, and I think I have found exactly
how the error came out.
In fact, if I remove the whole free energy calculation section from my mdp
file, the system will run perfectly in parallel, both HBN in vacuum and HBN
in solution.
So I suspect somewhere in FEC section is causing the trouble.
Then I looked up each command in GROMACS documentation, and found this
"couple-intramol" flag:

"couple-intramol:
no
All intra-molecular non-bonded interactions for moleculetype couple-moltype
are replaced by exclusions and explicit pair interactions. In this manner
the decoupled state of the molecule corresponds to the proper vacuum state
without periodicity effects.
yes
The intra-molecular Van der Waals and Coulomb interactions are also turned
on/off. This can be useful for partitioning free-energies of relatively
large molecules, where the intra-molecular non-bonded interactions might
lead to kinetically trapped vacuum conformations. The 1-4 pair interactions
are not turned off."

I wonder if the "couple-intramol = no" works periodic molecules.

In my mdp file, I set the flag to "no", because both reference tutorials
set it to "no".
However, I am dealing with a relative large molecule with periodicity, so I
should have set this flag to "yes".

Then I ran several more test with "couple-intramol = yes". Now, for
EM/NVT/NPT/PROD_MD, the system could finally run in parallel. No more
disturbing errors. Hooray!!

But I am still confused about the setting of "couple-intramol". Why it can
run parallelly if "couple-intramol" is set "yes". The hBN sheet in my
simulation is infinite with PBC.

I wonder if the "couple-intramol = yes" is a must. Does it have any
influence on the output results if we turn off the intra-molecular
non-bonded interactions of a large infinite molecule?

Please note that "periodic-molecules = yes" in all my simulations.


2. What GROMACS forcefield are you using and what are the nonbonded
types for boron and nitrogen?


There is no universal force field for HBN, so I am using a modified
gromos54a7_atb force field, i.e., manually adding the parameters for boron
and nitrogen to the bonded & nonbonded .itp files.
The parameters are obtained from literature.

 
3. How was the topology obtained? If x2top was used, was the BN sample
in a box with in-plane PBC?


I did use x2top to generate the topology of hBN sheet by using the
following command:

g_x2top_mpi -f HBN-z-2.gro -o hbn-x2top.top -ff bnt_oplsaa -name BNT
-noparam -pbc -alldih

PBC is considered here. And some additional modifications are carried out
manually.


Best,
Xuliang
-- 
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] Periodic Molecule's Free Energy Calculation Error

[Jason Zhu]

Hello Gromacs Community,

I am trying to calculate the solvation free energy of a hBN sheet following
Justin Lemkul and Alchemistry's tutorials.
Since the hBN sheet is infinitely large, I turned the periodic molecules
flag on.
This runs all fine on one core, but when I try to run NVT in parallel (e.g.
4 ranks), the job would throw the following error:

A list of missing interactions:
LJC Pairs NB of 890400 missing 338688

---
Program gmx mdrun, VERSION 5.1.4
Source code file:
/gpfs/runtime/opt/gromacs/5.1.4/src/gromacs-5.1.4/src/gromacs/domdec/domdec_topology.cpp,
line: 242

Software inconsistency error:
Some interactions seem to be assigned multiple times
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---

Halting parallel program gmx mdrun on rank 1 out of 4
In: PMI_Abort(1, application called MPI_Abort(MPI_COMM_WORLD, 1) - process
1)

---
Program gmx mdrun, VERSION 5.1.4
Source code file:
/gpfs/runtime/opt/gromacs/5.1.4/src/gromacs-5.1.4/src/gromacs/domdec/domdec_topology.cpp,
line: 242

Software inconsistency error:
Some interactions seem to be assigned multiple times
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---

It seems like a domain decomposition error. My first thought was that the
system "explode".
However, when I check my topology and pbc condition carefully, there is no
sign of anything wrong.
I also tried NPT & PROD MD. The same error when I ran on multiple MPI
threads.

My question is: why the system could run fine on one MPI, but not if I
increased the number of MPI threads?

Any help on this issue will be really appreciated.


Here below is my .mdp file:

; RUN CONTROL—NVT
;——
define   = -DPOSRES_HBN
integrator   = sd; stochastic leap-frog integrator
nsteps   = 5000  ; 2 * 5,000 fs = 10 ps
dt   = 0.002 ; 2 fs
comm-mode= Linear; remove center of mass translation
nstcomm  = 100   ; frequency for center of mass motion removal

;——
; OUTPUT CONTROL
;——
nstxout= 0  ; don't save coordinates to .trr
nstvout= 0  ; don't save velocities to .trr
nstfout= 0  ; don't save forces to .trr
nstxout-compressed = 5000   ; xtc compressed trajectory output
every 5000 steps
compressed-x-precision = 1000   ; precision with which to write to the
compressed trajectory file
nstlog = 5000   ; update log file every 10 ps
nstenergy  = 5000   ; save energies every 10 ps
nstcalcenergy  = 100; calculate energies every 100 steps

;——
; BONDS
;——
constraint_algorithm   = lincs  ; holonomic constraints
constraints= all-bonds  ; hydrogens only are constrained
lincs_iter = 1  ; accuracy of LINCS (1 is default)
lincs_order= 4  ; also related to accuracy (4 is
default)
lincs-warnangle= 30 ; maximum angle that a bond can rotate
before LINCS will complain (30 is default)
continuation   = no ; formerly known as
'unconstrained-start' - useful for exact continuations and reruns

;——
; NEIGHBOR SEARCHING
;——
cutoff-scheme   = Verlet
ns-type = grid   ; search neighboring grid cells
nstlist = 10 ; 20 fs (default is 10)
rlist   = 1.2; short-range neighborlist cutoff (in nm)
pbc = xyz; 3D PBC
; PBC: grp is infinite
periodic-molecules = yes

;——
; ELECTROSTATICS
;——
coulombtype  = PME  ; Particle Mesh Ewald for long-range
electrostatics
rcoulomb = 1.2  ; short-range electrostatic cutoff (in nm)
ewald_geometry   = 3d   ; Ewald sum is performed in all three dimensions
pme-order= 4; interpolation order for PME (default is 4)
fourierspacing   = 0.16 ; grid spacing for FFT
ewald-rtol   = 1e-6 ; relative strength of the Ewald-shifted direct
potential at rcoulomb

;——
; VDW
;——
vdw-type= PME
rvdw= 1.2
vdw-modifier= Potential-Shift
ewald-rtol-lj   = 1e-3
lj-pme-comb-rule= Geometric
DispCorr= EnerPres

;——
; TEMPERATURE & PRESSURE COUPL
;——
tc_grps=  System
tau_t  =  0.1
ref_t  =  300
pcoupl =  no

;——
; VELOCITY GENERATION
;——

Re: [gmx-users] The Rationality of Position Restrain in Umbrella Sampling

[Jason Zhu]

Dear Justin,

Thank you for your prompt response. Your information is very helpful.

Now I understand that the position restrain is used to to mimic the
stability of a much larger structure of A-beta fibrils and make the pulling
direction coincident with the fibril axis.

The system I am studying is a complex of two same molecules (in linear
shapes). I tried to separate them and calculate the binding free energy. My
question is whether the free energy or energy curve is dependent of pulling
direction with respect to the complex orientation. If I removed the
position restrain, the two linear molecules preferred to rotate first until
their interface parallel to the pulling direction and slide away from each
other. Is it reasonable and valid to get the correct binding free energy
and energy curve? Do I need change the initial orientation to the preferred
one?

Best,
Wenpeng Zhu

On 5/2/17 9:54 AM, Jason Zhu wrote:
> Dear All,
>
> I am modeling the free binding energy between two molecules following the
> Gromacs Tutorial 3: Umbrella Sampling (
>
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html
> ).
>
> In this tutorial, the Chain_B of proteins is fixed as an immobile
reference
> by position restrain during not only the stage of Steered MD but also the
> stage of Umbrella Sampling.
>
> I am wondering what the rationality of the position restrain is for the
> calculation of free energy. Does this reduce of degree of freedom due to
> the position restrain have any artificial influence on the calculated free
> energy?
>

Please read the paper I linked from the tutorial, and references therein,
which
explain this rather unique case.  Normally one does not need additional
restraints during US.

-Justin

> Why don't we just use "COM motion removal" method to fix the COM of the
> protein complex and use pulling method to separate the Chain_A from other
> protein chains? Is it because of the excessively strong binding between
> them? If we use position restrain, do we need to make some corrections to
> the calculated binding free energy, just like the below paper?
>
> Boresch, S., Tettinger, F., Leitgeb, M., and Karplus, M. (2003) Absolute
> binding free energies: A quantitative approach for their calculation. J.
> Phys. Chem. A 107, 9535?9551
> (http://pubs.acs.org/doi/abs/10.1021/jp0217839)
>
> Best,
> Wenpeng Zhu
>

--
==

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

==========


--

2017-05-02 9:54 GMT-04:00 Jason Zhu <jasonzhu...@gmail.com>:

> Dear All,
>
> I am modeling the free binding energy between two molecules following the
> Gromacs Tutorial 3: Umbrella Sampling (http://www.bevanlab.biochem.
> vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html).
>
> In this tutorial, the Chain_B of proteins is fixed as an immobile
> reference by position restrain during not only the stage of Steered MD but
> also the stage of Umbrella Sampling.
>
> I am wondering what the rationality of the position restrain is for the
> calculation of free energy. Does this reduce of degree of freedom due to
> the position restrain have any artificial influence on the calculated free
> energy?
>
> Why don't we just use "COM motion removal" method to fix the COM of the
> protein complex and use pulling method to separate the Chain_A from other
> protein chains? Is it because of the excessively strong binding between
> them? If we use position restrain, do we need to make some corrections to
> the calculated binding free energy, just like the below paper?
>
> Boresch, S., Tettinger, F., Leitgeb, M., and Karplus, M. (2003) Absolute
> binding free energies: A quantitative approach for their calculation. J.
> Phys. Chem. A 107, 9535–9551
> (http://pubs.acs.org/doi/abs/10.1021/jp0217839)
>
> Best,
> Wenpeng Zhu
>
-- 
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] The Rationality of Position Restrain in Umbrella Sampling

[Jason Zhu]

Dear All,

I am modeling the free binding energy between two molecules following the
Gromacs Tutorial 3: Umbrella Sampling (
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html
).

In this tutorial, the Chain_B of proteins is fixed as an immobile reference
by position restrain during not only the stage of Steered MD but also the
stage of Umbrella Sampling.

I am wondering what the rationality of the position restrain is for the
calculation of free energy. Does this reduce of degree of freedom due to
the position restrain have any artificial influence on the calculated free
energy?

Why don't we just use "COM motion removal" method to fix the COM of the
protein complex and use pulling method to separate the Chain_A from other
protein chains? Is it because of the excessively strong binding between
them? If we use position restrain, do we need to make some corrections to
the calculated binding free energy, just like the below paper?

Boresch, S., Tettinger, F., Leitgeb, M., and Karplus, M. (2003) Absolute
binding free energies: A quantitative approach for their calculation. J.
Phys. Chem. A 107, 9535–9551
(http://pubs.acs.org/doi/abs/10.1021/jp0217839)

Best,
Wenpeng Zhu
-- 
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] Non-overlapping region appears in histograms when using SMD+Umbrella Sampling for strong binding complex

[Jason Zhu]

Dear Justin,

Thank you very much for your prompt response.

It is really a good idea to get the desired structures by "setting the
reference distance to the one you want, setting the pull rate to zero, and
just letting the system equilibrate under the biasing potential". I will
try what you suggest.

Only one small question about the "SMD/PMF" you are referring to in the end
of the letter. Is it the same method (SMD+Umbrella Sampling) I used, or
calculating PMF directly by integrating the force in SMD? I also tried the
latter one. But the result is dependent of the loading rate.

Best,
Wenpeng Zhu

On 4/20/17 2:30 PM, Jason Zhu wrote:
> Dear All,
>
> I am modeling the binding free energy between a molecule of sodium cholate
> and a 2D nanosheet by using steered MD and umbrella sampling as guided by
> Gromacs Tutorial 3 (
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/
gmx-tutorials/umbrella/index.html
> ).
>
> The small molecule has a strong binding with the surface of nanosheet. I
> first performed a Steered MD to pull the molecule away from the surface.
> The loading rate is low as 0.01 nm/ns, the force constant for pulling is
> 1000 kJ/nm2 and the output frequency of configuration is 1 frame very 100
> ps.
>
> However, due to the strong binding, the COM distance between the molecule
> and nanosheet has a big jump (~0.4 nm) at the detachment. I cannot extract
> successive sufficiently configurations at the moment of detachment for
> umbrella sampling (10 ns sampling). Even if I choose every configurations
> near the detachment for sampling, the histograms still have a wide gap and
> the curve in the energy profile is noncontinuous.
>
> I have tried to increase the force constant for pulling and even decreased
> the loading rate. But the problem still exists.
>
> I wonder if you have any experience to solve this problem. It would be
> appreciated if you could provide any suggestion.
>

Force reasonable starting structures by taking the closest snapshot you
have to
the desired distance, set the reference distance to the one you want, set
the
pull rate to zero, and just let the system equilibrate under the biasing
potential.  After a bit of time, you should have a reasonable starting
structure.

> Should I use the other method for free energy calculation as in the
Gromacs
> Tutorial 6 (
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/
gmx-tutorials/free_energy/index.htm),
> instead of using SMD and umbrella sampling? But there are also some
> questions when I use this method. First, the molecule is negatively
> charged. When I uncouple the electrostatic interactions, should I uncouple
> the sodium ions? If so, I am calculating the free energy of the molecule
> and its neutralizing ion, not only the molecule. Is it correct? If I am
> only interested in the molecule, what should I do? In this method, unlike
> the umbrella sampling, we only can get the free energy change, not the
> energy curve that includes the information of energy barrier. Is it
correct?
>

You're unlikely to get good convergence with this approach.  SMD/PMF is the
way
to go.

-Justin
-- 
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] Non-overlapping region appears in histograms when using SMD+Umbrella Sampling for strong binding complex

[Jason Zhu]

Dear All,

I am modeling the binding free energy between a molecule of sodium cholate
and a 2D nanosheet by using steered MD and umbrella sampling as guided by
Gromacs Tutorial 3 (
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html
).

The small molecule has a strong binding with the surface of nanosheet. I
first performed a Steered MD to pull the molecule away from the surface.
The loading rate is low as 0.01 nm/ns, the force constant for pulling is
1000 kJ/nm2 and the output frequency of configuration is 1 frame very 100
ps.

However, due to the strong binding, the COM distance between the molecule
and nanosheet has a big jump (~0.4 nm) at the detachment. I cannot extract
successive sufficiently configurations at the moment of detachment for
umbrella sampling (10 ns sampling). Even if I choose every configurations
near the detachment for sampling, the histograms still have a wide gap and
the curve in the energy profile is noncontinuous.

I have tried to increase the force constant for pulling and even decreased
the loading rate. But the problem still exists.

I wonder if you have any experience to solve this problem. It would be
appreciated if you could provide any suggestion.

Should I use the other method for free energy calculation as in the Gromacs
Tutorial 6 (
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.htm),
instead of using SMD and umbrella sampling? But there are also some
questions when I use this method. First, the molecule is negatively
charged. When I uncouple the electrostatic interactions, should I uncouple
the sodium ions? If so, I am calculating the free energy of the molecule
and its neutralizing ion, not only the molecule. Is it correct? If I am
only interested in the molecule, what should I do? In this method, unlike
the umbrella sampling, we only can get the free energy change, not the
energy curve that includes the information of energy barrier. Is it correct?

Many thanks for your helps in advance!

Best,
Wenpeng Zhu
-- 
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.