Force Field Tutorial¶

Learn how to work with force fields in MolPy! Force fields define the interaction parameters for molecular simulations - atom types, bond styles, and more.

What is a Force Field?¶

A ForceField in MolPy represents the collection of interaction parameters:

• Atom Types: Define atom properties (mass, charge, etc.)
• Bond Styles: Define how bonds interact
• Angle Styles: Define angle interactions
• Dihedral Styles: Define dihedral interactions
• Pair Styles: Define non-bonded interactions

Let's explore:

In [1]:

from molpy.core.forcefield import AtomType, BondStyle, ForceField

from molpy.core.forcefield import AtomType, BondStyle, ForceField

Creating a Force Field¶

You can create a force field from scratch or load from files:

In [2]:

# Create an empty force field
ff = ForceField(name="my_forcefield", units="real")
print(f"Force field: {ff}")

# Create an empty force field ff = ForceField(name="my_forcefield", units="real") print(f"Force field: {ff}")

Force field: <ForceField: my_forcefield>

Defining Atom Types¶

Atom types define the basic properties of atoms:

In [3]:

from molpy.core.forcefield import AtomStyle

# Define an atom style
atom_style = ff.def_style(AtomStyle, "atomic")

# Add atom types to the style
carbon_type = AtomType("C", mass=12.01, charge=0.0)
hydrogen_type = AtomType("H", mass=1.008, charge=0.0)

atom_style.types.add(carbon_type)
atom_style.types.add(hydrogen_type)

print(f"Atom types: {ff.get_types(AtomType)}")

from molpy.core.forcefield import AtomStyle # Define an atom style atom_style = ff.def_style(AtomStyle, "atomic") # Add atom types to the style carbon_type = AtomType("C", mass=12.01, charge=0.0) hydrogen_type = AtomType("H", mass=1.008, charge=0.0) atom_style.types.add(carbon_type) atom_style.types.add(hydrogen_type) print(f"Atom types: {ff.get_types(AtomType)}")

Atom types: [<AtomType: H>, <AtomType: C>]

Defining Bond Styles¶

Bond styles define how bonded atoms interact:

In [4]:

from molpy.core.forcefield import BondType

# Define a bond style (e.g., harmonic)
bond_style = ff.def_style(BondStyle, "harmonic")

# Add bond type parameters
ch_bond = BondType("C-H", carbon_type, hydrogen_type, k=340.0, r0=1.09)
bond_style.types.add(ch_bond)

print(f"Bond styles: {ff.get_styles(BondStyle)}")

from molpy.core.forcefield import BondType # Define a bond style (e.g., harmonic) bond_style = ff.def_style(BondStyle, "harmonic") # Add bond type parameters ch_bond = BondType("C-H", carbon_type, hydrogen_type, k=340.0, r0=1.09) bond_style.types.add(ch_bond) print(f"Bond styles: {ff.get_styles(BondStyle)}")

Bond styles: [<BondStyle: harmonic>]

Using Typifiers¶

Typifiers automatically assign atom types to structures based on SMARTS patterns:

In [5]:

# Example: Using a typifier to assign types
# First, create a simple molecule
from molpy.core.atomistic import Atomistic

# Create a simple methane molecule
methane = Atomistic()
c = methane.def_atom(symbol="C", xyz=[0.0, 0.0, 0.0])
h1 = methane.def_atom(symbol="H", xyz=[1.0, 0.0, 0.0])
h2 = methane.def_atom(symbol="H", xyz=[-1.0, 0.0, 0.0])
h3 = methane.def_atom(symbol="H", xyz=[0.0, 1.0, 0.0])
h4 = methane.def_atom(symbol="H", xyz=[0.0, 0.0, 1.0])

methane.def_bond(c, h1)
methane.def_bond(c, h2)
methane.def_bond(c, h3)
methane.def_bond(c, h4)

print(f"Created molecule: {methane}")
print(f"Atoms: {len(methane.atoms)}, Bonds: {len(methane.bonds)}")

# Note: To use OplsAtomisticTypifier, you need a force field with SMARTS patterns
# from molpy.typifier.atomistic import OplsAtomisticTypifier
# typifier = OplsAtomisticTypifier(ff, skip_atom_typing=False)
# typifier.typify(methane)  # modifies methane in place
# print("Types assigned to atoms")

# Example: Using a typifier to assign types # First, create a simple molecule from molpy.core.atomistic import Atomistic # Create a simple methane molecule methane = Atomistic() c = methane.def_atom(symbol="C", xyz=[0.0, 0.0, 0.0]) h1 = methane.def_atom(symbol="H", xyz=[1.0, 0.0, 0.0]) h2 = methane.def_atom(symbol="H", xyz=[-1.0, 0.0, 0.0]) h3 = methane.def_atom(symbol="H", xyz=[0.0, 1.0, 0.0]) h4 = methane.def_atom(symbol="H", xyz=[0.0, 0.0, 1.0]) methane.def_bond(c, h1) methane.def_bond(c, h2) methane.def_bond(c, h3) methane.def_bond(c, h4) print(f"Created molecule: {methane}") print(f"Atoms: {len(methane.atoms)}, Bonds: {len(methane.bonds)}") # Note: To use OplsAtomisticTypifier, you need a force field with SMARTS patterns # from molpy.typifier.atomistic import OplsAtomisticTypifier # typifier = OplsAtomisticTypifier(ff, skip_atom_typing=False) # typifier.typify(methane) # modifies methane in place # print("Types assigned to atoms")

Created molecule: <Atomistic, 5 atoms (C:1 H:4), 4 bonds, with coords>
Atoms: 5, Bonds: 4

Merging Force Fields¶

You can combine multiple force fields:

In [6]:

# Create two force fields
ff1 = ForceField(name="ff1")
ff2 = ForceField(name="ff2")

# Merge ff2 into ff1
ff1.merge(ff2)

print(f"Merged force field: {ff1}")

# Create two force fields ff1 = ForceField(name="ff1") ff2 = ForceField(name="ff2") # Merge ff2 into ff1 ff1.merge(ff2) print(f"Merged force field: {ff1}")

Merged force field: <ForceField: ff1>