Getting Started
Constructing the System object
The System
object of PyOpenCAP contains the geometry and basis set information, as well
as the overlap matrix. The constructor takes in a Python dictionary as an argument,
and understands a specific set of keywords . There are four equivalent ways of specifying
the geometry and basis set: qchem_fchk, rassi_h5, molden, and inline. Here, we’ll use the rassi_h5 file.
sys_dict = {"molecule": "rassi_h5","basis_file": "path/to/rassi/file.h5"}
s = pyopencap.System(sys_dict)
smat = s.get_overlap_mat()
Constructing the CAP object
The CAP matrix is computed by the CAP
object. The constructor
requires a System
object, a dictionary containing the CAP parameters,
and the number of states.
nstates = 10
cap_dict = {"cap_type": "box",
"cap_x":"2.76",
"cap_y":"2.76",
"cap_z":"4.88",
"Radial_precision": "14",
"angular_points": "110"}
pc = pyopencap.CAP(s,cap_dict,nstates)
Parsing electronic structure data from file
The read_data()
function can read in the effective Hamiltonian
and densities in one-shot when passed a Python dictionary with the right keywords.
es_dict = {"method" : "ms-caspt2",
"package": "openmolcas",
"molcas_output":"path/to/output/file.out",
"rassi_h5": "path/to/rassi/file.h5"}
pc.read_data(es_dict)
h0 = pc.get_H()
Passing densities in RAM
Alternatively, one can load in the densities one at a time using the add_tdms()
or add_tdm()
functions. The density matrices should be in atomic orbital basis,
with the same atomic orbital ordering as the System
(which can be verify using
check_overlap_matrix
). The example below shows how one might pass the densities
from a PySCF calculation:
s.check_overlap_mat(pyscf_smat,"pyscf")
pc = pyopencap.CAP(s,cap_dict,10)
for i in range(0,10):
for j in range(i,10):
dm1 = fs.trans_rdm1(fs.ci[i],fs.ci[j],myhf.mo_coeff.shape[1],mol.nelec)
dm1_ao = np.einsum('pi,ij,qj->pq', myhf.mo_coeff, dm1, myhf.mo_coeff.conj())
pc.add_tdm(dm1_ao,i,j,"pyscf")
if i!=j:
pc.add_tdm(dm1_ao,j,i,"pyscf")
Once all of the densities are loaded, the CAP matrix is computed
using the compute_projected_cap()
function. The matrix can be retrieved using the
get_projected_cap()
function.
>>> pc.compute_projected_cap()
>>> W_mat=pc.get_projected_cap()
We now have our zeroth order Hamiltonian (stored in h0) and our CAP matrix(W_mat) in the state basis. Extracting resonance position and width requires analysis of the eigenvalue trajectories.
Analysis
PyOpenCAP provides user friendly tools for analysis of eigenvalue trajectories.
The CAPHamiltonian
contains functions aimed at diagonalization
of the CAP Hamiltonian over a range of eta values. Assuming one has already obtained H0 and
W in the state basis as numpy matrices, it can be constructed as such:
from pyopencap.analysis.CAPTrajectory import CAPHamiltonian
eta_list = np.linspace(0,5000,101)
eta_list = np.around(eta_list * 1E-5,decimals=5)
CAPH = CAPHamiltonian(H0=h0,W=mat)
# equivalently
CAPH = CAPHamiltonian(pc=pc)
CAPH.run_trajectory(eta_list,cap_lambda=0.0)
One can easily plot the eigenvalue spectrum in au or eV (relative to a given reference energy) as follows:
# total energies
plt.plot(np.real(CAPH.total_energies),np.imag(CAPH.total_energies),'ro')
plt.show()
# excitation energies
plt.plot(np.real(CAPH.energies_ev(ref_energy)),np.imag(CAPH.energies_ev(ref_energy)),'ro')
plt.show()
To analyze a given trajectory, use track_state()
traj = CAPH.track_state(1,tracking="overlap")
traj is now a EigenvalueTrajectory
object, which
contains helpful functions for analysis. One can plot raw and corrected trajectories:
plt.plot(np.real(traj.energies_ev(ref_energy)),np.imag(traj.energies_ev(ref_energy)),'-ro')
plt.plot(np.real(traj.energies_ev(ref_energy,corrected=True)),np.imag(traj.energies_ev(ref_energy,corrected=True)),'-bo')
There are also functions to help find the optimal value of the CAP strength parameter (and therefore, best estimate of resonance position and width) for uncorrected/corrected trajectories:
uc_energy,uc_eta_opt = traj.find_eta_opt()
corr_energy,corr_eta_opt = traj.find_eta_opt(corrected=True)
For more information, please see the documentation for the CAPHamiltonian
and EigenvalueTrajectory
classes.
See more
Please see the notebooks in our repository for detailed examples which demonstrate the full functionality of PyOpenCAP.