Tutorial on the usage of intrinsic_density_micelle tool

Here, two examples on how to use the intrinsic_density_micelle tool for the calculation of the distribution of components of a micelle. First we will show how to use intrinsic_density to calculate the density of the hydrophobic block of the polymer. Then we will show how to use intrinsic_density_water to calculate the spherical density of the water in the simulation

spherical_density tool

Before starting any analysis, load the neccesary modules for this class.

from  utils_mda import MDA_input
#from pysoftk.pol_analysis.tools.utils_mda import MDA_input
from utils_tools import *
#from pysoftk.pol_analysis.tools.utils_tools import *
from clustering import SCP
#from pysoftk.pol_analysis.clustering import SCP
from make_micelle_whole import micelle_whole
#from pysoftk.pol_analysis.make_micelle_whole import micelle_whole
from intrinsic_density_micelle import intrinsic_density
#from pysoftk.pol_analysis.intrinsic_density_micelle import intrinsic_density

import numpy as np
import pandas as pd
  1. Select your simulation files, it is recommended to use a tpr file for the topology and xtc file for the trajectory.

topology='data/intrinsic_density_movie.tpr'
trajectory='data/intrinsic_density_movie.xtc'
  1. Run the clustering

results_name='data/pictures_tutorial/cyclic_intrinsic_scp_result'

cluster_cutoff = 12

start=0
stop=10001
step=1

atom_names=['C02T']

c = SCP(topology, trajectory).spatial_clustering_run(start, stop, step, atom_names, cluster_cutoff, results_name)


#load the results from SCP
resids_total='data/pictures_tutorial/cyclic_intrinsic_scp_result.parquet'
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [00:00<00:00, 22.20it/s]
The file data/pictures_tutorial/cyclic_intrinsic_scp_result.parquet has been successfully created.
Function spatial_clustering_run Took 2.2583 seconds
  1. Obtain the largest micelle from the clustering pandas dataframe

largest_micelle_resids = micelle_whole(topology, trajectory).obtain_largest_micelle_resids(resids_total)
  1. ‘largest_micelle_resids’ is a np.array with the resids of the molecules that belong to the same cluster of the steps of the trajectory where SCP was ran.

largest_micelle_resids
[array([ 1,  2,  3,  4,  6,  7,  8, 11, 12, 13, 14, 15, 17, 18, 19]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 11, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20]),
 array([ 2,  4,  6,  7, 13, 14, 15, 17, 18, 19, 20])]

Now, let’s obtain the coordinates of the largest micelle made whole across the periodic boundary conditions (PBC).

  1. Let’s define the resname of the molecules that we want to reconstruct with the make_whole function. More than one resname can be inputted. Note that it should be the resname of the molecules of the largest_micelle array.

resname=['LIG']
  1. Also, the user needs to define the start, step and final step of frames in which the user wants to run the analysis on. Note that they need to be the same as the ones you ran the SCP clustering on.

start=0
stop=10000
step=1
  1. Obtain the whole coordinates of the micelle. This is key to perform the densitity calculation taking properly into account the PBC

atom_pos = micelle_whole(topology, trajectory).running_make_cluster_whole(resname, largest_micelle_resids, start, stop, step)
  0%|                                                                                                                                                   | 0/11 [00:00<?, ?it/s]/home/raquellrdc/Desktop/PhD/pysoftk/alejandro_newest_releast_check/pysoftk_analysis_code/test_final/make_micelle_whole.py:347: FutureWarning: arrays to stack must be passed as a "sequence" type such as list or tuple. Support for non-sequence iterables such as generators is deprecated as of NumPy 1.16 and will raise an error in the future.
  atom_positions_over_trajectory = list(tqdm(map(self.make_cluster_whole, frames, resname, cluster_resids_f[0],
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [00:16<00:00,  1.47s/it]
Elapsed time for matrix calculation: 17.7199 seconds

Now, atom_pos contains the coordinates of all the atoms of the micelle made whole at each time step selected. In each array, the first element is the time frame of the analysis, and the second the positions array

  1. We are almost ready for the density calculation. But first, we need to define the atoms for which you want to calculate the density. In this case, we are going to calculate the density of the hydrophobic block of the polymer, so we are going to select the main backbone atoms of the hydrophobic block.

MA_names = ['C028', 'C024', 'C022', 'C02E', 'C02I', 'C02M', 'C02Q', 'C02U', 'C02Y', 'C00S', 'C00Q', 'C00M',
             'C00L', 'C027',  'C023', 'C021',  'C02H',  'C02L',  'C02P',  'C02T', 'C02X', 'C00U',  'C00R',  'C00P',
             'O00A', 'O00D', 'O00C', 'O01G', 'O01F', 'O01D', 'O01A', 'O019', 'O016', 'O015', 'O010', 'O012',
             'O00B', 'O00E', 'O01I', 'O01H', 'O01C', 'O01B', 'O01E', 'O018', 'O017', 'O014', 'O00Z', 'O013',
             'C00O', 'C00W', 'C031', 'C030', 'C02W', 'C02S', 'C02O', 'C02K', 'C02G', 'C02D', 'C026', 'C02B']
  1. Also, we need to define the heavy atoms of the solvent. This is neccesary even if we are calculating the density of the polymer components. This means that we need a trajectory with the solvent for the intrinsic density calculation.

water_atom=['OW'] #for water the heavy atom is the Oxygen
  1. That’s it! We have all the neccesary input files for the spherical density calculation of the hydrophobic component. This function outputs two arrays, the first one are the density values and the second one are the binned space values. This function outputs these two arrays for easier plotting.

density_whole_MA, binned_space = intrinsic_density(topology, trajectory).run_intrinsic_density(largest_micelle_resids, atom_pos, MA_names, MA_names, water_atom, start, stop, step)
  0%|                                                                                                                                                   | 0/11 [00:00<?, ?it/s]/home/raquellrdc/Desktop/PhD/polymer/pysoftwhere/pySoftWhere/pySoftWhere/ICSI/icsi.py:85: RuntimeWarning: Mean of empty slice
  interface_vals[i][j]=np.nanmean([interface_vals[divmod(ip,len(interface_vals))[1]][divmod(jp,len(interface_vals))[1]] for ip,jp in product(n_i,n_j)])
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [00:34<00:00,  3.12s/it]
/home/raquellrdc/Desktop/PhD/pysoftk/alejandro_newest_releast_check/pysoftk_analysis_code/test_final/intrinsic_density_micelle.py:280: RuntimeWarning: divide by zero encountered in divide
  final_density=intrinsic_density_profile/S_bar
/home/raquellrdc/Desktop/PhD/pysoftk/alejandro_newest_releast_check/pysoftk_analysis_code/test_final/intrinsic_density_micelle.py:280: RuntimeWarning: invalid value encountered in divide
  final_density=intrinsic_density_profile/S_bar

Note that MA_names is specified twice. This is because for the intrinsic density you need to define the atoms that belong to the (hydrophobic) core of the micelle. So the first specification is used to describe the atoms in the core, and the second time to define the atoms you want to calculate the density of. You do not need ot include all atoms in the core, you can only input at least one atom of the backbone chain per hydrophobic monomer

Also, note that we have left the values to be binned as the default one, but you may want to change them depending on your system

Now, let’s plot the density to see how it looks.

import matplotlib.pyplot as plt

plt.plot(binned_space, density_whole_MA)

plt.xlabel('r($\AA$)')
plt.ylabel('$\widetilde{\u03C1}(r)$($\AA^{-3}$)')
plt.ylim(0, 0.5) #the peak at the interface can be quite high, so this limit is to properly observe the distribution of the component within the micelle
plt.title('Intrinsic density of MA')
Text(0.5, 1.0, 'Intrinsic density of MA')
../_images/intrinsic_density_tutorial_30_1.png

Remember that the more data (frames and atoms) that you input, the smoother your density will be

intrinsic_density_water tool

Now, let’s calculate the density of the water in the same simulation. For this we need a different tool, intrinsic_density_water. We are going to use the same trajectory, so no need to load it again. We just need to start by importing the intrinsic_density_water_micelle class

#from pysoftk.pol_analysis.intrinsic_density_water_micelle  import intrinsic_density_water
from intrinsic_density_water_micelle import intrinsic_density_water

We have already obtained the resids of the polymer belonging to the largest micelle as well as the atom positions of the micelle made whole across the PBC. So the only inputs needed for this function are the name of the solvant atoms that we want to calculate the density of. In this case the solvent is water, so we select the oxygen atoms.

water_atoms=['OW']

We are ready to run the density of the water! - This will take longer, since there are many water atoms in the system

density_whole_water, binned_space_water = intrinsic_density_water(topology, trajectory).run_intrinsic_density(largest_micelle_resids, atom_pos, water_atoms, MA_names, water_atoms, start, stop, step)
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [02:58<00:00, 16.22s/it]
/home/raquellrdc/Desktop/PhD/pysoftk/alejandro_newest_releast_check/pysoftk_analysis_code/test_final/intrinsic_density_water_micelle.py:332: RuntimeWarning: divide by zero encountered in divide
  final_density=intrinsic_density_profile/S_bar
/home/raquellrdc/Desktop/PhD/pysoftk/alejandro_newest_releast_check/pysoftk_analysis_code/test_final/intrinsic_density_water_micelle.py:332: RuntimeWarning: invalid value encountered in divide
  final_density=intrinsic_density_profile/S_bar

At this point, it is also needed to specify the core of the micelle, since the intrinsic density is estimated with respect to the core-shell interface. That is why we need to again input ‘MA_names’ that contains some of the atoms of the hydrophobic block of the polymer.

Now let’s plot it to see how it looks

plt.plot(binned_space_water, density_whole_water)

plt.xlabel('r($\AA$)')
plt.ylabel('$\widetilde{\u03C1}(r)$($\AA^{-3}$)')
plt.ylim(0, 0.1)
plt.xlim(-40, 75) # the densities of the far end reach very high values that are just artifacts of the calculations of the intrinsic density performed very far away from the micelle.
plt.title('Intrinsic density of Water')
Text(0.5, 1.0, 'Intrinsic density of Water')
../_images/intrinsic_density_tutorial_41_1.png

Looks like this micelle is not very stable since it has a lot of water in the core!