Example 9.1: Aqueous Polymer (PEO) MD in TIP3P Water using GROMACS

Example 9.1: Aqueous Polymer (PEO) MD in TIP3P Water using GROMACS

In this example, we demonstrate how to use fbtk to construct a system where a single polymer chain (PEO) is solvated by a large number of TIP3P water molecules. We then apply the latest OpenFF (Sage 2.2.1) force field along with the TIP3P water model and export the system to GROMACS format (.gro, .top) for simulation and analysis.

Water Model Handling and Simulation Stability

The combination of FBTK + OpenFF + GROMACS handles water models (specifically TIP3P in this case) in an extremely rational and robust manner.

1. Placement by FBTK: When packing a large number of solvent molecules around a polymer, the most critical factor is avoiding “atomic overlaps.” FBTK uses preliminary UFF-based geometries to create high-quality initial configurations without overlaps at high speed.
2. Blueprinting by OpenFF: When to_openff() is called, the OpenFF Toolkit integrates the specified force field files (e.g., ["openff-2.2.1.offxml", "tip3p.offxml"]). It assigns Sage 2.2.1 parameters to the PEO molecules and strict TIP3P water model parameters (O-H distance 0.9572 Å, H-O-H angle 104.52°, charges, etc.) as the physical “blueprint.”
3. Geometric Enforcement and Maintenance by GROMACS: The GROMACS .top file includes the [ settles ] directive for water molecules (SOL). Consequently, when performing Energy Minimization (EM) in GROMACS, the preliminary structure created by FBTK rapidly converges to the ideal TIP3P geometry. Furthermore, during the MD run, water molecules are maintained as strict rigid bodies by settles, ensuring physically correct behavior. Constraints (e.g., constraints = h-bonds) are applied to other hydrogen bonds, such as those in the PEO molecules, during MD.

Construction Script (run.py)

import fbtk
from openff.toolkit import ForceField # Import the ForceField class

# 1. Define the polymer (PEO, 20mer)
monomer = fbtk.Molecule.from_smiles("*COC*")
peo = fbtk.Molecule.from_polymer(monomer, degree=20, name="PEO")

# 2. Define water molecules (SOL)
water = fbtk.Molecule.from_smiles("O", name="SOL")

# 3. Build and relax the mixture
builder = fbtk.Builder(density=1.0)
builder.add_molecule(peo, count=1)
builder.add_molecule(water, count=1000)
system = builder.build()

# Structural relaxation (resolve overlaps before passing to OpenFF)
system.relax()

# 4. Convert to OpenFF format and save as GROMACS
# Apply both OpenFF 2.2.1 and the TIP3P water model
combined_forcefield = ForceField("openff-2.2.1.offxml", "tip3p.offxml")
interchange = system.to_openff(forcefield=combined_forcefield)
interchange.to_gromacs("peo_water_tip3p") # Updated output filename

GROMACS Configuration Files (.mdp)

Below are the configuration settings required to reproduce the simulation.

Energy Minimization (em.mdp)

integrator  = steep
emtol       = 1000.0
emstep      = 0.01
nsteps      = 5000

cutoff-scheme   = Verlet
nstlist         = 20
vdw-type        = Cut-off
vdw-modifier    = Force-switch
rvdw-switch     = 1.0
rvdw            = 1.2
coulombtype     = PME
rcoulomb        = 1.2

; Do not apply constraints during EM to allow the FBTK structure
; to naturally converge to the force field's ideal geometry.
constraints     = none
pbc             = xyz

NVT Equilibration MD (md.mdp)

integrator      = md
dt              = 0.002
nsteps          = 10000

nstxout-compressed = 1000
nstlog             = 1000
nstenergy          = 1000

cutoff-scheme   = Verlet
nstlist         = 20
vdw-type        = Cut-off
vdw-modifier    = Force-switch
rvdw-switch     = 1.0
rvdw            = 1.2
coulombtype     = PME
rcoulomb        = 1.2

tcoupl          = v-rescale
tc-grps         = System
tau-t           = 0.1
ref-t           = 300

; Constraints (PEO H-bonds; Water rigidity is handled by 'settles' in .top file)
constraints     = h-bonds
constraint-algorithm = LINCS
pbc             = xyz

Execution Commands

Run the simulation using the generated peo_water_tip3p.gro, peo_water_tip3p.top, and the MDP files above.

# 1. Energy Minimization
gmx grompp -f em.mdp -c peo_water_tip3p.gro -p peo_water_tip3p.top -o em.tpr
gmx mdrun -v -deffnm em

# 2. MD Execution
gmx grompp -f md.mdp -c em.gro -p peo_water_tip3p.top -o md.tpr -maxwarn 1
gmx mdrun -v -deffnm md

# 3. Convert to multi-frame gro (preparation for analysis)
echo 0 | gmx trjconv -f md.xtc -s md.tpr -o trajectory.gro

Application to Analysis (fbtk-analyze)

Using fbtk-analyze, intuitive analysis is possible without the need for index files.

Example: RDF between PEO oxygen and water oxygen

fbtk-analyze rdf -q "resname PEO and element O with resname SOL and element O" trajectory.gro -o peo_sol_rdf.dat

Results and Validation

• Fast Convergence: Energy minimization was completed in just 76 steps. Skipping constraints during EM allowed the FBTK packing structure to quickly snap into the force field’s ideal geometry.
• Thermodynamic Data: A stable simulation at 300K was maintained, and the residue names PEO and SOL were correctly managed (Performance: ~78.5 ns/day).

Summary

The fbtk + OpenFF + GROMACS toolchain significantly streamlines all stages of MD from aqueous system preparation using the latest force fields to analysis.