1
#!/usr/bin/env python
2 100
"""
3
basic_plot.py
4

5
Basic plots (lines or bars) of information in SD tags.
6

7
Version:    Apr 12 2019
8
By:         Victoria T. Lim
9

10

11
TODO:
12
- plot diff tags in same file
13

14
"""
15 100
import numpy as np
16 100
import itertools
17

18 100
import matplotlib.pyplot as plt
19

20 100
import quanformer.proc_tags as pt
21 100
import quanformer.reader as reader
22

23

24 100
def basic_plot(infile, tag, style, molname=None, take_relative=False, har_to_kcal=False):
25
    """
26
    TODO
27

28
    Parameters
29
    ----------
30
    infile : string
31
        Name of SDF file with information in SD tags.
32
    tag : string
33
        Full tag string directly as listed in the SD file.
34
    style : string
35
        plot style. can be 'scatter', 'line', or 'bar'
36
        TODO
37
    take_relative : Boolean
38
        subtract lowest value
39
    har_to_kcal : Boolean
40
        multiply data in Hartrees by 627.5095 to yield kcal/mol
41

42
    """
43
    # Open molecule file.
44 100
    mols = reader.read_mols(infile)
45

46 100
    for i, mol_i in enumerate(mols):
47 100
        if molname is not None and mol_i.GetTitle() != molname:
48 0
            continue
49

50
        # get array of all conformer data of this mol
51 100
        try:
52 100
            data_array = np.fromiter(pt.get_sd_list(mol_i, datum='', taglabel=tag), dtype=np.float64)
53 0
        except ValueError:
54 0
            data_array = np.asarray([np.nan])*mol_i.NumConfs()
55

56
        # exclude conformers for which job did not finish (nan)
57 100
        nanIndices = np.argwhere(np.isnan(data_array))
58 100
        for j in reversed(nanIndices):  # loop in reverse to delete correctly
59 0
            data_array = np.delete(data_array, j)
60

61 100
        if take_relative:
62 100
            data_array = data_array - np.amin(data_array)
63 100
        if har_to_kcal:
64 100
            data_array = 627.5095*data_array
65

66
        # generate plot
67 100
        plt.plot(data_array)
68 100
        plt.grid()
69 100
        plt.title(mol_i.GetTitle()+'\n'+tag, fontsize=14)
70 100
        plt.savefig(f'output_{i}.png', bbox_inches='tight')
71 100
        plt.show()
72

73

74 100
def combine_files_plot(infile, figname='combined.png', molname=None, verbose=False, take_relative=False, har_to_kcal=False):
75
    """
76
    TODO
77

78
    This only supports plotting of ONE specified molecule across different files.
79

80
    Note on take_relative:
81
        [1] Subtracting global minimum (single value) from all energies
82
        doesn't work since everything is still on different scale.
83
    subtract: (1) first conformer of each?, (2) global minimum?, (3) minimum of each?
84

85
    Parameters
86
    ----------
87
    infile : str
88
        Filename with information on the files to read in, and
89
        the SDF tags to be extracted from each. Columns are:
90
        (1) QM method/basis, (2) sdf file, (3) tag key in sdf (like 'QM spe'),
91
        (4) arbitrary label for plotting. Separate columns by comma.
92
    molname
93
    verbose
94

95
    """
96 100
    wholedict = reader.read_text_input(infile)
97

98 100
    numFiles = len(wholedict)
99 100
    xarray = []
100 100
    yarray = []
101 100
    labels = []
102 100
    titles = []
103 100
    for i in wholedict:
104 100
        print("Reading molecule(s) from file: ", wholedict[i]['fname'])
105 100
        mols = reader.read_mols(wholedict[i]['fname'])
106 100
        qmethod, qbasis = reader.separated_theory(wholedict[i]['theory'])
107 100
        short_tag = wholedict[i]['tagkey']
108

109 100
        for j, mol_j in enumerate(mols):
110 100
            if molname is not None and mol_j.GetTitle() != molname:
111 100
                continue
112 100
            data_array = np.array(list(map(float,
113
                pt.get_sd_list(mol_j, short_tag, 'Psi4', qmethod, qbasis))))
114

115 100
        if take_relative:
116 100
            data_array = data_array - data_array[0]
117
            #data_array = data_array/data_array[0]
118 100
        if har_to_kcal:
119 100
            data_array = 627.5095*data_array
120

121 100
        titles.append(mol_j.GetTitle())
122 100
        labels.append(wholedict[i]['label'])
123 100
        yarray.append(data_array)
124 100
        xarray.append(range(len(data_array)))
125

126 100
    if verbose:
127 100
        header = '{}\n'.format(molname)
128 100
        for l in labels:
129 100
            header += ("%s\n" % l)
130 100
        xydata = np.vstack((xarray[0], yarray)).T
131 100
        np.savetxt('combined.dat', xydata, delimiter='\t', header=header,
132
            fmt=' '.join(['%i'] + ['%10.4f']*numFiles))
133

134
    # letter labels for x-axis
135 100
    num_confs = len(xarray[0])
136 100
    letters = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
137 100
    rpt = int((num_confs / 26) + 1)
138 100
    xlabs = [''.join(i)
139
             for i in itertools.product(letters, repeat=rpt)][:num_confs]
140

141 100
    fig = plt.figure()
142 100
    ax = fig.add_subplot(111)
143 100
    xlabel='conformer'
144 100
    ylabel="energy"
145

146
    # vtl print max range of relative energies
147 100
    conf_then_file = np.array(yarray).T
148 100
    ranges = []
149 100
    for c in conf_then_file:
150 100
        c_spread = max(c)-min(c)
151 100
        ranges.append(c_spread)
152 100
    print(f'mol {molname} max range: {max(ranges)}')
153

154 100
    ax.set_prop_cycle(plt.cycler('color', plt.cm.rainbow(np.linspace(0, 1, len(yarray)))))
155 100
    for i, (xs, ys) in enumerate(zip(xarray,yarray)):
156 100
        plt.plot(xs, ys, '-o', lw=0.8, label=labels[i])
157

158
    # publication view
159
#    plt.ylabel(ylabel,fontsize=8)
160
#    plt.xlabel(xlabel,fontsize=8)
161
#    plt.legend(bbox_to_anchor=(0.08,1.05),loc=3,fontsize=8)
162
#    fig.set_size_inches(3.37,1.7)
163

164
    # standard view
165 100
    plt.ylabel(ylabel,fontsize=14)
166 100
    plt.xlabel(xlabel,fontsize=14)
167 100
    plt.xticks(list(range(num_confs)), xlabs)
168 100
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
169

170 100
    plt.title(molname)
171 100
    plt.grid()
172 100
    plt.savefig(figname, bbox_inches='tight',dpi=300)
173 100
    plt.show()

Read our documentation on viewing source code .

Loading