"""HOOMD simulation format."""
import itertools
import operator
import warnings
from collections import namedtuple
import numpy as np
import parmed as pmd
import mbuild as mb
from mbuild.utils.conversion import RB_to_OPLS
from mbuild.utils.io import import_
from mbuild.utils.sorting import natural_sort
from .hoomd_snapshot import to_hoomdsnapshot
gsd = import_("gsd")
gsd.hoomd = import_("gsd.hoomd")
hoomd = import_("hoomd")
hoomd.md = import_("hoomd.md")
[docs]def create_hoomd_simulation(
structure,
r_cut,
ref_distance=1.0,
ref_mass=1.0,
ref_energy=1.0,
auto_scale=False,
snapshot_kwargs={},
pppm_kwargs={"Nx": 8, "Ny": 8, "Nz": 8, "order": 4},
init_snap=None,
restart=None,
nlist=hoomd.md.nlist.cell,
):
"""Convert a parametrized pmd.Structure to hoomd.SimulationContext.
Parameters
----------
structure : parmed.Structure
ParmEd Structure object
r_cut : float
Cutoff radius in simulation units
ref_distance : float, optional, default=1.0
Reference distance for unit conversion (from Angstrom)
ref_mass : float, optional, default=1.0
Reference mass for unit conversion (from Dalton)
ref_energy : float, optional, default=1.0
Reference energy for unit conversion (from kcal/mol)
auto_scale : bool, optional, default=False
Scale to reduced units by automatically using the largest sigma value
as ref_distance, largest mass value as ref_mass, and largest epsilon
value as ref_energy
snapshot_kwargs : dict
Kwargs to pass to to_hoomdsnapshot
pppm_kwargs : dict
Kwargs to pass to hoomd's pppm function
init_snap : hoomd.data.SnapshotParticleData, optional, default=None
Initial snapshot to which to add the ParmEd structure object
(useful for rigid bodies)
restart : str, optional, default=None
Path to the gsd file from which to restart the simulation.
https://hoomd-blue.readthedocs.io/en/v2.9.4/restartable-jobs.html
Note: It is assumed that the ParmEd structure and the system in
restart.gsd contain the same types. The ParmEd structure is still used
to initialize the forces, but restart.gsd is used to initialize the
system state (e.g., particle positions, momenta, etc).
nlist : hoomd.md.nlist, default=hoomd.md.nlist.cell
Type of neighborlist to use, see
https://hoomd-blue.readthedocs.io/en/stable/module-md-nlist.html
for more information.
Returns
-------
hoomd_objects : list
List of hoomd objects created during conversion
ReferenceValues : namedtuple
Values used in scaling
Note
----
While the hoomd objects are returned, the hoomd.SimulationContext is
accessible via `hoomd.context.current`. If you pass a non-parametrized
pmd.Structure, you will not have angle, dihedral, or force field
information. You may be better off creating a hoomd.Snapshot.
Reference units should be expected to convert parmed Structure units:
--- angstroms, kcal/mol, and daltons
"""
if isinstance(structure, mb.Compound):
raise ValueError(
"You passed mb.Compound to create_hoomd_simulation, there will be "
"no angles, dihedrals, or force field parameters. Please use "
"hoomd_snapshot.to_hoomdsnapshot to create a hoomd.Snapshot, then "
"create your own hoomd context and pass your hoomd.Snapshot to "
"hoomd.init.read_snapshot()"
)
elif not isinstance(structure, pmd.Structure):
raise ValueError(
"Please pass a parmed.Structure to create_hoomd_simulation"
)
_check_hoomd_version()
version_numbers = _check_hoomd_version()
if float(version_numbers[0]) >= 3:
raise RuntimeError("This method does not support HOOMD-blue v3.x.")
hoomd_objects = [] # Potential adaptation for Hoomd v3 API
if auto_scale:
if not all([i == 1 for i in (ref_distance, ref_energy, ref_mass)]):
warnings.warn(
"Autoscale option selected--provided reference values will not "
"be used."
)
pair_coeffs = list(
set((a.type, a.epsilon, a.sigma) for a in structure.atoms)
)
ref_mass = max([atom.mass for atom in structure.atoms])
ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1]
ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2]
ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"])
ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy)
if not hoomd.context.current:
hoomd.context.initialize("")
if restart is None:
snapshot, _ = to_hoomdsnapshot(
structure,
ref_distance=ref_distance,
ref_mass=ref_mass,
ref_energy=ref_energy,
**snapshot_kwargs,
hoomd_snapshot=init_snap,
)
hoomd_objects.append(snapshot)
hoomd_system = hoomd.init.read_snapshot(snapshot)
hoomd_objects.append(hoomd_system)
else:
with gsd.hoomd.open(restart) as f:
snapshot = f[-1]
hoomd_objects.append(snapshot)
hoomd_system = hoomd.init.read_gsd(restart, restart=restart)
hoomd_objects.append(hoomd_system)
print("Simulation initialized from restart file")
nl = nlist()
nl.reset_exclusions(exclusions=["1-2", "1-3"])
hoomd_objects.append(nl)
if structure.atoms[0].type != "":
print("Processing LJ and QQ")
lj = _init_hoomd_lj(
structure,
nl,
r_cut,
ref_distance=ref_distance,
ref_energy=ref_energy,
)
qq = _init_hoomd_qq(structure, nl, r_cut, **pppm_kwargs)
hoomd_objects.append(lj)
hoomd_objects.append(qq)
if structure.adjusts:
print("Processing 1-4 interactions, adjusting neighborlist exclusions")
lj_14, qq_14 = _init_hoomd_14_pairs(
structure,
nl,
r_cut,
ref_distance=ref_distance,
ref_energy=ref_energy,
)
hoomd_objects.append(lj_14)
hoomd_objects.append(qq_14)
if structure.bond_types:
print("Processing harmonic bonds")
harmonic_bond = _init_hoomd_bonds(
structure, ref_distance=ref_distance, ref_energy=ref_energy
)
hoomd_objects.append(harmonic_bond)
if structure.angle_types:
print("Processing harmonic angles")
harmonic_angle = _init_hoomd_angles(structure, ref_energy=ref_energy)
hoomd_objects.append(harmonic_angle)
if structure.dihedral_types:
print("Processing periodic torsions")
periodic_torsions = _init_hoomd_dihedrals(
structure, ref_energy=ref_energy
)
hoomd_objects.append(periodic_torsions)
if structure.rb_torsion_types:
print("Processing RB torsions")
rb_torsions = _init_hoomd_rb_torsions(structure, ref_energy=ref_energy)
hoomd_objects.append(rb_torsions)
print("HOOMD SimulationContext updated from ParmEd Structure")
return hoomd_objects, ref_values
def _init_hoomd_lj(structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0):
"""LJ parameters."""
# Identify the unique atom types before setting
atom_type_params = {}
for atom in structure.atoms:
if atom.type not in atom_type_params:
atom_type_params[atom.type] = atom.atom_type
# Set the hoomd parameters for self-interactions
lj = hoomd.md.pair.lj(r_cut, nl)
for name, atom_type in atom_type_params.items():
lj.pair_coeff.set(
name,
name,
sigma=atom_type.sigma / ref_distance,
epsilon=atom_type.epsilon / ref_energy,
)
if atom_type.epsilon / ref_energy == 0:
lj.pair_coeff.set(name, name, r_cut=0)
# Cross interactions, mixing rules, NBfixes
all_atomtypes = sorted(atom_type_params.keys())
for a1, a2 in itertools.combinations_with_replacement(all_atomtypes, 2):
nb_fix_info = atom_type_params[a1].nbfix.get(a2, None)
# nb_fix_info = (rmin, eps, rmin14, eps14)
if nb_fix_info is None:
# No nbfix means use mixing rule to find cross-interaction
if structure.combining_rule == "lorentz":
sigma = (
atom_type_params[a1].sigma + atom_type_params[a2].sigma
) / (2 * ref_distance)
epsilon = (
(
atom_type_params[a1].epsilon
* atom_type_params[a2].epsilon
)
/ ref_energy**2
) ** 0.5
elif structure.combining_rule == "geometric":
sigma = (
(atom_type_params[a1].sigma * atom_type_params[a2].sigma)
/ ref_distance**2
) ** 0.5
epsilon = (
(
atom_type_params[a1].epsilon
* atom_type_params[a2].epsilon
)
/ ref_energy**2
) ** 0.5
else:
raise ValueError(
f"Mixing rule {structure.combining_rule} not supported, "
"use lorentz"
)
else:
# If we have nbfix info, use it
sigma = nb_fix_info[0] / (ref_distance * (2 ** (1 / 6)))
epsilon = nb_fix_info[1] / ref_energy
lj.pair_coeff.set(a1, a2, sigma=sigma, epsilon=epsilon)
if epsilon == 0:
lj.pair_coeff.set(a1, a2, r_cut=0)
return lj
def _init_hoomd_qq(structure, nl, r_cut, Nx=1, Ny=1, Nz=1, order=4):
"""Charge interactions."""
charged = hoomd.group.charged()
if len(charged) == 0:
print("No charged groups found, ignoring electrostatics")
return None
else:
qq = hoomd.md.charge.pppm(charged, nl)
qq.set_params(Nx, Ny, Nz, order, r_cut)
return qq
def _init_hoomd_14_pairs(
structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0
):
"""Special_pairs to handle 14 scalings.
See discussion: https://groups.google.com/forum/#!topic/hoomd-users/
iZ9WCpHczg0
"""
# Update neighborlist to exclude 1-4 interactions,
# but impose a special_pair force to handle these pairs
nl.exclusions.append("1-4")
if hoomd.context.current.system_definition.getPairData().getN() == 0:
print("No 1,4 pairs found in hoomd snapshot")
return None, None
lj_14 = hoomd.md.special_pair.lj()
qq_14 = hoomd.md.special_pair.coulomb()
params_14 = {}
# Identify unique 14 scalings
for adjust in structure.adjusts:
t1 = adjust.atom1.type
t2 = adjust.atom2.type
ps = "-".join(sorted([t1, t2]))
if ps not in params_14:
params_14[ps] = adjust.type
for name, adjust_type in params_14.items():
lj_14.pair_coeff.set(
name,
sigma=adjust_type.sigma / ref_distance,
# The adjust epsilon already carries the scaling
epsilon=adjust_type.epsilon / ref_energy,
# Do NOT use hoomd's alpha to modify any LJ terms
alpha=1,
r_cut=r_cut,
)
qq_14.pair_coeff.set(name, alpha=adjust_type.chgscale, r_cut=r_cut)
return lj_14, qq_14
def _init_hoomd_bonds(structure, ref_distance=1.0, ref_energy=1.0):
"""Harmonic bonds."""
# Identify the unique bond types before setting
bond_type_params = {}
for bond in structure.bonds:
t1, t2 = bond.atom1.type, bond.atom2.type
t1, t2 = sorted([t1, t2], key=natural_sort)
if t1 != "" and t2 != "":
bond_type = "-".join((t1, t2))
if bond_type not in bond_type_params:
bond_type_params[bond_type] = bond.type
# Set the hoomd parameters
harmonic_bond = hoomd.md.bond.harmonic()
for name, bond_type in bond_type_params.items():
# A (paramerized) parmed structure with no bondtype
# is because of constraints
if bond_type is None:
print("Bond with no bondtype detected, setting coefficients to 0")
harmonic_bond.bond_coeff.set(name, k=0, r0=0)
else:
harmonic_bond.bond_coeff.set(
name,
k=2 * bond_type.k * ref_distance**2 / ref_energy,
r0=bond_type.req / ref_distance,
)
return harmonic_bond
def _init_hoomd_angles(structure, ref_energy=1.0):
"""Harmonic angles."""
# Identify the unique angle types before setting
angle_type_params = {}
for angle in structure.angles:
t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
t1, t3 = sorted([t1, t3], key=natural_sort)
angle_type = "-".join((t1, t2, t3))
if angle_type not in angle_type_params:
angle_type_params[angle_type] = angle.type
# set the hoomd parameters
harmonic_angle = hoomd.md.angle.harmonic()
for name, angle_type in angle_type_params.items():
harmonic_angle.angle_coeff.set(
name,
t0=np.deg2rad(angle_type.theteq),
k=2 * angle_type.k / ref_energy,
)
return harmonic_angle
def _init_hoomd_dihedrals(structure, ref_energy=1.0):
"""Periodic dihedrals (dubbed harmonic dihedrals in HOOMD)."""
# Identify the unique dihedral types before setting
# need Hoomd 2.8.0 to use proper dihedral implemtnation
# from this PR https://github.com/glotzerlab/hoomd-blue/pull/492
version_numbers = _check_hoomd_version()
if float(version_numbers[0]) < 2 or float(version_numbers[1]) < 8:
from mbuild.exceptions import MBuildError
raise MBuildError("Please upgrade Hoomd to at least 2.8.0")
dihedral_type_params = {}
for dihedral in structure.dihedrals:
t1, t2 = dihedral.atom1.type, dihedral.atom2.type
t3, t4 = dihedral.atom3.type, dihedral.atom4.type
if [t2, t3] == sorted([t2, t3], key=natural_sort):
dihedral_type = "-".join((t1, t2, t3, t4))
else:
dihedral_type = "-".join((t4, t3, t2, t1))
if dihedral_type not in dihedral_type_params:
if isinstance(dihedral.type, pmd.DihedralType):
dihedral_type_params[dihedral_type] = dihedral.type
elif isinstance(dihedral.type, pmd.DihedralTypeList):
if len(dihedral.type) > 1:
warnings.warn(
"Multiple dihedral types detected for single dihedral, "
"will ignore all except first dihedral type. First "
"dihedral type: {}".format(dihedral.type[0])
)
dihedral_type_params[dihedral_type] = dihedral.type[0]
# Set the hoomd parameters
# These are periodic torsions
periodic_torsion = hoomd.md.dihedral.harmonic()
for name, dihedral_type in dihedral_type_params.items():
periodic_torsion.dihedral_coeff.set(
name,
k=2 * dihedral_type.phi_k / ref_energy,
d=1,
n=dihedral_type.per,
phi_0=np.deg2rad(dihedral_type.phase),
)
return periodic_torsion
def _init_hoomd_rb_torsions(structure, ref_energy=1.0):
"""RB dihedrals (implemented as OPLS dihedrals in HOOMD)."""
# Identify the unique dihedral types before setting
dihedral_type_params = {}
for dihedral in structure.rb_torsions:
t1, t2 = dihedral.atom1.type, dihedral.atom2.type
t3, t4 = dihedral.atom3.type, dihedral.atom4.type
if [t2, t3] == sorted([t2, t3], key=natural_sort):
dihedral_type = "-".join((t1, t2, t3, t4))
else:
dihedral_type = "-".join((t4, t3, t2, t1))
if dihedral_type not in dihedral_type_params:
dihedral_type_params[dihedral_type] = dihedral.type
# Set the hoomd parameter
rb_torsion = hoomd.md.dihedral.opls()
for name, dihedral_type in dihedral_type_params.items():
F_coeffs = RB_to_OPLS(
dihedral_type.c0 / ref_energy,
dihedral_type.c1 / ref_energy,
dihedral_type.c2 / ref_energy,
dihedral_type.c3 / ref_energy,
dihedral_type.c4 / ref_energy,
dihedral_type.c5 / ref_energy,
error_if_outside_tolerance=False,
)
rb_torsion.dihedral_coeff.set(
name, k1=F_coeffs[1], k2=F_coeffs[2], k3=F_coeffs[3], k4=F_coeffs[4]
)
return rb_torsion
def _check_hoomd_version():
version = hoomd.__version__
version_numbers = version.split(".")
return version_numbers