In [ ]:
import itertools
import string
import os
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from msibi import MSIBI, State, Pair, mie
import mdtraj as md
In [ ]:
t = md.load('traj_unwrapped.dcd', top='start_aa.hoomdxml')
In [ ]:
cg_idx = 0
start_idx = 0
n_propane = 1024 #passed later
propane_map = {0: [0, 1, 2]}
system_mapping = {}
for n in range(n_propane):
for bead, atoms in propane_map.items():
system_mapping[cg_idx] = [x + start_idx for x in atoms]
start_idx += len(atoms)
cg_idx += 1
# print(system_mapping)
With mapping for whole system, apply to all atom trajectory
In [ ]:
from mdtraj.core import element
list(t.top.atoms)[0].element = element.carbon
list(t.top.atoms)[0].element.mass
for atom in t.top.atoms:
atom.element = element.carbon
In [ ]:
cg_xyz = np.empty((t.n_frames, len(system_mapping), 3))
for cg_bead, aa_indices in system_mapping.items():
cg_xyz[:, cg_bead, :] = md.compute_center_of_mass(t.atom_slice(aa_indices))
# print(cg_xyz)
In [ ]:
cg_top = md.Topology()
for cg_bead in system_mapping.keys():
cg_top.add_atom('carbon', element.virtual_site, cg_top.add_residue('A', cg_top.add_chain()))
cg_traj = md.Trajectory(cg_xyz, cg_top, time=None, unitcell_lengths=t.unitcell_lengths, unitcell_angles=t.unitcell_angles)
cg_traj.save_dcd('cg_traj.dcd')
# print(cg_traj)
# print(cg_top)
# print(cg_xyz)
In [ ]:
pairs = cg_traj.top.select_pairs(selection1='name "carbon"', selection2='name "carbon"')
# mdtraj.compute_rdf(traj, pairs=None, r_range=None, bin_width=0.005, n_bins=None, periodic=True, opt=True)
r, g_r = md.compute_rdf(cg_traj, pairs=pairs, r_range=(0, 1.2), bin_width=0.005)
np.savetxt('rdfs_aa.txt', np.transpose([r, g_r]))
In [ ]:
fig, ax = plt.subplots()
ax.plot(r, g_r)
ax.set_xlabel("r")
ax.set_ylabel("g(r)")