%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 | 11374 | 3629 | 3592 | 2 | 3592 | 3 | 3589 | 2 | 1 | 1 |
msm.info(molecular_system_cub, element='entity')
| index | name | type | n atoms | n groups | n components | n chains | n molecules |
|---|---|---|---|---|---|---|---|
| 0 | VILLIN | peptide | 605 | 38 | 1 | 1 | 1 |
| 1 | water | water | 10767 | 3589 | 3589 | 1 | 3589 |
| 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.943507344884426 0.0 0.0] [0.0 4.943507344884426 0.0] [0.0 0.0 4.943507344884426]]] |
|---|---|
| Units | nanometer |
box_angles
| Magnitude | [[89.99998127603166 89.99998127603166 89.99998127603166]] |
|---|---|
| Units | degree |
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 | 5068 | 1527 | 1490 | 2 | 1490 | 3 | 1487 | 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: 2037 water (100.0% -cubic reference-)
Truncated octahedral box: 1487 water (73.0% -cubic reference-)
Rhombic dodecahedron box: 1332 water (65.39% -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.40742661931107 nanometer ** 3 volue (100% -cubic reference-)
Truncated octahedral box: 56.5090633589026 nanometer ** 3 volume (76.98% -cubic reference-)
Rhombic dodecahedron box: 51.906889151968734 nanometer ** 3 volume (70.71% -cubic reference-)