mBuild fill box not able to save to gro file via parmed or mdtraj

[rsdefever]

Bug summary

I am trying to create a system and save a .gro file with mbuild. The process appears to work OK until the saving step. In the case of saving through parmed, it crashes part way through writing the file. If saving through mdtraj via system.to_trajectory().save("test.gro"), it completes saving the file but the .gro file is corrupted.

Code to reproduce the behavior

import mbuild

k = mbuild.Compound(name="K")
li = mbuild.Compound(name="Li")
cl = mbuild.Compound(name="Cl")

# Taken from Wang et al PIM results
# at 1 bar and 1100 K
licl_density=1431 #kg/m^3 
kcl_density=1497 #kg/m^3 

# Use 1000 total ions
# x_li = 0.5
n_ions = 1000
x_li = 0.5

# Compute mixture density with ideal mixing
mixed_density = x_li * licl_density + (1.0-x_li) * kcl_density

# Compute the number of ions of each type
n_cl = int(0.5 * n_ions)
n_li = int(0.5 * x_li * n_ions)
n_k = n_ions - n_cl - n_li

# Create the system
system = mbuild.fill_box(
    [k, li, cl],
    n_compounds=[n_k, n_li, n_cl],
    density=mixed_density
)

system.to_trajectory().save("test-mdtraj.gro")
system.save("test-pmd.gro", overwrite=True)

Software versions

• Which version of mBuild are you using? (python -c "import mbuild as mb; print(mb.version)") 8a28926
• Which version of Python (python --version)? 3.7.6
• Which operating system? MacOS

[daico007]

I will look into the packing functions and see if I can fix it.

[rsdefever]

@daico007 thanks! I'm surprised that such a simple case is causing problems...the error from parmed and the corruption in the mdtraj gro file seem to perhaps point towards residue numbering? But I'm unsure...

[rsdefever]

@daico007 update: writing to .gro from MDTraj seems to be working...though I'm not sure what changed from earlier...parmed is still throwing the same error.

[ahy3nz]

system.save("test-pmd.gro", overwrite=True, combine='all')

Parmed, by default, does some logic to try to re-index atoms to keep molecules/residues together. Passing the kwarg for combine, I was able to get the mdtraj and pmd gro file to have the same contents, besides some residue off-by-one indexing and atom off-by-one indexing

[rsdefever]

@ahy3nz thanks for looking into this and for the pointer. Is this documented anywhere or is there anything we can do to try to prevent this issue from cropping up?

[ahy3nz]

https://github.com/mosdef-hub/mbuild/issues/525

[bc118]

The error with the system.save("test-pmd_combine_None.gro", overwrite=True, combine=None ) function, appears to be in the parmed folder.
Path that the system.save() function takes to to the parmed function:

1. to the mbuild.compound.py file and conversion.save() (line 1857)
2. to the mbuild.conversion.py file and structure.save()(line 867)
3. to the parmed.structure.py file and gromacs.GromacsGroFile.write() (line 1484)
4. to the parmed.gromacs.gromacsgro.py file and _write_atom_line() (line 225)

In the parmed.gromacs.gromacsgro.py file and _write_atom_line() function, there are 2 paths to take:

1. combine =’all’; which allows the atoms to be written in order and does not have any if/else statements adding complexity and problems. This is likely why it works, as opposed to the None option
2. combine = None (i.e., combine !=’all’): it has many if statements, one of which only allows the next atom to be written if has the same molecule type, name and residue as the original atom. The real failure occurs in lines 269 to 272 in the _write_atom_line() function, where it looks the original molecule and checks for the same type, name, and residue.

This is quite odd and not explained well in their description below: “combine : 'all', None, or list of iterables, optional Equivalent to the combine argument of the GromacsTopologyFile.write method. If None, system atom order may be changed to meet the need for contiguously bonded groups of atoms to be part of a single moleculetype. All other values leave the atom order unchanged. Default is None.”

Please let me know if you want me to fix this functionality in parmed, or fix it in the new GMSO, or other. It appear they may need and else statement in there if it is not the same type, name and residue. Make changes to parmed code, help writing, in MoSDeF notes… ?

[bc118]

Additional notes:

1. when using combine =’all’: the residue number for all molecules is the same (=1)

2)when using combine = None (i.e., combine !=’all’): the residue number for every molecule is different (+1 in order)

[bc118]

I built the water and salt system below, and it seems to run fine despite the atom and residue ID indexing starting at 0 instead of 1 (this indexing is fixable in mdtraj). The following can be changed in the mdtraj file, or changed before the input to this function to change the file indexing:

• mdtraj.formats.gro file and the _write_frame() function (line 410). The residue.resSeq, and serial values need to have a +1 added to them to fix to an index starting at 1.

However, this is the first time I used GROMACS, so I'm not 100% confident. The system builder (build.py) then GROMACS energy minimization files are included below, so please advise:

water_NaCl.zip

[CalCraven]

Is this still an issue? I haven't seen it when writing out .gro files, but not sure of any specific PR's that addressed it.

[bc118]

I reran my example here (from forever ago). It seems the .gro file is still starting with a 0 index and the others starting with 1