Use Aqua to clusterize and calculate lifetimes of Hbond networks in pure-water simulation

lifetime.py

from aqua import *
import matplotlib.pyplot as plt
import operator
# import misc

# Importing the topology of the system from a pdb
system=msystem('solvate.pdb',verbose=True)
# Encoding the info about donors and acceptors for the whole system. -needed
# to compute hbonds-
acc_don_all=system.selection_hbonds('all')
# Opening the dcd file
system.load_traj('solvate2.dcd',verbose=True)

# Creating an empty array to store the microstates
mss_exp=numpy.empty((system.traj[0].total_frames,
  system.num_waters,17),dtype=int,order='F')
# loop running over the total number of frames in the dcd
for ii in range(system.traj[0].total_frames):
    # computing hbonds for the entire system in frame ii
    hbout=system.hbonds(definition='R(o,o)-Ang(o,o,h)',
      acc_don_A=acc_don_all,frame=ii,infile=True,optimize=True)
    # Creating microstates
    mss,mss_ind=system.mss_hbonds_wat_prot(hbonds=hbout,verbose=True)
    # Storing microstates
    mss_exp[ii,:,:]=mss[:,:]

# Initialization of class kinetic_analysis with the trajectory of water
# microstates
kinlab=kinetic_analysis(mss_exp)
# Computing the kinetic network 
kinlab.kinetic_network(verbose=True)

# Run MCL
kinlab.network.symmetrize(new=False,verbose=True)
kinlab.network.mcl()
# Access clusters: kinlab.network.cluster...
# or run gradient clusters: kinlab.network.gradient_clusters()
num_nodes = kinlab.network.num_nodes
aux_list = numpy.empty(num_nodes,dtype=int,order='F')
for ii in range(num_nodes):
  aux_list[ii] = kinlab.network.node[ii].cluster
new_num_frames = kinlab.traj_nodes.shape[0]
kinlab.traj_clusters = f_kin_anal.trajnodes2trajclusters(
  aux_list,kinlab.traj_nodes,num_nodes,new_num_frames,kinlab.particles)

# Output numc most important clusters
kinlab.dimensions = 1
# Need the following line to store links between clusters.
kinlab.network.clusters_links()
numc = 5
print "\nLargest clusters:"
for i in range(numc):
  print "cluster",i,\
    kinlab.network.cluster[i].label, \
    kinlab.network.cluster[i].weight / kinlab.network.weight
  lt_x,lt_y = kinlab.life_time(
    traj='clusters',state=i,norm=True,verbose=True)
  # misc.writeOutToFile2D(lt_x,lt_y,'results.lifetime.'+str(i)+'.dat')
  clusterWeight = kinlab.network.cluster[i].weight
  for otherCluster,otherWeight in sorted(
    kinlab.network.cluster[i].link.iteritems(), 
    key=operator.itemgetter(1), reverse=True)[:numc]:
    if otherCluster == i:
      print "  stationary        : %7.4f" % \
        float(otherWeight / clusterWeight)
    else:
      print "  transition to %4d: %7.4f" % (int(otherCluster), 
        float(otherWeight / clusterWeight))
      lt_x,lt_y = kinlab.first_passage_time(traj='clusters',from_state=i,
        to_state=otherCluster,norm=False,verbose=True)
      # misc.writeOutToFile2D(lt_x,lt_y,'results.fpt.'+
      #   str(i)+'-'+str(otherCluster)+'.dat')


# Close trajectory
system.delete_traj()