From PDB to Solvated Box#

Preparing, capping, and solvating experimental structures in a neutral periodic box.

Setting up a molecular system for molecular dynamics simulations involves multiple manual steps: adding missing hydrogens, choosing box geometry, determining clearance, calculating counterions for electroneutrality, and adding salt ions to reach physiological ionic strength.

In this recipe, we build an automated one-pass setup workflow using molsysmt.build.solvate() and verify the geometric and chemical integrity of the resulting box.

Added in version 1.0.0.

Loading Structure#

We load the crystalline coordinates of the Villin Headpiece (HP35):

import molsysmt as msm

# Load dry protein structure
molsys = msm.convert(msm.systems['chicken villin HP35']['chicken_villin_HP35.h5msm'], to_form='molsysmt.MolSys')

print(f"Dry protein contains {msm.get(molsys, n_atoms=True)} atoms in {msm.get(molsys, n_groups=True)} residues.")

Dry protein contains 605 atoms in 38 residues.

Terminal Capping#

We cap terminal residues to stabilize charge boundaries before solvating:

# Cap terminal residues
molsys = msm.build.add_missing_terminal_cappings(molsys)

print(f"Capped protein contains {msm.get(molsys, n_atoms=True)} atoms.")

Capped protein contains 605 atoms.

Solvating Box#

We construct a cubic simulation box with 1.0 nm clearance around the solute and add water and neutral counterions with OpenMM engine:

# Solvate in a cubic box with AMBER14 and TIP3P
solvated_molsys = msm.build.solvate(
    [molsys, {'forcefield': 'AMBER14', 'water_model': 'TIP3P'}],
    box_shape='cubic',
    clearance='10.0 angstroms',
    to_form='molsysmt.MolSys',
    engine='OpenMM'
)

msm.info(solvated_molsys)
form n_atoms n_groups n_components n_chains n_molecules n_entities n_waters n_ions n_peptides n_structures
molsysmt.MolSys 7231 2248 2211 2 2211 3 2208 2 1 1

Validating Geometry#

We verify box lengths, angles, volume, and count the added ions and water molecules:

# Extract periodic boundary box parameters
box = msm.get(solvated_molsys, element='system', box=True)
lengths, angles = msm.pbc.get_lengths_and_angles_from_box(box)
volume = msm.pbc.get_volume_from_box(box)

n_waters = msm.get(solvated_molsys, element='molecule', selection='molecule_type=="water"', n_molecules=True)
n_cl = msm.get(solvated_molsys, element='atom', selection='atom_name=="Cl-"', n_atoms=True)

print(f"Box lengths: {lengths[0]}")
print(f"Box volume: {volume[0]:.3f}")
print(f"Added solvent: {n_waters} waters, {n_cl} Cl- counterions.")

Box lengths: [4.268523 4.268523 4.268523] nanometer
Box volume: 77.774 nanometer ** 3
Added solvent: 2208 waters, 0 Cl- counterions.

Viewing Solvated System#

We visualize the solvated system in 3D using MolSysViewer:

msm.view(solvated_molsys)

See also