Source code for gmso.formats.lammpsdata

"""Read and write LAMMPS data files."""

from __future__ import division

import copy
import datetime
import os
import warnings
from itertools import count
from pathlib import Path

import numpy as np
import unyt as u
from unyt.array import allclose_units

import gmso
from gmso.abc.abstract_site import Molecule
from gmso.core.angle import Angle
from gmso.core.atom import Atom
from gmso.core.atom_type import AtomType
from gmso.core.bond import Bond
from gmso.core.box import Box
from gmso.core.dihedral import Dihedral
from gmso.core.element import element_by_mass
from gmso.core.improper import Improper
from gmso.core.topology import Topology
from gmso.core.views import PotentialFilters
from gmso.formats.formats_registry import loads_as, saves_as
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.sorting import (
    reindex_molecules,
    sort_by_types,
    sort_connection_members,
)
from gmso.utils.units import LAMMPS_UnitSystems, write_out_parameter_and_units

pfilter = PotentialFilters.UNIQUE_SORTED_NAMES


# TODO: Write in header of each potential type any conversions that happened
# TODO: write in file header the source of the xml
[docs] @saves_as(".lammps", ".lammpsdata", ".data") def write_lammpsdata( top, filename, atom_style="full", unit_style="real", strict_potentials=False, strict_units=False, lj_cfactorsDict=None, ): """Output a LAMMPS data file. Outputs a LAMMPS data file in the 'full' atom style format. Assumes use of 'real' units. See http://lammps.sandia.gov/doc/atom_style.html for more information on atom styles. Parameters ---------- Topology : `Topology` A Topology Object filename : str Path of the output file atom_style : str, optional, default='full' Defines the style of atoms to be saved in a LAMMPS data file. The following atom styles are currently supported: 'full', 'atomic', 'charge', 'molecular' see http://lammps.sandia.gov/doc/atom_style.html for more information on atom styles. unit_style : str, optional, default='real' Can be any of "real", "lj", "metal", "si", "cgs", "electron", "micro", "nano". Otherwise an error will be thrown. These are defined in _unit_style_factory. See https://docs.lammps.org/units.html for LAMMPS documentation. strict_potentials : bool, optional, default False Tells the writer how to treat conversions. If False, then check for conversions to usable potential styles found in default_parameterMaps. If True, then error if potentials are not compatible. strict_units : bool, optional, default False Tells the writer how to treat unit conversions. If False, then check for conversions to unit styles defined in _unit_style_factory. If True, then error if parameter units do not match. lj_cfactorsDict : (None, dict), optional, default None If using unit_style="lj" only, can pass a dictionary with keys of ("mass", "energy", "length", "charge"), or any combination of these, and they will be used to non- dimensionalize all values in the topology. If any key is not passed, default values will be pulled from the topology (see _default_lj_val). These are the largest: sigma, epsilon, atomtype.mass, and atomtype.charge from the topology. Notes ----- See http://lammps.sandia.gov/doc/2001/data_format.html for a full description of the LAMMPS data format. This is a work in progress, as only a subset of everything LAMMPS supports is currently available. However, please raise issues as the current writer has been set up to eventually grow to support all LAMMPS styles. Some of this function has been adopted from `mdtraj`'s support of the LAMMPSTRJ trajectory format. See https://github.com/mdtraj/mdtraj/blob/master/mdtraj/formats/lammpstrj.py for details. """ if atom_style not in ["full", "atomic", "molecular", "charge"]: raise ValueError( 'Atom style "{}" is invalid or is not currently supported'.format( atom_style ) ) if unit_style not in [ "real", "lj", "metal", "si", "cgs", "electron", "micro", "nano", ]: raise ValueError( 'Unit style "{}" is invalid or is not currently supported'.format( unit_style ) ) if unit_style != "lj" and lj_cfactorsDict: raise ValueError("lj_cfactorsDict argument is only used if unit_style is lj.") base_unyts = LAMMPS_UnitSystems(unit_style) default_parameterMaps = { # TODO: sites are not checked currently because gmso # doesn't store pair potential eqn the same way as the connections. "impropers": [ "HarmonicImproperPotential", "HarmonicTorsionPotential", "PeriodicTorsionPotential", ], "dihedrals": ["OPLSTorsionPotential", "PeriodicTorsionPotential"], "angles": ["LAMMPSHarmonicAnglePotential"], "bonds": ["LAMMPSHarmonicBondPotential"], # "sites":"LennardJonesPotential", # "sites":"CoulombicPotential" } # TODO: Use strict_x, (e.g. x=bonds) to validate what topology attrs to convert if not strict_potentials: _try_default_potential_conversions(top, default_parameterMaps) potentialsMap = _validate_potential_compatibility(top) potential_typesDict = {} for potential in potentialsMap: pot_container = potential.__class__.__name__ potStr = pot_container.lower() + "s" potStr = potStr[:-5] + "_" + potStr[-5:] if not potential_typesDict.get(potStr): potential_typesDict[potStr] = {potentialsMap[potential]} else: potential_typesDict[potStr].add(potentialsMap[potential]) dihedral_parser = _identify_dihedral_parser(top, potential_typesDict) improper_parser = _identify_improper_parser(top, potential_typesDict) if strict_units: _validate_unit_compatibility(top, base_unyts) else: if base_unyts and unit_style != "lj": lj_cfactorsDict = None else: # LJ unit styles if lj_cfactorsDict is None: lj_cfactorsDict = {} source_factorsList = list(lj_cfactorsDict.keys()) defaultsList = ["length", "energy", "mass", "charge"] for source_factor in defaultsList + source_factorsList: if source_factor not in defaultsList: raise ValueError( f"Conversion factor {source_factor} is not used. Pleas only provide some of {defaultsList}" ) if lj_cfactorsDict.get(source_factor): continue default_val_from_topology = _default_lj_val(top, source_factor) lj_cfactorsDict[source_factor] = lj_cfactorsDict.get( source_factor, default_val_from_topology ) reindex_molecules(top) # reset the topology molecule index to match with lammps path = Path(filename) if not path.parent.exists(): msg = "Provided path to file that does not exist" raise FileNotFoundError(msg) with open(path, "w") as out_file: _write_header(out_file, top, atom_style, dihedral_parser) _write_box(out_file, top, base_unyts, lj_cfactorsDict) all_ordered_typesDict = {} if top.is_fully_typed(): _write_atomtypes(out_file, top, base_unyts, lj_cfactorsDict) _write_pairtypes(out_file, top, base_unyts, lj_cfactorsDict) if top.bond_types: sorted_bondsList = _write_bondtypes( out_file, top, base_unyts, lj_cfactorsDict ) all_ordered_typesDict["bonds"] = sorted_bondsList if top.angle_types: sorted_anglesList = _write_angletypes( out_file, top, base_unyts, lj_cfactorsDict ) all_ordered_typesDict["angles"] = sorted_anglesList if top.dihedral_types: sorted_dihedralsList = ( _write_dihedraltypes( # return a list of dihedraltypes out_file, top, base_unyts, dihedral_parser, lj_cfactorsDict, ) ) all_ordered_typesDict["dihedrals"] = sorted_dihedralsList if top.improper_types: sorted_impropersList = _write_impropertypes( out_file, top, base_unyts, improper_parser, lj_cfactorsDict ) all_ordered_typesDict["impropers"] = sorted_impropersList _write_site_data(out_file, top, atom_style, base_unyts, lj_cfactorsDict) for conn in ["bonds", "angles", "dihedrals", "impropers"]: connIter = getattr(top, conn) conn_typesList = all_ordered_typesDict.get(conn) if connIter and conn_typesList: _write_conn_data(out_file, top, conn, conn_typesList)
@loads_as(".lammps", ".lammpsdata", ".data") def read_lammpsdata( filename, atom_style="full", unit_style="real", ): """Read in a lammps data file as a GMSO topology. Parameters ---------- filename : str LAMMPS data file atom_style : str, optional, default='full' Inferred atom style defined by LAMMPS, be certain that this is provided accurately. unit_style : str, optional, default='real LAMMPS unit style used for writing the datafile. Can be "real", "lj", "metal", "si", "cgs", "electron", "micro", "nano". Returns ------- top : GMSO Topology A GMSO Topology object Notes ----- See http://lammps.sandia.gov/doc/2001/data_format.html for a full description of the LAMMPS data format. This is a work in progress, as only several atom styles, potential styles, and unit styes are currently supported. Currently supporting the following atom styles: 'full' Currently supporting the following unit styles: 'real', "real", "lj", "metal", "si", "cgs", "electron", "micro", "nano". Currently supporting the following potential styles: 'lj' Currently supporting the following bond styles: 'harmonic' Currently supporting the following angle styles: 'harmonic' Currently supporting the following dihedral styles: 'opls' Currently supporting the following improper styles: 'harmonic' """ top = Topology() # Validate 'atom_style' if atom_style not in ["full"]: raise ValueError( 'Atom Style "{}" is invalid or is not currently supported'.format( atom_style ) ) # Validate 'unit_style' if unit_style not in [ "real", "lj", "metal", "si", "cgs", "electron", "micro", "nano", ]: raise ValueError( 'Unit Style "{}" is invalid or is not currently supported'.format( unit_style ) ) base_unyts = LAMMPS_UnitSystems(unit_style) # Parse box information _get_box_coordinates(filename, base_unyts, top) # Parse atom type information top, type_list = _get_ff_information(filename, base_unyts, top) # Parse atom information _get_atoms(filename, top, base_unyts, type_list) # Parse connection (bonds, angles, dihedrals, impropers) information # TODO: Add more atom styles if atom_style in ["full"]: _get_connection(filename, top, base_unyts, connection_type="bond") _get_connection(filename, top, base_unyts, connection_type="angle") _get_connection(filename, top, base_unyts, connection_type="dihedral") _get_connection(filename, top, base_unyts, connection_type="improper") top.update_topology() return top def get_units(base_unyts, dimension): """Get u.Unit for specific LAMMPS unit style with given dimension.""" # Need separate angle units for harmonic force constant and angle if base_unyts.usystem.name == "lj": if dimension == "angle": return u.radian return u.dimensionless if dimension == "angle_eq": return u.degree # LAMMPS specifies different units for some angles, such as equilibrium angles return u.Unit(base_unyts.usystem[dimension], registry=base_unyts.reg) def _get_connection(filename, topology, base_unyts, connection_type): """Parse connection types.""" # TODO: check for other connection types besides the defaults with open(filename, "r") as lammps_file: types = False for i, line in enumerate(lammps_file): if connection_type in line.split(): n_connection_types = int(line.split()[0]) types = True if connection_type.capitalize() in line.split(): break if types is False: return topology templates = PotentialTemplateLibrary() connection_type_lines = open(filename, "r").readlines()[ i + 2 : i + n_connection_types + 2 ] connection_type_list = list() for line in connection_type_lines: if connection_type == "bond": template_potential = templates["LAMMPSHarmonicBondPotential"] # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term conn_params = { "k": float(line.split()[1]) * get_units(base_unyts, "energy") / get_units(base_unyts, "length") ** 2 * 2, "r_eq": float(line.split()[2]) * get_units(base_unyts, "length"), } name = template_potential.name expression = template_potential.expression variables = template_potential.independent_variables c_type = getattr(gmso, "BondType")( name=name, parameters=conn_params, expression=expression, independent_variables=variables, ) elif connection_type == "angle": template_potential = templates["LAMMPSHarmonicAnglePotential"] # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term conn_params = { "k": float(line.split()[1]) * get_units(base_unyts, "energy") / get_units(base_unyts, "angle") ** 2 * 2, "theta_eq": float(line.split()[2]) * get_units(base_unyts, "angle_eq"), } name = template_potential.name expression = template_potential.expression variables = template_potential.independent_variables c_type = getattr(gmso, "AngleType")( name=name, parameters=conn_params, expression=expression, independent_variables=variables, ) elif connection_type == "dihedral": template_potential = templates["OPLSTorsionPotential"] conn_params = { "k1": float(line.split()[1]) * get_units(base_unyts, "energy"), "k2": float(line.split()[2]) * get_units(base_unyts, "energy"), "k3": float(line.split()[3]) * get_units(base_unyts, "energy"), "k4": float(line.split()[4]) * get_units(base_unyts, "energy"), } name = template_potential.name expression = template_potential.expression variables = template_potential.independent_variables c_type = getattr(gmso, "DihedralType")( name=name, parameters=conn_params, expression=expression, independent_variables=variables, ) elif connection_type == "improper": template_potential = templates["HarmonicImproperPotential"] conn_params = { "k": float(line.split()[2]) * get_units(base_unyts, "energy") / get_units(base_unyts, "energy") ** 2 * 2, "phi_eq": float(line.split()[3]) * get_units(base_unyts, "angle_eq"), } name = template_potential.name expression = template_potential.expression variables = template_potential.independent_variables c_type = getattr(gmso, "ImproperType")( name=name, parameters=conn_params, expression=expression, independent_variables=variables, ) connection_type_list.append(c_type) with open(filename, "r") as lammps_file: for i, line in enumerate(lammps_file): if connection_type + "s" in line.split(): n_connections = int(line.split()[0]) if connection_type.capitalize() + "s" in line.split(): break connection_lines = open(filename, "r").readlines()[i + 2 : i + n_connections + 2] # Determine number of sites to generate if connection_type == "bond": n_sites = 2 elif connection_type == "angle": n_sites = 3 else: n_sites = 4 for i, line in enumerate(connection_lines): site_list = list() for j in range(n_sites): site = topology.sites[int(line.split()[j + 2]) - 1] site_list.append(site) ctype = copy.copy(connection_type_list[int(line.split()[1]) - 1]) ctype.member_types = tuple(map(lambda x: x.atom_type.name, site_list)) ctype.member_classes = ctype.member_types if connection_type == "bond": connection = Bond( connection_members=site_list, bond_type=ctype, ) elif connection_type == "angle": connection = Angle( connection_members=site_list, angle_type=ctype, ) elif connection_type == "dihedral": connection = Dihedral( connection_members=site_list, dihedral_type=ctype, ) elif connection_type == "improper": connection = Improper( connection_members=site_list, improper_type=ctype, ) topology.add_connection(connection) return topology def _get_atoms(filename, topology, base_unyts, type_list): """Parse the atom information in the LAMMPS data file.""" with open(filename, "r") as lammps_file: for i, line in enumerate(lammps_file): if "atoms" in line.split(): n_atoms = int(line.split()[0]) if "Atoms" in line.split(): break atom_lines = open(filename, "r").readlines()[i + 2 : i + n_atoms + 2] for line in atom_lines: atom_line = line.split() atom_type = atom_line[2] charge = u.unyt_quantity(float(atom_line[3]), get_units(base_unyts, "charge")) coord = u.unyt_array( [float(atom_line[4]), float(atom_line[5]), float(atom_line[6])] ) * get_units(base_unyts, "length") site = Atom( charge=charge, position=coord, atom_type=copy.deepcopy(type_list[int(atom_type) - 1]), # 0-index molecule=Molecule( name=atom_line[1], number=int(atom_line[1]) - 1 ), # 0-index ) element = element_by_mass(site.atom_type.mass.value) site.name = element.name if element else site.atom_type.name site.element = element topology.add_site(site) return topology def _get_box_coordinates(filename, base_unyts, topology): """Parse box information.""" with open(filename, "r") as lammps_file: for line in lammps_file: if "xlo" in line.split(): break x_line = line.split() y_line = lammps_file.readline().split() z_line = lammps_file.readline().split() x = float(x_line[1]) - float(x_line[0]) y = float(y_line[1]) - float(y_line[0]) z = float(z_line[1]) - float(z_line[0]) # Check if box is triclinic tilts = lammps_file.readline().split() if "xy" in tilts: xy = float(tilts[0]) xz = float(tilts[1]) yz = float(tilts[2]) xhi = float(x_line[1]) - np.max([0.0, xy, xz, xy + xz]) xlo = float(x_line[0]) - np.min([0.0, xy, xz, xy + xz]) yhi = float(y_line[1]) - np.max([0.0, yz]) ylo = float(y_line[0]) - np.min([0.0, yz]) zhi = float(z_line[1]) zlo = float(z_line[0]) lx = xhi - xlo ly = yhi - ylo lz = zhi - zlo c = np.sqrt(lz**2 + xz**2 + yz**2) b = np.sqrt(ly**2 + xy**2) a = lx alpha = np.arccos((yz * ly + xy * xz) / (b * c)) beta = np.arccos(xz / c) gamma = np.arccos(xy / b) # Box Information lengths = u.unyt_array([a, b, c], get_units(base_unyts, "length")) angles = u.unyt_array([alpha, beta, gamma], get_units(base_unyts, "angle")) topology.box = Box(lengths, angles) else: # Box Information lengths = u.unyt_array([x, y, z], get_units(base_unyts, "length")) topology.box = Box(lengths) return topology def _get_ff_information(filename, base_unyts, topology): """Parse atom-type information.""" with open(filename, "r") as lammps_file: types = False for i, line in enumerate(lammps_file): if "atom" in line: n_atomtypes = int(line.split()[0]) types = True elif "Masses" in line: break if types is False: return topology mass_lines = open(filename, "r").readlines()[i + 2 : i + n_atomtypes + 2] type_list = list() for line in mass_lines: atom_type = AtomType( name=line.split()[0], mass=float(line.split()[1]) * get_units(base_unyts, "mass"), ) type_list.append(atom_type) with open(filename, "r") as lammps_file: for i, line in enumerate(lammps_file): if "Pair" in line: break # Need to figure out if we're going have mixing rules printed out # Currently only reading in LJ params warn_ljcutBool = False pair_lines = open(filename, "r").readlines()[i + 2 : i + n_atomtypes + 2] for i, pair in enumerate(pair_lines): if len(pair.split()) == 3: type_list[i].parameters["sigma"] = float(pair.split()[2]) * get_units( base_unyts, "length" ) type_list[i].parameters["epsilon"] = float(pair.split()[1]) * get_units( base_unyts, "energy" ) elif len(pair.split()) == 4: warn_ljcutBool = True if warn_ljcutBool: warnings.warn( "Currently not reading in LJ cutoff values." "These should be specified in the engine run files." ) return topology, type_list def _accepted_potentials(): """List of accepted potentials that LAMMPS can support.""" templates = PotentialTemplateLibrary() lennard_jones_potential = templates["LennardJonesPotential"] harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"] harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"] periodic_torsion_potential = templates["PeriodicTorsionPotential"] harmonic_improper_potential = templates["HarmonicImproperPotential"] opls_torsion_potential = templates["OPLSTorsionPotential"] accepted_potentialsList = [ lennard_jones_potential, harmonic_bond_potential, harmonic_angle_potential, periodic_torsion_potential, harmonic_improper_potential, opls_torsion_potential, ] return accepted_potentialsList def _validate_potential_compatibility(top): """Check compatability of topology object potentials with LAMMPSDATA format.""" pfilter = PotentialFilters.UNIQUE_EXPRESSION pot_types = check_compatibility( top, _accepted_potentials(), site_pfilter=pfilter, conn_pfilter=pfilter ) return pot_types def _validate_unit_compatibility(top, base_unyts): """Check compatability of topology object units with LAMMPSDATA format.""" for attribute in ["sites", "bonds", "angles", "dihedrals", "impropers"]: if attribute == "sites": atype = "atom_types" else: atype = attribute[:-1] + "_types" parametersList = [ (parameter, name) for attr_type in getattr(top, atype) for name, parameter in attr_type.parameters.items() ] for parameter, name in parametersList: assert np.isclose( float(base_unyts.convert_parameter(parameter, n_decimals=6, name=name)), parameter.value, atol=1e-3, ), f"Units System {base_unyts.usystem} is not compatible with {atype} with value {parameter}" def _write_header(out_file, top, atom_style, dihedral_parser): """Write Lammps file header.""" out_file.write( "{} written by {} at {} using the GMSO LAMMPS Writer\n\n\n".format( top.name if top.name is not None else "Topology", os.environ.get("USER"), str(datetime.datetime.now()), ) ) out_file.write("{:d} atoms\n".format(top.n_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)) if dihedral_parser in [parse_opls_style_dihedral]: # no layered dihedrals n_dihedrals = top.n_dihedrals elif dihedral_parser in [parse_charmm_style_dihedral]: # layered dihedrals n_dihedrals = 0 for dihedral in top.dihedrals: param = next(iter(dihedral.dihedral_type.parameters.values())) if isinstance(param, u.unyt_quantity): n_dihedrals += 1 else: n_dihedrals += len(param) elif dihedral_parser is None: n_dihedrals = 0 out_file.write("{:d} dihedrals\n".format(n_dihedrals)) 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)))) 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))) ) if top.n_angles > 0 and atom_style in ["full", "molecular"]: out_file.write( "{:d} angle types\n".format(len(top.angle_types(filter_by=pfilter))) ) if top.n_dihedrals > 0 and atom_style in ["full", "molecular"]: unique_dtypes = top.dihedral_types(filter_by=pfilter) nkeys = len(next(iter(unique_dtypes)).parameters.keys()) nparams = 0 # write out the total number of found for dihedrals for potential in unique_dtypes: for param in potential.parameters.values(): paramList = param.tolist() if isinstance(paramList, float): nparams += 1 else: for _ in param.tolist(): nparams += 1 ntypes = int( nparams / nkeys ) # allows us to count multiples for ones stored in a single object out_file.write("{:d} dihedral types\n".format(ntypes)) if top.n_impropers > 0 and atom_style in ["full", "molecular"]: out_file.write( "{:d} improper types\n".format(len(top.improper_types(filter_by=pfilter))) ) out_file.write("\n") def _write_box(out_file, top, base_unyts, cfactorsDict): """Write GMSO Topology box to LAMMPS file.""" if allclose_units( top.box.angles, u.unyt_array([90, 90, 90], "degree"), rtol=1e-5, atol=1e-8, ): box_lengths = [ float(base_unyts.convert_parameter(top.box.lengths[i], cfactorsDict)) for i in range(3) ] for i, dim in enumerate(["x", "y", "z"]): out_file.write( "{0:.6f} {1:.6f} {2}lo {2}hi\n".format(0, box_lengths[i], dim) ) out_file.write("0.000000 0.000000 0.000000 xy xz yz\n") else: box_lengths = [ float(base_unyts.convert_parameter(top.box.lengths[i], cfactorsDict)) for i in range(3) ] vectors = (box_lengths * top.box.get_unit_vectors().T).T xhi = vectors[0][0] yhi = vectors[1][1] zhi = vectors[2][2] xy = vectors[1][0] xz = vectors[2][0] yz = vectors[2][1] xlo = u.unyt_array(0, xy.units) ylo = u.unyt_array(0, xy.units) zlo = u.unyt_array(0, xy.units) xlo_bound = xlo + u.unyt_array(np.min([0.0, xy, xz, xy + xz]), xy.units) xhi_bound = xhi + u.unyt_array(np.max([0.0, xy, xz, xy + xz]), xy.units) ylo_bound = ylo + u.unyt_array(np.min([0.0, yz]), xy.units) yhi_bound = yhi + u.unyt_array(np.max([0.0, yz]), xy.units) zlo_bound = zlo zhi_bound = zhi out_file.write( "{0:.6f} {1:.6f} xlo xhi\n".format(xlo_bound.value, xhi_bound.value) ) out_file.write( "{0:.6f} {1:.6f} ylo yhi\n".format(ylo_bound.value, yhi_bound.value) ) out_file.write( "{0:.6f} {1:.6f} zlo zhi\n".format(zlo_bound.value, zhi_bound.value) ) out_file.write( "{0:.6f} {1:.6f} {2:.6f} xy xz yz\n".format(xy.value, xz.value, yz.value) ) def _write_atomtypes(out_file, top, base_unyts, cfactorsDict): """Write out atomtypes in GMSO topology to LAMMPS file.""" 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) for atom_type in atypesView: out_file.write( "{:d}\t{}\t# {}\n".format( atypesView.index(atom_type) + 1, base_unyts.convert_parameter(atom_type.mass, cfactorsDict), atom_type.name, ) ) 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 # Pair coefficients test_atomtype = top.sites[0].atom_type out_file.write(f"\nPair Coeffs # {test_atomtype.expression}\n") nb_style_orderTuple = ( "epsilon", "sigma", ) # this will vary with new pair styles param_labels = [ write_out_parameter_and_units( key, convert_kelvin_to_energy_units(test_atomtype.parameters[key], "kJ"), base_unyts, ) for key in nb_style_orderTuple ] out_file.write("#\t" + "\t".join(param_labels) + "\n") sorted_atomtypes = sorted(top.atom_types(filter_by=pfilter), key=lambda x: x.name) for idx, param in enumerate(sorted_atomtypes): out_file.write( "{}\t{:7}\t\t{:7}\t\t# {}\n".format( idx + 1, *[ base_unyts.convert_parameter( convert_kelvin_to_energy_units(param.parameters[key], "kJ"), cfactorsDict, n_decimals=5, ) for key in nb_style_orderTuple ], param.name, ) ) def _write_bondtypes(out_file, top, base_unyts, cfactorsDict): """Write out bonds to LAMMPS file.""" # TODO: Use any accepted lammps styles (only takes harmonic now) test_bondtype = top.bonds[0].bond_type out_file.write(f"\nBond Coeffs #{test_bondtype.name}\n") bond_style_orderTuple = ("k", "r_eq") param_labels = [ write_out_parameter_and_units( key, convert_kelvin_to_energy_units(test_bondtype.parameters[key], "kJ"), base_unyts, ) for key in bond_style_orderTuple ] out_file.write("#\t" + "\t".join(param_labels) + "\n") bond_types = list(top.bond_types(filter_by=pfilter)) bond_types.sort(key=lambda x: sorted(x.member_types)) for idx, bond_type in enumerate(bond_types): member_types = sorted([bond_type.member_types[0], bond_type.member_types[1]]) out_file.write( "{}\t{:7}\t{:7}\t\t# {}\t{}\n".format( idx + 1, *[ base_unyts.convert_parameter( convert_kelvin_to_energy_units(bond_type.parameters[key], "kJ"), cfactorsDict, n_decimals=6, ) for key in bond_style_orderTuple ], *member_types, ) ) return bond_types def _write_angletypes(out_file, top, base_unyts, cfactorsDict): """Write out angles to LAMMPS file.""" # TODO: Use any accepted lammps parameters, only harmonic now test_angletype = top.angles[0].angle_type out_file.write(f"\nAngle Coeffs #{test_angletype.name}\n") angle_style_orderTuple = ( "k", "theta_eq", ) # this will vary with new angle styles param_labels = [ write_out_parameter_and_units( key, convert_kelvin_to_energy_units(test_angletype.parameters[key], "kJ"), base_unyts, ) for key in angle_style_orderTuple ] out_file.write("#\t" + "\t".join(param_labels) + "\n") angle_types = list(top.angle_types(filter_by=pfilter)) angle_types.sort( key=lambda x: ( x.member_types[1], min(x.member_types[::2]), max(x.member_types[::2]), ) ) for idx, angle_type in enumerate(angle_types): out_file.write( "{}\t{:7}\t{:7}\t#{:11s}\t{:11s}\t{:11s}\n".format( idx + 1, *[ base_unyts.convert_parameter( convert_kelvin_to_energy_units( angle_type.parameters[key], "kJ" ), cfactorsDict, n_decimals=6, name=key, ) for key in angle_style_orderTuple ], *angle_type.member_types, ) ) return angle_types def _write_dihedraltypes(out_file, top, base_unyts, parser, cfactorsDict): """Write out dihedrals to LAMMPS file.""" test_dihedraltype = top.dihedrals[0].dihedral_type out_file.write(f"\nDihedral Coeffs #{test_dihedraltype.name}\n") param_labels0 = parser(test_dihedraltype) # tuple (paramsList, params_namesList) if isinstance( param_labels0[0][0], list ): # check for parsing out multiple instances from the dihedral param_labels = [ write_out_parameter_and_units( name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts ) for param, name in zip(param_labels0[0][0], param_labels0[1]) ] else: param_labels = [ write_out_parameter_and_units( name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts ) for param, name in zip(param_labels0[0], param_labels0[1]) ] out_file.write("#\t" + "\t".join(param_labels) + "\n") indexList = list(top.dihedral_types(filter_by=pfilter)) index_membersList = [ (dihedral_type, sort_by_types(dihedral_type)) for dihedral_type in indexList ] index_membersList.sort(key=lambda x: ([x[1][i] for i in [1, 2, 0, 3]])) # handle variable lengths for parameters base_msg = "{}\t" # handles index end_msg = "# {}\t{}\t{}\t{}\n" if ( parser.__name__ == "parse_opls_style_dihedral" ): # one opls parameter per dihedral type dihedral_typesList = [] for idx, (dihedral_type, members) in enumerate(index_membersList): param_labels = parser(dihedral_type) variable_msg = "{:8}\t" * len(param_labels[1]) full_msg = base_msg + variable_msg + end_msg out_file.write( full_msg.format( idx + 1, *[ base_unyts.convert_parameter( convert_kelvin_to_energy_units(parameter, "kJ"), cfactorsDict, n_decimals=6, name=parameterStr, ) for parameter, parameterStr in zip(*parser(dihedral_type)) ], *members, ) ) dihedral_typesList.append(dihedral_type) elif parser.__name__ == "parse_charmm_style_dihedral": ndecimalsDict = {"k": 6, "n": 0, "phi_eq": 0, "weights": 1} idx = 0 dihedral_typesList = [] for dihedral_type, members in index_membersList: parameter_termList, parameterStrList = parser(dihedral_type) variable_msg = "{:8}\t" * len(parameterStrList) full_msg = base_msg + variable_msg + end_msg for parameter_terms in parameter_termList: # list of params on each line out_file.write( full_msg.format( idx + 1, *[ base_unyts.convert_parameter( convert_kelvin_to_energy_units(parameter, "kJ"), cfactorsDict, n_decimals=ndecimalsDict[parameterStr], name=parameterStr, ) for parameter, parameterStr in zip( parameter_terms, parameterStrList ) ], *members, ) ) dihedral_typesList.append( dihedral_type ) # add dihedral type multiple times if it is layered idx += 1 return dihedral_typesList def parse_opls_style_dihedral(dihedral_type): """Take a dihedral type and list parameters as expected in lammps outputs.""" parametersList = [] namesList = ["k1", "k2", "k3", "k4"] for k in namesList: parametersList.append(dihedral_type.parameters[k]) return parametersList, namesList def parse_charmm_style_dihedral(dihedral_type, weightsArray=None): """Take a dihedral type and list parameters as expected in lammps outputs.""" kArray = dihedral_type.parameters["k"].flatten() nArray = dihedral_type.parameters["n"].flatten() phi_eqArray = dihedral_type.parameters["phi_eq"].flatten() if not weightsArray: # used for amber forcefield weights are 0 weightsArray = np.zeros(kArray.size) * u.dimensionless allParamsList = [] for a, b, c, d in zip(kArray, nArray, phi_eqArray, weightsArray): allParamsList.append([a, b, c, d]) return allParamsList, ["k", "n", "phi_eq", "weights"] def _write_impropertypes(out_file, top, base_unyts, parser, cfactorsDict): """Write out impropers to LAMMPS file.""" test_impropertype = top.impropers[0].improper_type out_file.write(f"\nImproper Coeffs #{test_impropertype.name}\n") param_labels0 = parser(test_impropertype) # tuple (paramsList, params_namesList) if isinstance( param_labels0[0][0], list ): # check for parsing out multiple instances from the dihedral param_labels = [ write_out_parameter_and_units( name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts ) for param, name in zip(param_labels0[0][0], param_labels0[1]) ] else: param_labels = [ write_out_parameter_and_units( name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts ) for param, name in zip(param_labels0[0], param_labels0[1]) ] out_file.write("#\t" + "\t".join(param_labels) + "\n") indexList = list(top.improper_types(filter_by=pfilter)) index_membersList = [ (improper_type, sort_by_types(improper_type)) for improper_type in indexList ] index_membersList.sort(key=lambda x: ([x[1][i] for i in [0, 1, 2, 3]])) # handle variable lengths for parameters base_msg = "{}\t" # handles index end_msg = "# {}\t{}\t{}\t{}\n" if ( parser.__name__ == "parse_cvff_style_improper" or "parse_harmonic_style_improper" ): # one cvff parameter per dihedral type ndecimalsDict = {"k": 6, "n": 0, "phi_eq": 0} idx = 0 for improper_type, members in index_membersList: param_labels = parser(improper_type) variable_msg = "{:8}\t" * len(param_labels[1]) full_msg = base_msg + variable_msg + end_msg out_file.write( full_msg.format( idx + 1, *[ base_unyts.convert_parameter( convert_kelvin_to_energy_units(parameter, "kJ"), cfactorsDict, n_decimals=ndecimalsDict[parameterStr], name=parameterStr, ) for parameter, parameterStr in zip(*parser(improper_type)) ], *members, ) ) idx += 1 return index_membersList # cvff is not layered, so no added to list def parse_cvff_style_improper(improper_type): """Take a dihedral type and list parameters as expected in lammps outputs.""" parametersList = [] namesList = ["k", "n", "phi_eq"] for k in namesList: parametersList.append(improper_type.parameters[k]) return parametersList, namesList def parse_harmonic_style_improper(improper_type): """Take a dihedral type and list parameters as expected in lammps outputs.""" parametersList = [] namesList = ["k", "phi_eq"] for k in namesList: parametersList.append(improper_type.parameters[k]) return parametersList, namesList 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") 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": atom_line = "{index:d}\t{type_index:d}\t{charge:.8}\t{x:.8}\t{y:.8}\t{z:.8}\n" elif atom_style == "molecular": atom_line = ( "{index:d}\t{moleculeid:d}\t{type_index:d}\t{x:.8}\t{y:.8}\t{z:.8}\n" ) 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" unique_sorted_typesList = sorted( top.atom_types(filter_by=pfilter), key=lambda x: x.name ) for i, site in enumerate(top.sites): 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, charge=base_unyts.convert_parameter( site.charge, cfactorsDict, n_decimals=6, ), x=base_unyts.convert_parameter( site.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 ), ) ) def _angle_order_sorter(angle_typesList): return [angle_typesList[i] for i in [1, 0, 2]] def _dihedral_order_sorter(dihedral_typesList): return [dihedral_typesList[i] for i in [1, 2, 0, 3]] def _improper_order_sorter(improper_typesList): return [improper_typesList[i] for i in [0, 1, 2, 3]] sorting_funcDict = { "bonds": None, "angles": _angle_order_sorter, "dihedrals": _dihedral_order_sorter, "impropers": _improper_order_sorter, } def _write_conn_data(out_file, top, connStr, sorted_typesList): """Write all connections to LAMMPS datafile.""" out_file.write(f"\n{connStr.capitalize()}\n\n") i = 0 for conn in getattr(top, connStr): ctype_members = sort_by_types(getattr(conn, connStr[:-1] + "_type")) indexList = [ ind for ind, ele in zip(count(), sorted_typesList) if sort_by_types(ele) == ctype_members ] for index in indexList: typeStr = f"{i+1:<6d}\t{index+1:<6d}\t" sorted_membersList = sort_connection_members(conn, sort_by="index", top=top) indexStr = "\t".join( [ str(top.get_index(member) + 1).ljust(6) for member in sorted_membersList ] ) out_file.write(typeStr + indexStr + "\n") i += 1 def _try_default_potential_conversions(top, potentialsDict): """Take a topology a convert all potentials to the style in potentialDict.""" for pot_container in potentialsDict: containerStr = pot_container[:-1] + "_types" if getattr(top, containerStr): for potential in potentialsDict[pot_container]: top.convert_potential_styles({pot_container: potential}) elif getattr(top, pot_container): raise AttributeError( f"Missing parameters in {pot_container} for {top.get_untyped(pot_container)}" ) def _default_lj_val(top, source): """Generate default lj non-dimensional values from topology.""" if source == "length": return copy.deepcopy( max(list(map(lambda x: x.parameters["sigma"], top.atom_types))) ) elif source == "energy": return copy.deepcopy( max(list(map(lambda x: x.parameters["epsilon"], top.atom_types))) ) elif source == "mass": return copy.deepcopy(max(list(map(lambda x: x.mass, top.atom_types)))) elif source == "charge": return copy.deepcopy(max(list(map(lambda x: x.charge, top.atom_types)))) else: raise ValueError( f"Provided {source} for default LJ cannot be found in the topology." ) def _identify_dihedral_parser(top, potential_typesDict): if not getattr(top, "dihedral_types"): return None # This is where dihedral_parser should get found parserDict = { "PeriodicTorsionPotential": parse_charmm_style_dihedral, "OPLSTorsionPotential": parse_opls_style_dihedral, } assert ( len(potential_typesDict["dihedral_types"]) == 1 ) # only allowing one potential type atm dihedralparser = parserDict[potential_typesDict["dihedral_types"].pop()] return dihedralparser def _identify_improper_parser(top, potential_typesDict): if not getattr(top, "improper_types"): return None # This is where improper_parser should be stored parserDict = { "PeriodicTorsionPotential": parse_cvff_style_improper, "HarmonicTorsionPotential": parse_harmonic_style_improper, "HarmonicImproperPotential": parse_harmonic_style_improper, } assert ( len(potential_typesDict["improper_types"]) == 1 ) # only allowing one potential type atm improper_parser = parserDict[potential_typesDict["improper_types"].pop()] return improper_parser