diff --git a/mbuild/formats/gsdwriter.py b/mbuild/formats/gsdwriter.py index c67620c0..92d00252 100644 --- a/mbuild/formats/gsdwriter.py +++ b/mbuild/formats/gsdwriter.py @@ -17,7 +17,8 @@ def write_gsd(structure, filename, ref_distance=1.0, ref_mass=1.0, - ref_energy=1.0, rigid_bodies=None, shift_coords=True): + ref_energy=1.0, rigid_bodies=None, shift_coords=True, + write_special_pairs=True): """Output a GSD file (HOOMD v2 default data format). Parameters @@ -39,6 +40,9 @@ def write_gsd(structure, filename, ref_distance=1.0, ref_mass=1.0, part of a rigid body. shift_coords : bool, optional, default=True Shift coordinates from (0, L) to (-L/2, L/2) if necessary. + write_special_pairs : bool, optional, default=True + Writes out special pair information necessary to correctly use the OPLS fudged 1,4 interactions + in HOOMD. Notes ----- @@ -79,6 +83,8 @@ def write_gsd(structure, filename, ref_distance=1.0, ref_mass=1.0, _write_particle_information(gsd_file, structure, xyz, ref_distance, ref_mass, ref_energy, rigid_bodies) + if write_special_pairs: + _write_pair_information(gsd_file, structure) if structure.bonds: _write_bond_information(gsd_file, structure) if structure.angles: @@ -86,6 +92,7 @@ def write_gsd(structure, filename, ref_distance=1.0, ref_mass=1.0, if structure.rb_torsions: _write_dihedral_information(gsd_file, structure) + gsd.hoomd.create(filename, gsd_file) def _write_particle_information(gsd_file, structure, xyz, ref_distance, @@ -124,6 +131,33 @@ def _write_particle_information(gsd_file, structure, xyz, ref_distance, rigid_bodies = [-1 if body is None else body for body in rigid_bodies] gsd_file.particles.body = rigid_bodies +def _write_pair_information(gsd_file, structure): + """Write the special pairs in the system. + + Parameters + ---------- + gsd_file : + The file object of the GSD file being written + structure : parmed.Structure + Parmed structure object holding system information + """ + pair_types = [] + pair_typeid = [] + pairs = [] + for ai in structure.atoms: + for aj in ai.dihedral_partners: + #make sure we don't double add + if ai.idx > aj.idx: + ps = '-'.join(sorted([ai.type, aj.type], key=natural_sort)) + if ps not in pair_types: + pair_types.append(ps) + pair_typeid.append(pair_types.index(ps)) + pairs.append((ai.idx, aj.idx)) + gsd_file.pairs.types = pair_types + gsd_file.pairs.typeid = pair_typeid + gsd_file.pairs.group = pairs + gsd_file.pairs.N = len(pairs) + def _write_bond_information(gsd_file, structure): """Write the bonds in the system. diff --git a/mbuild/tests/test_gsd.py b/mbuild/tests/test_gsd.py index 7d867a98..6b474d0e 100755 --- a/mbuild/tests/test_gsd.py +++ b/mbuild/tests/test_gsd.py @@ -190,6 +190,26 @@ def test_bonded(self, ethane): assert np.array_equal(dihedral_atoms, expected_dihedral_atoms) assert np.array_equal(dihedral_typeids, np.zeros(n_dihedrals)) + @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") + def test_pairs(self, benzene): + from foyer import Forcefield + import gsd, gsd.pygsd + + benzene.save(filename='benzene.gsd', forcefield_name='oplsaa') + gsd_file = gsd.pygsd.GSDFile(open('benzene.gsd', 'rb')) + frame = gsd.hoomd.HOOMDTrajectory(gsd_file).read_frame(0) + + structure = benzene.to_parmed() + forcefield = Forcefield(name='oplsaa') + structure = forcefield.apply(structure) + + # Pairs + assert len(frame.pairs.types) == 3 + assert frame.pairs.N == 21 + + + + @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") def test_units(self, ethane): import gsd, gsd.pygsd