Re: [gmx-users] Periodic Molecule's Free Energy Calculation Error

2017-07-14 Thread Alex

Hi Jason,

I didn't have time to look at all of those papers, but look for instance 
at eq. (2) in Hilder et al's paper and the definition of the f_dmp 
function on the right -- this tapers off short-range interactions for 
close neighbors. I am not sure if this particular function results in 
any serious difference (I think, given the value of R_r, it was designed 
pretty much exactly to remove all nonbonded interactions closer than 
second neighbor), but the point is that Gromacs doesn't use anything 
like this. As a result, many things cannot be translated into Gromacs 
FFs "as is," requiring a fairly broad range of effort from mild 
tinkering all the way to complete reparameterization.. Even then the 
best one can hope for with Gromacs when applied to the actual properties 
of crystals is structure stability and correct bond lengths -- nothing 
beyond that. This is really more of a general comment, because more and 
more people are trying to simulate solid-state structures with Gromacs.


Best of luck!

Alex


On 7/13/2017 6:43 PM, Jason Zhu wrote:

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 <nedoma...@gmail.com <mailto:nedoma...@gmail.com>>
To: Discussion list for GROMACS users <gmx-us...@gromacs.org 
<mailto:gmx-us...@gromacs.org>>

Subject: Re: [gmx-users] Periodic Molecule's Free Energy Calculation
Error
Message-ID: <4b927ef2-9a56-6655-8053-d284cf88b...@gmail.com 
<mailto: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 

Re: [gmx-users] Periodic Molecule's Free Energy Calculation Error

2017-07-13 Thread 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 <nedoma...@gmail.com>
To: Discussion list for GROMACS users <gmx-us...@gromacs.org>
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

2017-07-13 Thread 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

2017-07-13 Thread Justin Lemkul



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

2017-07-12 Thread Alex




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

2017-07-12 Thread 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.


Re: [gmx-users] Periodic Molecule's Free Energy Calculation Error

2017-07-11 Thread Alex
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?


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


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


Alex


On 7/11/2017 8:39 PM, Jason Zhu wrote:

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


[gmx-users] Periodic Molecule's Free Energy Calculation Error

2017-07-11 Thread 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
;——