Skip to content

Instantly share code, notes, and snippets.

@schneiderfelipe
Last active August 23, 2018 03:28
Show Gist options
  • Save schneiderfelipe/11e364daf626de6754f4f6cbd710af26 to your computer and use it in GitHub Desktop.
Save schneiderfelipe/11e364daf626de6754f4f6cbd710af26 to your computer and use it in GitHub Desktop.
Use cases for pyrrole
#!/usr/bin/python3
import pandas as pd
import pyrrole as pyr
from matplotlib import pyplot as plt
# if you want a single molecule
# pyr.Atoms(series)
# pyr.Atoms(dataframe["A"])
# pyr.atoms.from_cclib(ccopen("A.out")) # Atoms.data index will be jobfilename
# if you want all your molecules
# df = ccframe(...)
# molecules = pyr.atoms.from_dataframe(dataframe) # returns a list of pyr.Atoms, one for each row in df
# if you want an equation
# df = pd.DataFrame([m.to_series() for m in molecules])
# equation = pyr.ChemicalEquation("A -> B", data=dataframe)
# equation.calculate("freeenergy") # things are always calculated on demand such that equations are always stateless
# if you want a system
# system = pyr.ChemicalSystem(..., data=dataframe)
def main():
system = pyr.ChemicalSystem(["BH <- AH <=> A- + H+",
"A- -> B-",
"BH <=> B- + H+"])
# $ cat reactions.txt
# BH <- AH <=> A- + H+
# A- -> B-
# BH <=> B- + H+
# with open("reactions.txt") as f:
# system = pyr.ChemicalSystem(f.readlines()) # if you have a file for the model
# $ cat data.tsv
# label freeenergy
# A- 0.0
# B- 1.0
# AH 2.0
# BH 3.0
# H+ 4.0
# internal data is almost aways pd.Series or pd.DataFrame and we say shamelessly so:
# pass things that can be at least converted!
df = pd.read_table("data.tsv").set_index("label")
# if species is a dict of labels to paths, i.e., {"A": "A.out", ...}
# from cclib import ccopen
# df = pd.DataFrame({label: pyr.atoms.from_cclib(ccopen(path)).to_series()
# for label, path in species.items()})
# in the code, any attribute named "data" will be a species (row) by property (column) dataframe
# only accept dataframes, can also be passed to constructor as ChemicalSystem(..., data=df),
# delegates data to each ChemicalEquation.data
system.data = df
# https://en.wikipedia.org/wiki/Energy_profile_(chemistry)
# calls system.all_simple_paths(source, target), which in turn uses networkx,
# and assumes energies are as in cclib ("freeenergy", ..., and in hartree)
system.energy_profile("AH", "BH")
plt.show()
print(system.to_latex()) # calls system.to_dataframe(), which calls ChemicalEquation.to_series() to each equation
# x, y = system.spectrum("ir")
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment