

In [1]:

    
import numpy as np
import os
import pyspawn
import h5py
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display
from pyspawn.plotting import traj_plot
import MDAnalysis as mda
import nglview as nv

Single simulation analysis

Here we create a fafile object that pulls the data from the sim.hdf5 file and outputs the arrays for plotting.



In [2]:

    
print "Currently in directory:", os.getcwd()









    



Currently in directory: /Users/Dmitry/PycharmProjects/pySpawn17/examples



In [3]:

    
# THIS IS THE ONLY PART OF THE CODE THAT NEEDS TO BE CHANGED
dir_name = "/Users/Dmitry/Documents/Research/MSU/4tce/cis/"
h5filename = "sim.1.hdf5"



In [23]:

    
os.chdir(dir_name)
an = pyspawn.fafile(h5filename)

an.fill_electronic_state_populations(column_filename="N.dat")
an.fill_labels()
an.fill_istates()
an.get_numstates()

times = an.datasets["quantum_times"]
el_pop = an.datasets["electronic_state_populations"]
istates = an.datasets["istates"]
labels = an.datasets["labels"]
ntraj = len(an.datasets["labels"])
nstates = an.datasets['numstates']
an.fill_nuclear_bf_populations()

# write files with energy data for each trajectory
an.fill_trajectory_energies(column_file_prefix="E")

# write file with time derivative couplings for each trajectory
an.fill_trajectory_tdcs(column_file_prefix="tdc")

# compute Mulliken population of each trajectory
an.fill_mulliken_populations(column_filename="mull.dat")

mull_pop = an.datasets["electronic_state_populations"]

# istates dict
an.create_istate_dict()
istates_dict = an.datasets['istates_dict']









    



{'00b2': 0, '00b3': 0, '00': 1, '00b1': 0, '00b0': 0}

This part takes care of xyz trajectory files, bonds, angles (need to have atoms array for this to work)



In [5]:

    
# writing xyz files
an.write_xyzs()

# list of bonds to keep track of
bonds_list = np.array([[3, 11]])
# write datasets for bonds
an.fill_trajectory_bonds(bonds_list, column_file_prefix="bonds")

# dihedral angles list
diheds_list = np.array([[2, 6, 9, 10]])
# write datasets for dihedral angles
an.fill_trajectory_diheds(diheds_list, column_file_prefix="diheds")

Loading the arrays for plotting



In [7]:

    
arrays = ("poten", "pop", "toten", "aven", "kinen", "time", "tdc", "bonds", "diheds")
# creating dictionary for the datasets we want to plot
# keys are trajectory labels
for array in arrays:
    exec(array + "= dict()")

for traj in an.datasets["labels"]:
    
    poten[traj] = an.datasets[traj + "_poten"]
    toten[traj] = an.datasets[traj + "_toten"]
    kinen[traj] = an.datasets[traj + "_kinen"]
    time[traj] = an.datasets[traj + "_time"]
    tdc[traj] = an.datasets[traj + "_tdc"]
    bonds[traj] = an.datasets[traj + "_bonds"]
    diheds[traj] = an.datasets[traj + "_diheds"]

Setting plotting parameters (Perhaps there is a better way to do it, right now these hardcoded color and styles limit us to 7 electronic states and 16 trajectories. However, one could argue that more lines on a single plot would not be very informative anyway)



In [8]:

    
colors = ("r", "g", "b", "m", "y", "k", "k")
linestyles = ("-", "--", "-.", ":","-","-","-","-","-","-","-","-","-","-","-","-")
markers=("None","None","None","None","d","o","v","^","s","p","d","o","v","^","s","p")
large_size = 20
medium_size = 18
small_size = 16

This widget picks the trajectories we want to plot in case there are many of them



In [11]:

    
labels_to_plot_widget = widgets.SelectMultiple(
    options=labels,
    value=['00'],
    rows=10,
    description='Trajectories',
    disabled=False
)



In [12]:

    
display(labels_to_plot_widget)

Plotting Total Energies



In [11]:

    
%matplotlib notebook
traj_plot.plot_total_energies(time, toten, labels_to_plot_widget.value, istates_dict, colors, markers, linestyles)

Total Population



In [13]:

    
populated_states = np.amax(istates) + 1
traj_plot.plot_total_pop(times, mull_pop, populated_states, colors)

Plotting Potential Energies



In [14]:

    
display(labels_to_plot_widget)



In [15]:

    
%matplotlib notebook
traj_plot.plot_energies(labels_to_plot_widget.value, time, poten, nstates, colors, linestyles)

Plotting Energy gaps



In [16]:

    
display(labels_to_plot_widget)



In [22]:

    
%matplotlib notebook
# Gap between ground and first excited states
state1 = 0
state2 = 1
traj_plot.plot_e_gap(time, poten, labels_to_plot_widget.value, state1, state2, istates_dict,
                     colors, linestyles, markers)



In [ ]:

    
%matplotlib notebook
spawnthresh = 0.00785
# plot_tdc(labels, time, tdc, nstates, spawnthresh)
traj_plot.plot_tdc(time, tdc, labels_to_plot_widget.value,
                   nstates, istates_dict, spawnthresh, colors, linestyles, markers)

Bonds



In [19]:

    
display(labels_to_plot_widget)



In [21]:

    
%matplotlib notebook
traj_plot.plot_bonds(time, labels_to_plot_widget.value, bonds_list, bonds, colors, linestyles)

Dihedral Angles



In [ ]:

    
%matplotlib notebook
traj_plot.plot_diheds(time, labels_to_plot_widget.value, diheds_list, diheds, colors, linestyles)
plt.savefig('/Users/Dmitry/Documents/Research/MSU/4tce/cis/angles.png', dpi=300)

Trajectory visualization

In this widget we pick the trajectory label to visualize



In [ ]:

    
xyz_widget = widgets.RadioButtons(
    options=labels,
#     value='pineapple',
    description='Trajectory:',
    disabled=False
)
display(xyz_widget)



In [ ]:

    
print "Trajectory:", xyz_widget.value
path_to_xyz = dir_name + "/traj_" + xyz_widget.value + ".xyz"
print "Path to xyz file:", path_to_xyz
traj = mda.Universe(path_to_xyz)
w = nv.show_mdanalysis(traj)
w



In [ ]:



In [ ]: