Load some pymatgen specific classes and module
In [1]:
import pymatgen as mg
from pymatgen.io.gaussianio import GaussianInput, GaussianOutput
Image class of IPython for rendering pictures.
In [2]:
from IPython.display import Image
Hereafter is the docstring of the GaussianInput class.
In [3]:
help(GaussianInput)
In this example, we get a $PO_2N^{2-}$ molecule and look for the best place for two Li atoms in the plane of the molecule. This example shows how to set up a Gaussian input file for all possible combinations. Coordinates of this molecule are available in xyz format :
In [4]:
mol = mg.Molecule.from_file("PO2N.xyz")
print(mol.to("xyz"))
Image("PO2N.png")
Out[4]:
First we set up the three guessed positions for Li atoms in the plane of the molecule. Each position is along the bissector of two PO or PN bonds. We build a list of pymatgen.Site object which represent an atom with its symbol, coordinates ....
In [5]:
posLi = [mg.Site("Li", [-0.424, -0.001, 0.267]),
mg.Site("Li", [-3.789, -2.031, -0.046]),
mg.Site("Li", [-3.787, 2.032, -0.044])]
print(posLi)
Now for each combinaison of two Li atoms among three possibilities we will write a gaussian input file. Look at the combinations method of itertools module.
In [6]:
from itertools import combinations
for sites in combinations(posLi, 2):
print(sites, "\t indices : ", [posLi.index(site) for site in sites])
Hereafter is the loop over combinations which will write all input files.
In [7]:
# Gaussian keywords
route_parms = {"opt": "tight", "Freq": ""}
link0 = {"%Nproc":"5"}
DFT = "B3LYP"
# the loop
for sites in combinations(posLi, 2):
# load the molecule
mol = mg.Molecule.from_file("PO2N.xyz")
# set up a label according to the combination
title = "combination_" + "".join([str(posLi.index(site)) for site in sites])
print(title)
# add the two Li atoms to the molecule
for site in sites:
mol.append(site.specie, site.coords)
# set up the calculation
gau = GaussianInput(mol, charge=0, spin_multiplicity=1, title=title,
functional=DFT, route_parameters=route_parms, link0_parameters=link0)
gau.write_file(title + ".com", cart_coords=True)
We can also load an input file and get some information.
In [8]:
gau = GaussianInput.from_file("combination_01.com")
print(gau.basis_set, gau.functional)
In [9]:
print(gau.to_string(cart_coords=True))
In [10]:
help(GaussianOutput)
Used the GaussianOutput class to read a Gaussian output file.
In [11]:
logfile = GaussianOutput("O3P-O-PO3/config_0123.log")
Is termination Normal ?
In [12]:
logfile.properly_terminated
Out[12]:
Display final energy or mulliken charges :
In [13]:
logfile.final_energy
Out[13]:
In [14]:
logfile.Mulliken_charges
Out[14]:
You can extract the coordinates of the final structure. All structures are pymatgen.Molecule objects, see here or the reference guide. You can do lot of stuff with this object such as : neighbors list, compute distances or the distance matrix, symmetry operation, atomic substitution ...
In [15]:
print(logfile.final_structure.to(fmt="xyz"))
You would prefer a gaussian input format with a z-matrix of the 42th geometrical optimization step.
In [16]:
print(logfile.structures[41].to(fmt="gjf"))
You can plot the geometrical convergence (here using plotly).
In [17]:
import plotly.plotly as py
import plotly.graph_objs as go
In [21]:
fig = go.Figure(
data=go.Data([
go.Scatter(x=[i for i in range(len(logfile.energies))], y=logfile.energies)
]),
layout=go.Layout(xaxis=go.XAxis(title="step"), yaxis=go.YAxis(title="Energy (ua)"))
)
py.image.ishow(fig)