%load_ext autoreload
%autoreload 2
import molsysmt as msm
Solvate#
molecular_system = msm.convert('pdb_id:1vii', to_form='molsysmt.MolSys')
molecular_system = msm.basic.remove(molecular_system, selection='atom_type=="H"')
molecular_system = msm.build.add_missing_terminal_cappings(molecular_system, N_terminal='ACE',
C_terminal='NME')
molecular_system = msm.build.add_missing_hydrogens(molecular_system, pH=7.4)
msm.info(molecular_system)
| form | n_atoms | n_groups | n_components | n_chains | n_molecules | n_entities | n_peptides | n_structures |
|---|---|---|---|---|---|---|---|---|
| molsysmt.MolSys | 605 | 38 | 1 | 1 | 1 | 1 | 1 | 1 |
#msm.physchem.charge([molecular_system, {'forcefield':'AMBER14'}], target='system')
msm.build.is_solvated(molecular_system)
False
molecular_system_cub = msm.build.solvate(molecular_system,
box_shape='cubic', clearance='14.0 angstroms',
to_form='molsysmt.MolSys', engine="OpenMM", verbose=False)
msm.build.is_solvated(molecular_system_cub)
True
msm.info(molecular_system_cub)
| form | n_atoms | n_groups | n_components | n_chains | n_molecules | n_entities | n_waters | n_ions | n_peptides | n_structures |
|---|---|---|---|---|---|---|---|---|---|---|
| molsysmt.MolSys | 11242 | 3585 | 3548 | 2 | 3548 | 3 | 3545 | 2 | 1 | 1 |
msm.info(molecular_system_cub, element='entity')
| index | name | type | n atoms | n groups | n components | n chains | n molecules |
|---|---|---|---|---|---|---|---|
| 0 | peptide 0 | peptide | 605 | 38 | 1 | 1 | 1 |
| 1 | water | water | 10635 | 3545 | 3545 | 1 | 3545 |
| 2 | CL | ion | 2 | 2 | 2 | 1 | 2 |
#msm.physchem.charge([molecular_system_cub, {'forcefield':'AMBER14', 'water_model':'TIP3P'}], element='system')
box, box_angles, box_shape = msm.get(molecular_system_cub, element='system', box=True, box_angles=True,
box_shape=True)
box
| Magnitude | [[[4.929236526299059 0.0 0.0] [0.0 4.929236526299059 0.0] [0.0 0.0 4.929236526299059]]] |
|---|---|
| Units | nanometer |
box_angles
| Magnitude | [[1.570796 1.570796 1.570796]] |
|---|---|
| Units | radian |
box_shape
'cubic'
molecular_system_cub = msm.pbc.wrap_to_pbc(molecular_system_cub, center_of_selection='molecule_type=="peptide"')
msm.view(molecular_system_cub, standard=True, with_water_as='surface')
# msm.view(molecular_system_cub, standardize=True, water_as_surface=True)
Adding ions#
PBC box geometry#
All periodic boxes used in molecular dynamics simulations (cubic, triclinic, hexagonal, dodecahedral or octahedral) are equivalent equivalent. All of them can be transformed into a triclinic box with the proper angles and edge lengths. See: Bekker, H. “Unification of Box Shapes in Molecular Simulations.” Journal of Computational Chemistry 18, no. 15 (1997): 1930–42. https://doi.org/10.1002/(sici)1096-987x(19971130)18:15<1930::aid-jcc8>3.0.co;2-p.
molecular_system_oct = msm.build.solvate(molecular_system, box_shape='truncated octahedral',
clearance='14.0 angstroms', engine='PDBFixer')
msm.info(molecular_system_oct)
| form | n_atoms | n_groups | n_components | n_chains | n_molecules | n_entities | n_waters | n_ions | n_peptides | n_structures |
|---|---|---|---|---|---|---|---|---|---|---|
| molsysmt.MolSys | 5056 | 1523 | 1486 | 2 | 1486 | 3 | 1483 | 2 | 1 | 1 |
molecular_system_oct = msm.pbc.wrap_to_pbc(molecular_system_oct, center_of_selection='molecule_type=="peptide"')
#msm.view(molecular_system_oct, standard=True, with_water_as=True)
In a triclinic box it is not sure that all elements in the unit cell can be considered first neighbors. Some pairs of atoms minimize their distance when one of them are located in a neighbor unit cell. But ¿Which one? Finding the periodic image that minimizes the distance is not in general as straight forward as it is if the box is cubic. This problem is known as “the minimum image convention”. Actually, the distance between any two atoms in a periodic box is not computed removing the PBC, or centering a unit cell in any of those atoms. It is solved finding the minimum image convention. Then, let’s see what happens when we take only the image of every atom with minimal distance to the center of the protein:
molecular_system_oct = msm.pbc.wrap_to_mic(molecular_system_oct, center_of_selection='molecule_type=="peptide"')
#msm.view(molecular_system_oct, standard=True, with_water_as=True)
The equivalent geometry is now recovered. It is then “a proof” of the equivalency between the triclinic box and the truncated octahedral box.
But why do we need a non cubic periodic box? In general a case, we want to be sure that a molecule is “solvated”. What does this mean? It means that our molecule is surrounded by a thick enough layer of water molecules. ¿This can be accomplished by a cubic periodic box? Yes of course. But it can also be achieved with other geometries making use of a lower number of water molecules. Which means that running a molecular simulation with these other geometries will be computationally cheaper than with a periodic cube:
molecular_system_cub = msm.build.solvate(molecular_system, box_shape='cubic', clearance='14.0 angstroms',
engine='PDBFixer')
molecular_system_oct = msm.build.solvate(molecular_system, box_shape='truncated octahedral', clearance='14.0 angstroms',
engine='PDBFixer')
molecular_system_dod = msm.build.solvate(molecular_system, box_shape='rhombic dodecahedral', clearance='14.0 angstroms',
engine='PDBFixer')
n_waters_cub = msm.get(molecular_system_cub, element='system', n_waters=True)
n_waters_oct = msm.get(molecular_system_oct, element='system', n_waters=True)
n_waters_dod = msm.get(molecular_system_dod, element='system', n_waters=True)
n_waters_oct_to_cub = round(100.0* n_waters_oct/n_waters_cub, 2)
n_waters_dod_to_cub = round(100.0* n_waters_dod/n_waters_cub, 2)
print('Cubic box: {} water (100.0% -cubic reference-)'.format(n_waters_cub))
print('Truncated octahedral box: {} water ({}% -cubic reference-)'.format(n_waters_oct, n_waters_oct_to_cub))
print('Rhombic dodecahedron box: {} water ({}% -cubic reference-)'.format(n_waters_dod, n_waters_dod_to_cub))
Cubic box: 2035 water (100.0% -cubic reference-)
Truncated octahedral box: 1483 water (72.87% -cubic reference-)
Rhombic dodecahedron box: 1333 water (65.5% -cubic reference-)
volume_cub = msm.get(molecular_system_cub, element='system', box_volume=True)
volume_oct = msm.get(molecular_system_oct, element='system', box_volume=True)
volume_dod = msm.get(molecular_system_dod, element='system', box_volume=True)
volume_oct_to_cub = round(100.0* volume_oct[0]/volume_cub[0], 2).magnitude
volume_dod_to_cub = round(100.0* volume_dod[0]/volume_cub[0], 2).magnitude
print('Cubic box: {} volue (100% -cubic reference-)'.format(volume_cub[0]))
print('Truncated octahedral box: {} volume ({}% -cubic reference-)'.format(volume_oct[0], volume_oct_to_cub))
print('Rhombic dodecahedron box: {} volume ({}% -cubic reference-)'.format(volume_dod[0], volume_dod_to_cub))
Cubic box: 73.36535832264707 nanometer ** 3 volue (100% -cubic reference-)
Truncated octahedral box: 56.47667916903153 nanometer ** 3 volume (76.98% -cubic reference-)
Rhombic dodecahedron box: 51.87714237412465 nanometer ** 3 volume (70.71% -cubic reference-)