Mapping Water Molecules¶
In this notebook we will cover how to map molecules in different ways and look at some of the things we can do with them.
In this demonstration, we will load data from a GROMACS simulation and therefore, we need to define a set of units and a file reader object to use. For this reason, we have changed the imports a little bit to keep the code to minimum.
[1]:
import mdsuite as mds
import mdsuite.file_io.chemfiles_read
from mdsuite.utils import Units
from zinchub import DataHub
import shutil
import h5py as hf
import numpy as np
2022-07-27 15:46:18.687881: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:18.687908: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.
2022-07-27 15:46:22.721545: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcublas.so.11'; dlerror: libcublas.so.11: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.721903: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcublasLt.so.11'; dlerror: libcublasLt.so.11: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.722203: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcufft.so.10'; dlerror: libcufft.so.10: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.722501: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcurand.so.10'; dlerror: libcurand.so.10: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.722797: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcusolver.so.11'; dlerror: libcusolver.so.11: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.723129: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcusparse.so.11'; dlerror: libcusparse.so.11: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.723432: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudnn.so.8'; dlerror: libcudnn.so.8: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /opt/slurm/lib:/opt/slurm/lib:
2022-07-27 15:46:22.723466: W tensorflow/core/common_runtime/gpu/gpu_device.cc:1850] Cannot dlopen some GPU libraries. Please make sure the missing libraries mentioned above are installed properly if you would like to use GPU. Follow the guide at https://www.tensorflow.org/install/gpu for how to download and setup the required libraries for your platform.
Skipping registering GPU devices...
Loading the data¶
In this tutorial we are using 50 ns simulations of 14 water molecules in a continuum fluid performed with GROMACS. We will use pure atomistic naming as well as ligand naming, the topology files for which are contained on DataHub.
[2]:
water = DataHub(url="https://github.com/zincware/DataHub/tree/main/Water_14_Gromacs", tag="v0.1.0")
water.get_file('./')
file_paths = [
f for f in water.file_raw
]
Project definition¶
Here we create the project and define some custom units used by GROMACS.
[3]:
project = mds.Project("Mapping_Molecules")
gmx_units = Units(
time=1e-12,
length=1e-10,
energy=1.6022e-19,
NkTV2p=1.6021765e6,
boltzmann=8.617343e-5,
temperature=1,
pressure=100000,
)
2022-07-27 15:46:25,821 - INFO: Creating new project Mapping_Molecules
Mapping molecules with SMILES¶
In this section we take a look at how one can map molecules using SMILES strings.
[4]:
traj_path = file_paths[2]
topol_path = file_paths[0]
file_reader = mdsuite.file_io.chemfiles_read.ChemfilesRead(
traj_file_path=traj_path, topol_file_path=topol_path
)
water_chemical = project.add_experiment(
name="water_chemical",
timestep=0.002,
temperature=300.0,
units=gmx_units,
simulation_data=file_reader,
)
water_chemical.sample_rate=5000
water_chemical.run.CoordinateUnwrapper()
2022-07-27 15:46:26,276 - INFO: Creating a new experiment!
100%|███████████████████████████████████| 1/1 [00:00<00:00, 1.91it/s]
2022-07-27 15:46:27.439156: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations: AVX2 AVX512F FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
Applying transformation 'Unwrapped_Positions' to 'O': 100%|█| 1/1 [00:
Applying transformation 'Unwrapped_Positions' to 'H': 100%|█| 1/1 [00:
The first thing we need to do is define the molecule that will be mapped using the in-built MDSuite molecule data-class
[5]:
chemical_water = mds.Molecule(
name='water',
smiles="[H]O[H]",
cutoff=1.7,
amount=14,
mol_pbc=True
)
[6]:
water_chemical.run.MolecularMap(
molecules=[chemical_water]
)
2022-07-27 15:46:28,786 - INFO: Building molecular graph from configuration for water
100%|████████████████████████████████| 42/42 [00:00<00:00, 430.84it/s]
2022-07-27 15:46:29,062 - INFO: Performing molecule number isomorphism test.
2022-07-27 15:46:29,063 - INFO: Amount of molecules test passed.
2022-07-27 15:46:29,064 - INFO: Performing group equality isomorphism test.
2022-07-27 15:46:29,065 - INFO: Group equality isomorphism test passed.
2022-07-27 15:46:29,186 - INFO: Mapping molecule graphs onto trajectory for water
100%|███████████████████████████████████| 1/1 [00:00<00:00, 4.74it/s]
Applying transformation 'Unwrapped_Positions' to 'water': 100%|█| 1/1
Mapping Molecules with a reference dict¶
If you do not have particles with chemical names but you nonetheless wish to construct groups out of particles, this can be achieved by using a reference dict.
In this example, we use the ligand naming from GROMACS to construct water molecules.
[7]:
traj_path = file_paths[2]
topol_path = file_paths[1]
file_reader = mdsuite.file_io.chemfiles_read.ChemfilesRead(
traj_file_path=traj_path, topol_file_path=topol_path
)
water_ligand = project.add_experiment(
name="water_ligand",
timestep=0.002,
temperature=300.0,
units=gmx_units,
simulation_data=file_reader,
)
water_ligand.run.CoordinateUnwrapper()
2022-07-27 15:46:31,453 - INFO: Creating a new experiment!
2022-07-27 15:46:32,256 - WARNING: WARNING element OW has been assigned mass=0.0
100%|███████████████████████████████████| 1/1 [00:00<00:00, 1.73it/s]
Applying transformation 'Unwrapped_Positions' to 'OW': 100%|█| 1/1 [00
Applying transformation 'Unwrapped_Positions' to 'HW1': 100%|█| 1/1 [0
Applying transformation 'Unwrapped_Positions' to 'HW2': 100%|█| 1/1 [0
Keep in mind, as the particles are not named from the periodic tables, important properties such as mass will need to be filled in manually.
[8]:
water_ligand.species['OW'].mass = [15.999]
water_ligand.species['HW1'].mass = [1.00784]
water_ligand.species['HW2'].mass = [1.00784]
In this case, the molecule will be defined a little bit differently.
[9]:
ligand_water = mds.Molecule(
name='water',
cutoff=1.7,
amount=14,
species_dict={"HW1": 1, "OW": 1, "HW2": 1},
mol_pbc=True
)
[10]:
water_ligand.run.MolecularMap(
molecules=[ligand_water]
)
2022-07-27 15:46:35,910 - INFO: Building molecular graph from configuration for water
100%|████████████████████████████████| 42/42 [00:00<00:00, 541.50it/s]
2022-07-27 15:46:36,163 - INFO: Performing molecule number isomorphism test.
2022-07-27 15:46:36,165 - INFO: Amount of molecules test passed.
2022-07-27 15:46:36,166 - INFO: Performing group equality isomorphism test.
2022-07-27 15:46:36,167 - INFO: Group equality isomorphism test passed.
2022-07-27 15:46:36,289 - INFO: Mapping molecule graphs onto trajectory for water
100%|███████████████████████████████████| 1/1 [00:00<00:00, 5.89it/s]
Applying transformation 'Unwrapped_Positions' to 'water': 100%|█| 1/1
What information is stored?¶
So the molecule mapping itself was quick and easy, but what information has been stored along the way?
All meta-data about the molecules is stored in the experiment class under molecules. Let’s take a look at what this contains.
[11]:
water_chemical.molecules.keys()
[11]:
dict_keys(['water'])
This dict will contain all of the molecules that have been mapped, but this is not the information about the molecules, for that, we need to look at the water molecule.
[12]:
water_chemical.molecules['water']
[12]:
MoleculeInfo(name='water', n_particles=14, properties=[], mass=18.015, charge=0, groups={'0': {'H': [0, 1], 'O': [0]}, '1': {'H': [2, 3], 'O': [1]}, '2': {'H': [4, 5], 'O': [2]}, '3': {'H': [6, 7], 'O': [3]}, '4': {'H': [8, 9], 'O': [4]}, '5': {'H': [10, 11], 'O': [5]}, '6': {'H': [12, 13], 'O': [6]}, '7': {'H': [14, 15], 'O': [7]}, '8': {'H': [16, 17], 'O': [8]}, '9': {'H': [18, 19], 'O': [9]}, '10': {'H': [20, 21], 'O': [10]}, '11': {'H': [22, 23], 'O': [11]}, '12': {'H': [24, 25], 'O': [12]}, '13': {'H': [26, 27], 'O': [13]}})
Three of these are fairly trivial and we can look at them quickly, groups will require some more attention.
[13]:
print(f"n_particles: {water_chemical.molecules['water'].n_particles}")
print(f"mass: {water_chemical.molecules['water'].mass}")
n_particles: 14
mass: 18.015
Now let’s take a look at the groups key.
[14]:
print(water_chemical.molecules['water'].groups)
{'0': {'H': [0, 1], 'O': [0]}, '1': {'H': [2, 3], 'O': [1]}, '2': {'H': [4, 5], 'O': [2]}, '3': {'H': [6, 7], 'O': [3]}, '4': {'H': [8, 9], 'O': [4]}, '5': {'H': [10, 11], 'O': [5]}, '6': {'H': [12, 13], 'O': [6]}, '7': {'H': [14, 15], 'O': [7]}, '8': {'H': [16, 17], 'O': [8]}, '9': {'H': [18, 19], 'O': [9]}, '10': {'H': [20, 21], 'O': [10]}, '11': {'H': [22, 23], 'O': [11]}, '12': {'H': [24, 25], 'O': [12]}, '13': {'H': [26, 27], 'O': [13]}}
The groups key contains direct information about which atoms belong to which molecule, for example, the 10th molecule of water (id=9) consists of Hydrogen atoms 18 and 19 and oxygen atom 10.
[15]:
print(water_chemical.molecules['water'].groups['9'])
{'H': [18, 19], 'O': [9]}
With this information you can compute values, e.g. diffusion coefficients with only the atoms belonging to a single molecule using the atom_select arguments in the calculator.
Analysis with molecules¶
Now that we have seen how we can build molecules and what information this gives is, let’s look at what we can analyse using them.
Angular Distribution Functions (ADFs)¶
First things first, let’s confirm we are working with water by looking at the angular distribution function of the atoms.
[16]:
water_chemical.run.AngularDistributionFunction(
number_of_configurations=5000, number_of_bins=500, norm_power=8
)
0%| | 0/1 [00:00<?, ?it/s]2022-07-27 15:46:54.814097: W tensorflow/core/framework/cpu_allocator_impl.cc:82] Allocation of 4877315008 exceeds 10% of free system memory.
100%|███████████████████████████████████| 1/1 [00:35<00:00, 35.14s/it]
[16]:
Exp1_Angular_Distribution_Function_1
Looking at the O_H_H ADF in he top right we see a strong max peak at 109.591 degrees corresponding well with the bond angle of an SPCE model (109.47) as was used in the simulation. It is also worth noting that the oxygen triplet angle looks similar to that measured in QM and experimental studies.
When we want to study the molecular ADF we have two choices, we can either pass it as a species argument to the calculator if only one is desired, we we can call the calculator with the molecules=True keyword as we will do here.
[17]:
water_chemical.run.AngularDistributionFunction(
molecules=True, number_of_configurations=3000, number_of_bins=500, norm_power=8
)
100%|███████████████████████████████████| 1/1 [00:00<00:00, 1.02it/s]
[17]:
Exp1_Angular_Distribution_Function_2
In this case we have increased the norm power to suppress the noise floor and highlight only the most dominant peaks.
In the water molecule ADF it we do not see any clear stacking or structure suggesting there is not special organization of the molecules in these simulations.
Radial Distribution Functions (RDFs)¶
Now let’s look at the radial structure and distribution of particles in space of both the atomistic system and the molecules. This is where molecule mapping can be very helpful as often we are more interested in the positions of the molecules themselves and not necessarily those of the atoms.
[18]:
water_chemical.run.RadialDistributionFunction(
number_of_configurations=4000, start=100, stop=5100, number_of_bins=500
)
Running mini batch loop 1 / 7: 100%|████| 1/1 [00:01<00:00, 1.94s/it]
Running mini batch loop 2 / 7: 100%|████| 1/1 [00:00<00:00, 16.57it/s]
Running mini batch loop 3 / 7: 100%|████| 1/1 [00:00<00:00, 13.83it/s]
Running mini batch loop 4 / 7: 100%|████| 1/1 [00:00<00:00, 3.61it/s]
Running mini batch loop 5 / 7: 100%|████| 1/1 [00:00<00:00, 14.30it/s]
Running mini batch loop 6 / 7: 100%|████| 1/1 [00:00<00:00, 14.62it/s]
Running mini batch loop 7 / 7: 100%|████| 1/1 [00:00<00:00, 17.32it/s]
[18]:
Exp1_Radial_Distribution_Function_3
In the case of the hydrogen-hydrogen and the oxygen-hydrogen we can see clear peaks where the bond distance is fixed. Using the cursor to hover over the points in the plot we can identify a bond distance between hydrogens of approximately 0.163 nm, in good agreement with experimental values. The oxygen-hydrogen bond sits around 0.09 nm, also in good agreement with experiment values.
[19]:
water_ligand.run.RadialDistributionFunction(
number_of_configurations=5200, start=0, stop=5100, number_of_bins=500, molecules=True
)
Running mini batch loop 1 / 1: 100%|████| 1/1 [00:00<00:00, 3.33it/s]
[19]:
Exp2_Radial_Distribution_Function_4
Diffusion Coefficients¶
Now let’s start looking at the diffusion coefficients of the atoms and molecules.
[20]:
water_chemical.run.EinsteinDiffusionCoefficients(data_range=500)
2022-07-27 15:47:19,040 - INFO: starting Einstein Diffusion Computation
2022-07-27 15:47:19,042 - INFO: starting Einstein Diffusion Computation
O: 100%|████████████████████████████████| 1/1 [00:10<00:00, 10.47s/it]
H: 100%|████████████████████████████████| 1/1 [00:13<00:00, 13.40s/it]
[20]:
Exp1_Einstein Self-Diffusion Coefficients_5
[21]:
water_chemical.run.EinsteinDiffusionCoefficients(molecules=True, data_range=500)
2022-07-27 15:47:44,695 - INFO: starting Einstein Diffusion Computation
2022-07-27 15:47:44,697 - INFO: starting Einstein Diffusion Computation
water: 100%|████████████████████████████| 1/1 [00:11<00:00, 11.98s/it]
[21]:
Exp1_Einstein Self-Diffusion Coefficients_6
Group-wise analysis¶
Say we want to study a specific molecule. We only want the diffusion coefficients, ADFs, and RDFs of the atoms in that one molecule. This can be achieved with the MDSuite atom-selection command and is included here as a demonstration of the flexibility of the software.
First things first, let’s select a molecule group to study, say the first water molecule.
[22]:
water_group = water_chemical.molecules['water'].groups['0']
print(water_group)
{'H': [0, 1], 'O': [0]}
[23]:
water_chemical.run.RadialDistributionFunction(atom_selection={'H': [0, 1], 'O': [0]}, number_of_configurations=2517)
Running mini batch loop 1 / 5: 100%|████| 1/1 [00:00<00:00, 26.45it/s]
Running mini batch loop 2 / 5: 100%|████| 1/1 [00:00<00:00, 26.31it/s]
Running mini batch loop 3 / 5: 100%|████| 1/1 [00:00<00:00, 22.78it/s]
Running mini batch loop 4 / 5: 100%|████| 1/1 [00:00<00:00, 24.23it/s]
Running mini batch loop 5 / 5: 100%|████| 1/1 [00:00<00:00, 24.50it/s]
[23]:
Exp1_Radial_Distribution_Function_7
[24]:
water_chemical.run.AngularDistributionFunction(atom_selection=water_group, number_of_configurations=100)
100%|███████████████████████████████████| 1/1 [00:00<00:00, 10.64it/s]
[24]:
Exp1_Angular_Distribution_Function_8
[25]:
water_chemical.run.EinsteinDiffusionCoefficients(atom_selection=water_group, data_range=2500)
2022-07-27 15:47:59,776 - INFO: starting Einstein Diffusion Computation
2022-07-27 15:47:59,778 - INFO: starting Einstein Diffusion Computation
O: 100%|████████████████████████████████| 1/1 [00:07<00:00, 7.16s/it]
H: 100%|████████████████████████████████| 1/1 [00:07<00:00, 7.16s/it]
[25]:
Exp1_Einstein Self-Diffusion Coefficients_9
BMIM-BF4 RTIL¶
[26]:
bmim_bf4 = DataHub(
url="https://github.com/zincware/DataHub/tree/main/Bmim_BF4", tag="v0.1.0"
)
bmim_bf4.get_file()
bmim_file = bmim_bf4.file_raw
[27]:
project.add_experiment("bmim_bf4", simulation_data=bmim_file)
project.experiments.bmim_bf4.run.CoordinateUnwrapper()
2022-07-27 15:48:18,189 - INFO: Creating a new experiment!
100%|███████████████████████████████████| 1/1 [00:10<00:00, 10.56s/it]
Applying transformation 'Unwrapped_Positions' to 'N': 100%|█| 1/1 [00:
Applying transformation 'Unwrapped_Positions' to 'C': 100%|█| 1/1 [00:
Applying transformation 'Unwrapped_Positions' to 'H': 100%|█| 1/1 [00:
Applying transformation 'Unwrapped_Positions' to 'B': 100%|█| 1/1 [00:
Applying transformation 'Unwrapped_Positions' to 'F': 100%|█| 1/1 [00:
[28]:
bmim_molecule = mdsuite.Molecule(
name="bmim",
species_dict={"C": 8, "N": 2, "H": 15},
amount=50,
cutoff=1.9,
reference_configuration_idx=100,
)
bf_molecule = mdsuite.Molecule(
name="bf4",
smiles="[B-](F)(F)(F)F",
amount=50,
cutoff=2.4,
reference_configuration_idx=100,
)
project.experiments["bmim_bf4"].run.MolecularMap(
molecules=[bmim_molecule, bf_molecule]
)
2022-07-27 15:48:32,230 - INFO: Building molecular graph from configuration for bmim
100%|█████████████████████████████| 1250/1250 [00:18<00:00, 66.84it/s]
2022-07-27 15:49:27,661 - INFO: Performing molecule number isomorphism test.
2022-07-27 15:49:27,661 - INFO: Amount of molecules test passed.
2022-07-27 15:49:27,662 - INFO: Performing group equality isomorphism test.
2022-07-27 15:49:27,663 - INFO: Group equality isomorphism test passed.
2022-07-27 15:49:27,818 - INFO: Mapping molecule graphs onto trajectory for bmim
24it [00:11, 2.11it/s]
Applying transformation 'Positions' to 'bmim': 100%|█| 1/1 [00:00<00:0
2022-07-27 15:49:39,732 - INFO: Building molecular graph from configuration for bf4
100%|███████████████████████████████| 250/250 [00:02<00:00, 98.29it/s]
2022-07-27 15:49:45,479 - INFO: Performing molecule number isomorphism test.
2022-07-27 15:49:45,481 - INFO: Amount of molecules test passed.
2022-07-27 15:49:45,482 - INFO: Performing group equality isomorphism test.
2022-07-27 15:49:45,483 - INFO: Group equality isomorphism test passed.
2022-07-27 15:49:45,624 - INFO: Mapping molecule graphs onto trajectory for bf4
100%|███████████████████████████████████| 1/1 [00:00<00:00, 4.87it/s]
Applying transformation 'Positions' to 'bf4': 100%|█| 1/1 [00:00<00:00
[29]:
project.experiments.bmim_bf4.run.RadialDistributionFunction(
number_of_configurations=300, number_of_bins=100
)
100%|███████████████████████████████| 300/300 [03:10<00:00, 1.57it/s]
[29]:
Exp3_Radial_Distribution_Function_10
[30]:
project.experiments.bmim_bf4.run.RadialDistributionFunction(
number_of_configurations=500, number_of_bins=100, molecules=True
)
Running mini batch loop 1 / 5: 100%|████| 1/1 [00:00<00:00, 17.72it/s]
Running mini batch loop 2 / 5: 100%|████| 1/1 [00:00<00:00, 17.54it/s]
Running mini batch loop 3 / 5: 100%|████| 1/1 [00:00<00:00, 16.84it/s]
Running mini batch loop 4 / 5: 100%|████| 1/1 [00:00<00:00, 15.33it/s]
Running mini batch loop 5 / 5: 100%|████| 1/1 [00:00<00:00, 16.69it/s]
[30]:
Exp3_Radial_Distribution_Function_11