[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