In this tutorial, we show how to generate a CG mapping matrix for a molecule given a bead distribution.
import hoomd, hoomd.htf as htf, hoomd.md
import numpy as np
import MDAnalysis as mda
# 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)
[[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]]