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.