In this tutorial, we show how to generate a CG mapping matrix for a molecule given a bead mapping. The trajectory and topology file come from a AA simulation done in gromacs (see Molecules_CG_Mapping
folder). The protein is FF (diphenylalanine) and the solvent is a mixture of water and methanol.
# disable GPU
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
import hoomd, hoomd.htf as htf, hoomd.md
import matplotlib.pyplot as plt
import tensorflow as tf
import numpy as np
import MDAnalysis as mda
from os import path
import matplotlib.pyplot as plt,matplotlib
# Loading inputs
TPR = 'Molecules_CG_Mapping/nvt_prod.tpr'
tpr = mda.Universe(TPR)
TRAJECTORY = 'Molecules_CG_Mapping/traj.trr'
u = mda.Universe(TPR, TRAJECTORY)
# Generating Mapping Matrix for FF
protein_FF = u.select_atoms("resname PHE and resid 0:1")
Beads_distribution = [['N','H1','H2','H3'],
['CA','HA','CB','HB1','HB2'],
['CG','CD1','HD1','CD2','HD2','CE1','HE1','CE2','HE2','CZ','HZ'],
['C','O'],
['N','H'],
['CA','HA','CB','HB1','HB2'],
['CG','CD1','HD1','CD2','HD2','CE1','HE1','CE2','HE2','CZ','HZ'],
['C','O1','O2']]
mapping_FF = htf.matrix_mapping(protein_FF,Beads_distribution)
print (mapping_FF)
/home/mgholiza/.conda/envs/hoomd-tf2/lib/python3.7/site-packages/MDAnalysis/core/universe.py:171: UserWarning: No coordinate reader found for Molecules_CG_Mapping/nvt_prod.tpr. Skipping this file. 'this file.'.format(filename))
[[0.8224383 0.05918723 0.05918723 0.05918723 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0.44409524 0.03726984 0.44409524 0.03726984 0.03726984 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.15577257 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.42880501 0.57119499 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.93286579 0.06713421 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.44409524 0.03726984 0.44409524 0.03726984 0.03726984 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.15577257 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291 0. 0. 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.27291648 0.36354176 0.36354176]]
# Generating Mapping Matrix for Water
water = u.select_atoms("resname SOL and resid 500")
Beads_distribution = [['OW','HW1','HW2']]
mapping_water = htf.matrix_mapping(water,Beads_distribution)
print (mapping_water)
[[0.88809574 0.05595213 0.05595213]]
# Generating Mapping Matrix for Methanol
methanol = u.select_atoms("resname MET and resid 11665 ")
Beads_distribution_methanol = [['C','H','H','H','OA','HO']]
mapping_methanol = htf.matrix_mapping(methanol,Beads_distribution_methanol)
print (mapping_methanol)
[[0.37484707 0.03145832 0.03145832 0.03145832 0.49931966 0.03145832]]
# Getting the segment id of each molecule in topology
_,idx = np.unique(u.select_atoms('all').segids,return_index=True)
seg_id_list = u.select_atoms('all').segids[np.sort(idx)].tolist()
# Getting the list of every molecule type name in topology
_,idx = np.unique(u.atoms.resnames,return_index=True)
resname_list = u.atoms.resnames[np.sort(idx)].tolist()
# Getting list of atoms in each type of molecule
atoms_in_molecule_list = [protein_FF.names,
water.names,
methanol.names]
atoms_in_molecule_list = [protein_FF.names]
molecule_mapping_index = htf.find_molecules_from_topology(u,atoms_in_molecule_list, selection = "resname PHE")
molecule_mapping = mapping_FF
# get number of atoms
N = sum([len(m) for m in molecule_mapping_index])
# get number of molecules
M = len(molecule_mapping_index)
# get number of beads
B = molecule_mapping.shape[0]
#create a mass-weighted (M * bead_number) x N mapping operator
cg_mapping = htf.sparse_mapping([molecule_mapping for _ in molecule_mapping_index],
molecule_mapping_index)
assert cg_mapping.shape == (M * B, N)
print (M*B)
160