Setup


In [ ]:
%pylab --no-import-all inline
import matplotlib as mpl

Don't necessarily actually want the figures to be inline (not really necessary here as I really just want to save them)


In [ ]:
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

In [ ]:
import paratemp.coordinate_analysis as ca

from paratemp import cd
import pandas as pd

import os

import MDAnalysis
import MDAnalysis as mda

In [ ]:
configs = {'MaEn': 'major-endo/13-3htmf-etc/05',
           'MaEx': 'major-exo/13-3htmf-etc/05',
           'MiEn': 'minor-endo/13-3htmf-etc/05',
           'MiEx': 'minor-exo/13-3htmf-etc/05'
          }

Combined open/closed FESs

Production


In [ ]:
for i in xrange(16):
    comb_df = pd.DataFrame(columns=['O-O', 'O(l)-Cy', 'O(r)-Cy'])
    if os.path.exists('comb-o-data-{}.h5'.format(i)):
        with pd.HDFStore('comb-o-data-{}.h5'.format(i)) as store:
            try:
                comb_df = store['time_750ns']
            except KeyError:
                comb_df = store['750ns']
    else:
        for key in configs:
            with cd(configs[key]):
                with open('TOPO/temperatures.dat', 'r') as t_file:
                    temps = list(t_file.read()[1:-2].split(', '))
                temp = float(temps[i])
                print 'Now starting on {} {}...'.format(key, i)
                univ = ca.Taddol('../../../solutes.gro', 
                                 'npt-PT-{}-out{}.xtc'.format(key, i))
                try:
                    univ.read_data()
                except IOError:
                    save_data = True
                    pass
                else:
                    save_data = False
                comb_df = comb_df.append(univ.ox_dists, ignore_index=True)
                if save_data:
                    univ.save_data()
        with pd.HDFStore('comb-o-data-{}.h5'.format(i)) as store:
            store['time_750ns'] = comb_df
        print('Saved combined data as {}'.format('comb-o-data-{}.h5'.format(i)))
    fig = univ.fes_ox_dists(data=comb_df, temp=temp, linewidth=3)
    for axes in fig.axes[:3]:
        axes.set_ylim[0:8]
    fig.savefig('fes-ox-dists-rjm-PT-comb-{}.pdf'.format(i))
    plt.close('all')
    print('Saved figure as {}'.format('fes-ox-dists-rjm-PT-comb-{}.pdf'.format(i)))
    print('Done with temp {:.0f} K ({} of 16).'.format(temp, i+1) \
          ' Moving to next temperature...\n\n')

Create final figures

Production


In [ ]:
overwrite = True

for key in configs:
    with cd(configs[key]):
        with open('TOPO/temperatures.dat', 'r') as t_file:
            temps = list(t_file.read()[1:-2].split(', '))
        univ = ca.Taddol('../../../solutes.gro', 
                         'npt-PT-{}-out{}.xtc'.format(key, '0'))
        final_time = int(univ.data['Time'].iat[-1]/1000)
        if final_time >= 1000:
            final_time = str(final_time/1000)+'us'
        else:
            final_time = str(final_time) + 'ns'
        if os.path.isfile('fes-ox-dists-rjm-PT-{}-comb-{}.pdf'.format(key, final_time)) or not overwrite:
            print('\nAll of {} seems to be complete, moving on to next config...\n\n'.format(key))
            continue
        fig_fes_cvs, axes_fes_cvs = plt.subplots(4, 4, sharex=True, sharey=True)
        fig_fes_cvs_open, axes_fes_cvs_open = plt.subplots(4, 4, sharex=True, sharey=True)
        fig_fes_cvs_closed, axes_fes_cvs_closed = plt.subplots(4, 4, sharex=True, sharey=True)
        fig_fes_ox_dists, axes_fes_ox_dists = plt.subplots(8, 8, sharex=True, sharey=True)
        for i in xrange(16):
            print 'Now starting on {} {}...'.format(key, i)
            temp = float(temps[i])
            univ = ca.Taddol('../../../solutes.gro', 
                             'npt-PT-{}-out{}.xtc'.format(key, i))
            final_time = str(int(univ.data['Time'].iat[-1]/1000))+'ns'
            file_name_end = '-rjm-PT-{}-{}-{}.pdf'.format(key, i, final_time)
            try:
                univ.read_data()
            except IOError:
                univ.calculate_distances()
                univ.save_data()
            univ.calc_open_closed()  # This is fast, so doesn't really matter if repeated
            print 'Done importing and calculating, making figures...'
            
            ax = axes_fes_cvs.flat[i]
            univ.fes_2d_cvs(temp=temp, display=False, ax=ax)
            filename = 'fes-cvs'+file_name_end
            if (not os.path.isfile(filename)) or overwrite:
                fig = univ.fes_2d_cvs(temp=temp)
                fig.savefig(filename)
            
            if len(univ.data[univ.data['open_TAD']]['CV1']) > 1:
                ax = axes_fes_cvs_open.flat[i]
                univ.fes_2d_cvs(univ.data[univ.data['open_TAD']]['CV1'], 
                                univ.data[univ.data['open_TAD']]['CV2'],
                                temp=temp, display=False, ax=ax)
                filename = 'fes-cvs-open'+file_name_end
                if (not os.path.isfile(filename)) or overwrite:
                    fig = univ.fes_2d_cvs(univ.data[univ.data['open_TAD']]['CV1'], 
                                          univ.data[univ.data['open_TAD']]['CV2'],
                                          temp=temp)
                    fig.savefig(filename)
            else:
                print '\n\n!!!! Not enough data to plot for open ' \
                      '{} {} !!!!\n\n'.format(key, i)
            
            if len(univ.data[univ.data['closed_TAD']]['CV1']) > 1:
                ax = axes_fes_cvs_closed.flat[i]
                univ.fes_2d_cvs(univ.data[univ.data['closed_TAD']]['CV1'], 
                                univ.data[univ.data['closed_TAD']]['CV2'],
                                temp=temp, display=False, ax=ax)
                filename = 'fes-cvs-closed'+file_name_end
                if (not os.path.isfile(filename)) or overwrite:
                    fig = univ.fes_2d_cvs(univ.data[univ.data['closed_TAD']]['CV1'], 
                                          univ.data[univ.data['closed_TAD']]['CV2'],
                                          temp=temp)
                    fig.savefig(filename)
            else:
                print '\n\n!!!! Not enough data to plot for closed ' \
                      '{} {} !!!!\n\n'.format(key, i)
            
            axes = axes_fes_ox_dists.flat[4*i:4*i+4]
            univ.fes_ox_dists(temp=temp, save=False, display=False, linewidth=3,
                              axes=axes)
            filename = 'fes-ox-dists'+file_name_end
            if (not os.path.isfile(filename)) or overwrite:
                fig = univ.fes_ox_dists(temp=temp, save=False, display=True, 
                                        linewidth=3)
                fig.savefig(filename)
            print('Done making and saving figures for '
                  '{} {}, closing and moving on...'.format(key, i))
        for axes in axes_fes_cvs[0:3]:
            for ax in axes:
                ax.set_xlabel('')
        for axes in axes_fes_cvs[:,1:]:
            for ax in axes:
                ax.set_ylabel('')
        fig_fes_cvs.tight_layout()
        for axes in axes_fes_cvs_open[0:3]:
            for ax in axes:
                ax.set_xlabel('')
        for axes in axes_fes_cvs_open[:,1:]:
            for ax in axes:
                ax.set_ylabel('')
        fig_fes_cvs_open.tight_layout()
        for axes in axes_fes_cvs_closed[0:3]:
            for ax in axes:
                ax.set_xlabel('')
        for axes in axes_fes_cvs_closed[:,1:]:
            for ax in axes:
                ax.set_ylabel('')
        fig_fes_cvs_closed.tight_layout()
        print '\n\n---Done with all {}, moving to next config'.format(key)
        plt.close('all')

Combined Ox FESs Production


In [ ]:
for i in xrange(16):
    comb_df = pd.DataFrame(columns=['O-O', 'O(l)-Cy', 'O(r)-Cy'])
    if os.path.exists('comb-o-data-{}.h5'.format(i)):
        with pd.HDFStore('comb-o-data-{}.h5'.format(i)) as store:
            comb_df = store['time_1us']
    else:
        for key in configs:
            with cd(configs[key]):
                with open('TOPO/temperatures.dat', 'r') as t_file:
                    temps = list(t_file.read()[1:-2].split(', '))
                temp = float(temps[i])
                print 'Now starting on {} {}...'.format(key, i)
                univ = ca.Taddol('../../../solutes.gro', 
                                 'npt-PT-{}-out{}.xtc'.format(key, i))
                try:
                    univ.read_data()
                except IOError:
                    save_data = True
                    pass
                else:
                    save_data = False
                comb_df = comb_df.append(univ.ox_dists, ignore_index=True)
                if save_data:
                    univ.save_data()
        with pd.HDFStore('comb-o-data-{}.h5'.format(i)) as store:
            store['time_1us'] = comb_df
        print('Saved combined data as {}'.format('comb-o-data-{}.h5'.format(i)))
    fig = univ.fes_ox_dists(data=comb_df, temp=temp, linewidth=3)
    for axes in fig.axes[:3]:
        axes.set_ylim([0,8])
    fig.savefig('fes-ox-dists-rjm-PT-comb-{}.pdf'.format(i))
    plt.close('all')
    print('Saved figure as {}'.format('fes-ox-dists-rjm-PT-comb-{}.pdf'.format(i)))
    print('Done with temp {:.0f} K ({} of 16).'.format(temp, i+1) +
          ' Moving to next temperature...\n\n')

CV Cut Plots

Production


In [ ]:
r = 0.0019872  # kcal_th/(K mol)

for key in configs:
    with cd(configs[key]):
        i = 0
        print 'Now starting on {} {}...'.format(key, i)
        univ = ca.Taddol('../../../solutes.gro', 
                         'npt-PT-{}-out{}.xtc'.format(key, i))
        temp = 205
        final_time = int(univ.data['Time'].iat[-1]/1000)
        if final_time > 1000:
            final_time = str(final_time/1000)+'us'
        else:
            final_time = str(final_time) + 'ns'
        file_name_end = '-rjm-PT-{}-{}-{}.pdf'.format(key, i, final_time)
        try:
            univ.read_data()
        except IOError:
            univ.calculate_distances()
            univ.calc_open_closed()
            univ.save_data()
        univ.calc_open_closed()

        x = univ.data.loc[lambda x: (2.5 < x.CV2) & (x.CV2 < 3.5) & x.closed_TAD]['CV1']
        if len(x) > 2:
            n, bins = np.histogram(x, bins=20)
            n = [float(j) for j in n]
            prob = np.array([j / max(n) for j in n]) + 1e-40
            delta_g = np.array([-r * temp * np.log(p) for p in prob])
            delta_g
            fig, ax = plt.subplots()
            line, = ax.plot(bins[:-1], delta_g)
            ax.set_ylabel(r'$\Delta G$ / (kcal / mol)')
            ax.set_xlabel(r'CV 1 / $\mathrm{\AA}$')
            fig.tight_layout()
            fig.savefig('fes-CV1-closed-shortCV2'+file_name_end)
        else:
            print('Not enough closed frames for {} {}'.format(key, i))
            
        x = univ.data.loc[lambda x: (7.5 < x.CV2) & (x.CV2 < 8.5) & x.open_TAD]['CV1']
        if len(x) > 2:
            n, bins = np.histogram(x, bins=20)
            n = [float(j) for j in n]
            prob = np.array([j / max(n) for j in n]) + 1e-40
            delta_g = np.array([-r * temp * np.log(p) for p in prob])
            delta_g
            fig, ax = plt.subplots()
            line, = ax.plot(bins[:-1], delta_g)
            ax.set_ylabel(r'$\Delta G$ / (kcal / mol)')
            ax.set_xlabel(r'CV 1 / $\mathrm{\AA}$')
            fig.tight_layout()
            fig.savefig('fes-CV1-open-longCV2'+file_name_end)
        else:
            print('Not enough open frames for {} {}'.format(key, i))

Select geoms from minima

Production


In [ ]:
cv1_cuts = [6.5, 9.]
cv2_cuts = [1.5, 3.]
name_set = 'lCV1-sCV2'

univ = ca.Taddol('solutes.gro', 'major-endo/13-3htmf-etc/05/pbc-MaEn-0.xtc')
try:
    univ.data['CV1']
except KeyError:
    univ.read_data(filename='major-endo/13-3htmf-etc/05/npt-PT-MaEn-out0.h5')

bool_array = ((univ.data['CV1'] > cv1_cuts[0]) & (univ.data['CV1'] < cv1_cuts[1]) 
              & (univ.data['CV2'] > cv2_cuts[0]) & (univ.data['CV2'] < cv2_cuts[1]))
num = len(univ.data[bool_array])
print('These cutoffs include {} frames.'.format(num))

solutes = univ.select_atoms('resname is 3HT or resname is CIN or resname is TAD')

with mda.Writer('minim-structs-'+name_set+'-rjm-PT-MaEn-0.xtc', 
                solutes.n_atoms) as W:
    for ts in univ.trajectory:
        if bool_array[univ.trajectory.frame]:
            W.write(solutes)

In [ ]:
cv1_cuts = [1.5, 4.]
cv2_cuts = [6.75, 8.5]
name_set = 'sCV1-lCV2'

univ = ca.Taddol('solutes.gro', 'major-endo/13-3htmf-etc/05/pbc-MaEn-0.xtc')
try:
    univ.data['CV1']
except KeyError:
    univ.read_data(filename='major-endo/13-3htmf-etc/05/npt-PT-MaEn-out0.h5')

bool_array = ((univ.data['CV1'] > cv1_cuts[0]) & (univ.data['CV1'] < cv1_cuts[1]) 
              & (univ.data['CV2'] > cv2_cuts[0]) & (univ.data['CV2'] < cv2_cuts[1]))
num = len(univ.data[bool_array])
print('These cutoffs include {} frames.'.format(num))

solutes = univ.select_atoms('resname is 3HT or resname is CIN or resname is TAD')

with mda.Writer('minim-structs-'+name_set+'-rjm-PT-MaEn-0.xtc', 
                solutes.n_atoms) as W:
    for ts in univ.trajectory:
        if bool_array[univ.trajectory.frame]:
            W.write(solutes)

In [ ]:
cv1_cuts = [1.5, 4.]
cv2_cuts = [6.75, 8.5]
name_set = 'sCV1-lCV2'

univ = ca.Taddol('solutes.gro', 'major-exo/13-3htmf-etc/05/pbc-MaEx-0.xtc')
try:
    univ.data['CV1']
except KeyError:
    univ.read_data(filename='major-exo/13-3htmf-etc/05/npt-PT-MaEx-out0.h5')

bool_array = ((univ.data['CV1'] > cv1_cuts[0]) & (univ.data['CV1'] < cv1_cuts[1]) 
              & (univ.data['CV2'] > cv2_cuts[0]) & (univ.data['CV2'] < cv2_cuts[1]))
num = len(univ.data[bool_array])
print('These cutoffs include {} frames.'.format(num))

solutes = univ.select_atoms('resname is 3HT or resname is CIN or resname is TAD')

with mda.Writer('minim-structs-'+name_set+'-rjm-PT-MaEx-0.xtc', 
                solutes.n_atoms) as W:
    for ts in univ.trajectory:
        if bool_array[univ.trajectory.frame]:
            W.write(solutes)

In [ ]:
cv1_cuts = [6.5, 9.]
cv2_cuts = [1.5, 3.]
name_set = 'lCV1-sCV2'

univ = ca.Taddol('solutes.gro', 'major-exo/13-3htmf-etc/05/pbc-MaEx-0.xtc')
try:
    univ.data['CV1']
except KeyError:
    univ.read_data(filename='major-exo/13-3htmf-etc/05/npt-PT-MaEx-out0.h5')

bool_array = ((univ.data['CV1'] > cv1_cuts[0]) & (univ.data['CV1'] < cv1_cuts[1]) 
              & (univ.data['CV2'] > cv2_cuts[0]) & (univ.data['CV2'] < cv2_cuts[1]))
num = len(univ.data[bool_array])
print('These cutoffs include {} frames.'.format(num))

solutes = univ.select_atoms('resname is 3HT or resname is CIN or resname is TAD')

with mda.Writer('minim-structs-'+name_set+'-rjm-PT-MaEx-0.xtc', 
                solutes.n_atoms) as W:
    for ts in univ.trajectory:
        if bool_array[univ.trajectory.frame]:
            W.write(solutes)