Re: [gmx-users] Periodic Molecule's Free Energy Calculation Error
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
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
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
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
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
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
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
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 ;——