From e9b44a6bdb7bdef41999ba96c60e4d1aa3c2ebfb Mon Sep 17 00:00:00 2001 From: CalCraven Date: Thu, 16 Jul 2026 16:28:14 -0500 Subject: [PATCH] Add changes for gcpm potential loading and handling in lammps. Initial test for mixed gcpm system in lammps --- gmso/core/topology.py | 2 +- gmso/formats/lammpsdata.py | 135 ++++++++++++++++---- gmso/lib/jsons/GCPMPairPotential.json | 12 ++ gmso/tests/files/gcpm.xml | 58 +++++++++ gmso/tests/test_lammps.py | 174 ++++++++++++++++++++++++++ 5 files changed, 359 insertions(+), 22 deletions(-) create mode 100644 gmso/lib/jsons/GCPMPairPotential.json create mode 100644 gmso/tests/files/gcpm.xml diff --git a/gmso/core/topology.py b/gmso/core/topology.py index 1c2c699f2..dbde6522b 100644 --- a/gmso/core/topology.py +++ b/gmso/core/topology.py @@ -225,7 +225,7 @@ def combining_rule(self): @combining_rule.setter def combining_rule(self, rule): """Set the combining rule for the topology.""" - if rule not in ["lorentz", "geometric"]: + if rule not in ["lorentz", "geometric", "None"]: raise GMSOError("Combining rule must be `lorentz` or `geometric`") self._combining_rule = rule diff --git a/gmso/formats/lammpsdata.py b/gmso/formats/lammpsdata.py index a811128ee..6f7decad5 100644 --- a/gmso/formats/lammpsdata.py +++ b/gmso/formats/lammpsdata.py @@ -6,7 +6,7 @@ import datetime import logging import os -from itertools import count +from itertools import chain, count from pathlib import Path from typing import Optional, Union @@ -30,6 +30,7 @@ from gmso.lib.potential_templates import PotentialTemplateLibrary from gmso.utils.compatibility import check_compatibility from gmso.utils.conversions import convert_kelvin_to_energy_units +from gmso.utils.expression import NullPotentialExpression from gmso.utils.sorting import ( reindex_molecules, sort_by_types, @@ -64,7 +65,7 @@ def write_lammpsdata( Path of the output file. atom_style : str, optional, default='full' LAMMPS atom style. Supported values: ``'full'``, ``'atomic'``, - ``'charge'``, ``'molecular'``. + ``'charge'``, ``'molecular'``, ``'full + dipole + sphere'``. unit_style : str, optional, default='real' LAMMPS unit system. Supported values: ``'real'``, ``'lj'``, ``'metal'``, ``'si'``, ``'cgs'``, ``'electron'``, ``'micro'``, @@ -93,7 +94,13 @@ def write_lammpsdata( unit styles are currently supported. """ - if atom_style not in ["full", "atomic", "molecular", "charge"]: + if atom_style not in [ + "full", + "atomic", + "molecular", + "charge", + "full + dipole + sphere", + ]: raise ValueError( 'Atom style "{}" is invalid or is not currently supported'.format( atom_style @@ -596,6 +603,7 @@ def _accepted_potentials(): """List of accepted potentials that LAMMPS can support.""" templates = PotentialTemplateLibrary() lennard_jones_potential = templates["LennardJonesPotential"] + BuckinghamExp6_pair_potential = templates["GCPMPairPotential"] harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"] fene_bond_potential = templates["LAMMPSFENEBondPotential"] harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"] @@ -603,6 +611,7 @@ def _accepted_potentials(): periodic_torsion_potential = templates["PeriodicTorsionPotential"] harmonic_improper_potential = templates["HarmonicImproperPotential"] opls_torsion_potential = templates["OPLSTorsionPotential"] + null_atom_potential = NullPotentialExpression() accepted_potentialsList = [ lennard_jones_potential, harmonic_bond_potential, @@ -612,6 +621,8 @@ def _accepted_potentials(): periodic_torsion_potential, harmonic_improper_potential, opls_torsion_potential, + null_atom_potential, + BuckinghamExp6_pair_potential, ] return accepted_potentialsList @@ -656,7 +667,7 @@ def _write_header(out_file, top, atom_style, dihedral_parser): str(datetime.datetime.now()), ) ) - out_file.write("{:d} atoms\n".format(top.n_sites)) + out_file.write("{:d} atoms\n".format(top.n_sites + top.n_virtual_sites)) if atom_style in ["full", "molecular"]: out_file.write("{:d} bonds\n".format(top.n_bonds)) out_file.write("{:d} angles\n".format(top.n_angles)) @@ -676,7 +687,15 @@ def _write_header(out_file, top, atom_style, dihedral_parser): out_file.write("{:d} impropers\n\n".format(top.n_impropers)) # TODO: allow users to specify filter_by syntax - out_file.write("{:d} atom types\n".format(len(top.atom_types(filter_by=pfilter)))) + n_atomtypes = len(top.atom_types(filter_by=pfilter)) + n_atomtypes += len( + set( + { + site.virtual_type.name: site.virtual_type for site in top.virtual_sites + }.values() + ) + ) + out_file.write("{:d} atom types\n".format(n_atomtypes)) if top.n_bonds > 0 and atom_style in ["full", "molecular"]: out_file.write( "{:d} bond types\n".format(len(top.bond_types(filter_by=pfilter))) @@ -769,23 +788,76 @@ def _write_atomtypes(out_file, top, base_unyts, cfactorsDict): out_file.write("\nMasses\n") out_file.write(f"#\tmass ({base_unyts.usystem['mass']})\n") atypesView = sorted(top.atom_types(filter_by=pfilter), key=lambda x: x.name) + atypesView += sorted( + set([vsite.virtual_type for vsite in top.virtual_sites]), key=lambda x: x.name + ) for atom_type in atypesView: + if isinstance(atom_type, gmso.VirtualType): + mass = 1e-100 + else: + mass = float(base_unyts.convert_parameter(atom_type.mass, cfactorsDict)) out_file.write( - "{:d}\t{}\t# {}\n".format( - atypesView.index(atom_type) + 1, - base_unyts.convert_parameter(atom_type.mass, cfactorsDict), - atom_type.name, - ) + (f"{atypesView.index(atom_type) + 1:d}\t{mass:.6g}\t# {atom_type.name}\n") ) def _write_pairtypes(out_file, top, base_unyts, cfactorsDict): """Write out pair interaction to LAMMPS file.""" - # TODO: Handling of modified cross-interactions is not considered from top.pairpotential_types + worker_functions = { + "GCPMPairPotential": _write_gcpm_pairstyle, + "LennardJonesPotential": _write_lj_pairstyle, + } + # check for pairpotentials or atomtype potentials + if top.pairpotential_types: + first_pairtype = top.pairpotential_types[0] + return worker_functions[first_pairtype.name]( + out_file, top, first_pairtype, base_unyts, cfactorsDict + ) + elif top.atom_types: # By default try to write lj pairstyles + return _write_lj_pairstyle(out_file, top, base_unyts, cfactorsDict) + else: + raise ValueError(f"Unknown pairpotential type in {top}") + + +def _write_gcpm_pairstyle(out_file, top, first_pairtype, base_unyts, cfactorsDict): + test_pairtype = top.pairpotential_types[0] + out_file.write(f"\nPair Coeffs # {test_pairtype.expression}\n") + pair_style_orderTuple = ("epsilon", "sigma", "gamma_exp6", "alpha", "sigmai") + param_labels = [ + write_out_parameter_and_units( + key, + convert_kelvin_to_energy_units(first_pairtype.parameters[key], "kJ"), + base_unyts, + ) + for key in pair_style_orderTuple + ] + + out_file.write("#\t" + "\t".join(param_labels) + "\n") + pair_types = list(top.pairpotential_types) + pair_types.sort(key=lambda x: sorted(x.member_types)) + for idx, pair_type in enumerate(pair_types): + member_types = sorted([mem for mem in pair_type.member_types]) + out_file.write( + "{}\t{:7}\t{:7}\t\t# {}\t{}\n".format( + idx + 1, + *member_types, + *[ + base_unyts.convert_parameter( + convert_kelvin_to_energy_units(pair_type.parameters[key], "kJ"), + cfactorsDict, + n_decimals=6, + ) + for key in pair_style_orderTuple + ], + ) + ) + + +def _write_lj_pairstyle(out_file, top, base_unyts, cfactorsDict): # Pair coefficients test_atomtype = top.sites[0].atom_type out_file.write(f"\nPair Coeffs # {test_atomtype.expression}\n") - nb_style_orderTuple = ( + nb_style_orderTuple = ( # TODO: Handle different pair types "epsilon", "sigma", ) # this will vary with new pair styles @@ -1138,6 +1210,7 @@ def parse_harmonic_style_improper(improper_type): def _write_site_data(out_file, top, atom_style, base_unyts, cfactorsDict): """Write atomic positions and charges to LAMMPS file..""" out_file.write(f"\nAtoms #{atom_style}\n\n") + extraVarsDict = {} # dict for extra formatting if atom_style == "atomic": atom_line = "{index:d}\t{type_index:d}\t{x:.8}\t{y:.8}\t{z:.8}\n" elif atom_style == "charge": @@ -1148,34 +1221,54 @@ def _write_site_data(out_file, top, atom_style, base_unyts, cfactorsDict): ) elif atom_style == "full": atom_line = "{index:d}\t{moleculeid:d}\t{type_index:d}\t{charge:.8}\t{x:.8}\t{y:.8}\t{z:.8}\n" + elif atom_style == "full + dipole + sphere": + extraVarsDict = { + "mux": 0.0, + "muy": 0.0, + "muz": 0.0, + "radius": 0.0, + "rmass": 1.0, + } + atom_line = "{index:d}\t{type_index:d}\t{x}\t{y}\t{z}\t{moleculeid:d}\t{charge:.8}\t{mux:.8}\t{muy:.8}\t{muz:.8}\t{radius:.8}\t{rmass:.8}\n" unique_sorted_typesList = sorted( top.atom_types(filter_by=pfilter), key=lambda x: x.name ) - for i, site in enumerate(top.sites): + unique_sorted_typesList += sorted( # append virtual types + set(site.virtual_type for site in top.virtual_sites), key=lambda x: x.name + ) + for i, site in enumerate(chain(top._sites, top._virtual_sites)): + if isinstance(site, gmso.Atom): + position = site.position + type_index = unique_sorted_typesList.index(site.atom_type) + 1 + elif isinstance(site, gmso.VirtualSite): + position = site.position() # callable function + type_index = unique_sorted_typesList.index(site.virtual_type) + 1 + # set muz to 0.01 to allow for polarizeability of site + extraVarsDict["muz"] = 0.01 + out_file.write( atom_line.format( index=i + 1, moleculeid=site.molecule.number + 1, # index is 0-based in GMSO - type_index=unique_sorted_typesList.index(site.atom_type) + 1, + type_index=type_index, charge=base_unyts.convert_parameter( site.charge, cfactorsDict, n_decimals=6, ), x=base_unyts.convert_parameter( - site.position[0], + position[0], cfactorsDict, n_decimals=6, ), - y=base_unyts.convert_parameter( - site.position[1], cfactorsDict, n_decimals=6 - ), - z=base_unyts.convert_parameter( - site.position[2], cfactorsDict, n_decimals=6 - ), + y=base_unyts.convert_parameter(position[1], cfactorsDict, n_decimals=6), + z=base_unyts.convert_parameter(position[2], cfactorsDict, n_decimals=6), + **extraVarsDict, ) ) + if extraVarsDict.get("muz"): # reset muz for non virtual types + extraVarsDict["muz"] = 0.0 def _angle_order_sorter(angle_typesList): diff --git a/gmso/lib/jsons/GCPMPairPotential.json b/gmso/lib/jsons/GCPMPairPotential.json new file mode 100644 index 000000000..cca0ef0c9 --- /dev/null +++ b/gmso/lib/jsons/GCPMPairPotential.json @@ -0,0 +1,12 @@ +{ + "name": "GCPMPairPotential", + "expression": "epsilon/(1-6/gamma_exp6)*(6/gamma_exp6*exp(gamma_exp6*(1-rOO/sigma))-(sigma/rOO)**6)+alpha+sigmai", + "independent_variables": "rOO", + "expected_parameters_dimensions": { + "epsilon": "energy", + "sigma": "length", + "gamma_exp6": "dimensionless", + "alpha": "length**3", + "sigmai": "length" + } +} diff --git a/gmso/tests/files/gcpm.xml b/gmso/tests/files/gcpm.xml new file mode 100644 index 000000000..bc4e01f22 --- /dev/null +++ b/gmso/tests/files/gcpm.xml @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/gmso/tests/test_lammps.py b/gmso/tests/test_lammps.py index 1bb7062f2..737770547 100644 --- a/gmso/tests/test_lammps.py +++ b/gmso/tests/test_lammps.py @@ -621,3 +621,177 @@ def test_lammps_fene(self, typed_ethane): np.testing.assert_allclose(float(coeffs[2]), 10) np.testing.assert_allclose(float(coeffs[3]), 1) np.testing.assert_allclose(float(coeffs[4]), 10) + + def test_gcpm(self): + import mbuild as mb + + from gmso.parameterization import apply + + posO = np.array([2.5, 2.5, 2.5]) + theta = 104.52 / 180 * np.pi # rad + dOH = 0.09572 + dM = 0.027 + posH1 = np.array([dOH, 0, 0]) + posO + posH2 = np.array([dOH * np.cos(theta), dOH * np.sin(theta), 0]) + posO + posM = np.array([dM * np.cos(theta / 2), dM * np.sin(theta / 2), 0]) + posO + + # posO = np.array([1.9397813577, 5.6972356033, 22.4183552067]) / 10 # nm + # posH1 = np.array([1.6442699843, 4.9565922701, 21.8888721369]) / 10 + # posH2 = np.array([1.6248397438, 6.4692075712, 21.9481394044]) / 10 + # posM = np.array([1.7991197323, 5.7044544011, 22.1880028763]) / 10 + + cpd = mb.Compound(name="water") + particle1 = mb.Compound(name="O", element="O", pos=posO) + particle2 = mb.Compound(name="H", element="H", pos=posH1) + particle3 = mb.Compound(name="H", element="H", pos=posH2) + for part in [particle1, particle2, particle3]: + cpd.add(part) + cpd.add_bond((particle1, particle2)) + cpd.add_bond((particle1, particle3)) + cpd.box = mb.Box([5, 5, 5]) + + top = cpd.to_gmso() + top.set_rigid("water") + ff = gmso.ForceField(get_path("gcpm.xml")) + + ptop = apply(top, ff, ignore_params=["angle", "bond"]) + assert ptop.n_virtual_sites == 1 + assert len(ptop.pairpotential_types) == 3 + assert ptop.pairpotential_types[0].name == "GCPMPairPotential" + assert ptop.n_sites == 3 + + v_site = ptop.virtual_sites[0] + assert v_site.charge.value == -1.226 + assert np.allclose( + ptop.sites[0].position, posO + ) # hand validate position of O-site + assert np.allclose( + ptop.sites[1].position, posH1 + ) # hand validate position of H1-site + assert np.allclose( + ptop.sites[2].position, posH2 + ) # hand validate position of H2-site + assert np.allclose(v_site.position(), posM) # hand validate position of M-site + ptop.save( + "test.lammps", atom_style="full + dipole + sphere" + ) # should error, maybe should be gcpm + + expected_file = [ + "\n", + "\n", + "4 atoms\n", + "3 atom types\n", + "\n", + "0.000000 50.000000 xlo xhi\n", + "0.000000 50.000000 ylo yhi\n", + "0.000000 50.000000 zlo zhi\n", + "0.000000 0.000000 0.000000 xy xz yz\n", + "\n", + "Masses\n", + "#\tmass (amu)\n", + "1\t15.999\t# 1\n", + "2\t1.008\t# 2\n", + "3\t1e-100\t# 3\n", + "\n", + "Pair Coeffs # alpha + epsilon*(-sigma**6/rOO**6 + 6*exp(gamma_exp6*(-rOO/sigma + 1))/gamma_exp6)/(1 - 6/gamma_exp6) + sigmai\n", + "#\tepsilon (kcal/mol)\tsigma (Å)\tgamma_exp6 (dimensionless)\talpha (Å**3)\tsigmai (Å)\n", + "1\t1 \t1 \t\t# 0.218592\t3.690000\n", + "2\t2 \t2 \t\t# 0.000000\t1.000000\n", + "3\t3 \t3 \t\t# 0.000000\t1.000000\n", + "\n", + "Atoms #full + dipole + sphere\n", + "\n", + f"1\t1\t{posO[0] * 10:.6f}\t{posO[1] * 10:.6f}\t{round(posO[2] * 10, 6):.6f}\t1\t0.000000\t0.0\t0.0\t0.0\t0.0\t1.0\n", + f"2\t2\t{posH1[0] * 10:.6f}\t{posH1[1] * 10:.6f}\t{posH1[2] * 10:.6f}\t1\t0.611300\t0.0\t0.0\t0.0\t0.0\t1.0\n", + f"3\t2\t{posH2[0] * 10:.6f}\t{posH2[1] * 10:.6f}\t{posH2[2] * 10:.6f}\t1\t0.611300\t0.0\t0.0\t0.0\t0.0\t1.0\n", + f"4\t3\t{posM[0] * 10:.6f}\t{posM[1] * 10:.6f}\t{posM[2] * 10:.6f}\t1\t-1.22600\t0.0\t0.0\t0.01\t0.0\t1.0\n", + ] + with open("test.lammps", "r") as f: + readlines = f.readlines() + for l1, l2 in zip(readlines[1:], expected_file): + assert l1.strip() == l2.strip(), (l1, l2) + + def test_box_gcpm(self): + import mbuild as mb + + from gmso.parameterization import apply + + posO = np.array([0, 0, 0]) + theta = 104.52 / 180 * np.pi # rad + dOH = 0.09572 + # dM = 0.027 + posH1 = np.array([dOH, 0, 0]) + posH2 = np.array([dOH * np.cos(theta), dOH * np.sin(theta), 0]) + + posO = np.array([1.9397813577, 5.6972356033, 22.4183552067]) / 10 # nm + posH1 = np.array([1.6442699843, 4.9565922701, 21.8888721369]) / 10 + posH2 = np.array([1.6248397438, 6.4692075712, 21.9481394044]) / 10 + + cpd = mb.Compound(name="water") + particle1 = mb.Compound(name="O", element="O", pos=posO) + particle2 = mb.Compound(name="H", element="H", pos=posH1) + particle3 = mb.Compound(name="H", element="H", pos=posH2) + for part in [particle1, particle2, particle3]: + cpd.add(part) + cpd.add_bond((particle1, particle2)) + cpd.add_bond((particle1, particle3)) + box = mb.Box([10, 10, 10]) + n_molecules = 2 + filled_box = mb.fill_box(cpd, n_molecules, box) + + top = filled_box.to_gmso() + top.set_rigid("water") + ff = gmso.ForceField(get_path("gcpm.xml")) + + ptop = apply(top, ff, ignore_params=["angle", "bond"]) + assert ptop.n_virtual_sites == n_molecules + assert len(ptop.pairpotential_types) == 3 + assert ptop.n_sites == 3 * n_molecules + assert ptop.pairpotential_types[0].name == "GCPMPairPotential" + + ptop.save("test.lammps", atom_style="full + dipole + sphere") + assert os.path.exists("./test.lammps") # check that file was correctly run + + def test_lammps_gcpm_mixture(self, capsys): + import textwrap + + import mbuild as mb + + from gmso.parameterization import apply + + # create system with gcpm water and ethanol + water = mb.load("O", smiles=True) + water.name = "Water" + ethanol = mb.load("CCO", smiles=True) + ethanol.name = "Ethanol" + box = mb.fill_box([water, ethanol], [1, 1], box=[3, 3, 3]) + # forcefield with two different forcefields + ffgcpm = gmso.ForceField(get_path("gcpm.xml")) + ffgcpm.combining_rule = "geometric" + ffopls = gmso.ForceField("oplsaa") + + top = box.to_gmso() + top.set_rigid("Water") + ptop = apply( + top, + {"Water": ffgcpm, "Ethanol": ffopls}, + identify_connections=True, + ignore_params=["bond", "angle", "dihedral", "improper"], + ) + ptop.save("out.lammps") + # Should we just print out pair_styles for lammps input script? + output_script = """ + pair_style hybrid gcpm lj/cut/coul/cut {r_cut} + pair_modify mix {mixing} #arithmetic or geometric + pair_coeff 1 1 gcpm 0.218592 3.690000 + pair_coeff 2 2 gcpm 0.000000 1.000000 + pair_coeff 3 3 gcpm 0.000000 1.000000 + pair_coeff 4 4 lj/cut/coul/cut 0.276144 0.35 + pair_coeff 5 5 lj/cut/coul/cut 0.276144 0.35 + pair_coeff 6 6 lj/cut/coul/cut 0.71128 0.312 + pair_coeff 7 7 lj/cut/coul/cut 0.0 1.0 + pair_coeff 8 8 lj/cut/coul/cut 0.12552 0.25\n""" + output_script = textwrap.dedent(output_script).strip() + print(output_script) + captured_script = capsys.readouterr() + assert output_script + "\n" == captured_script.out