PySCF

PySCF is an ab initio computational chemistry program natively implemented in Python. The major advantage of using Pyscf in tandem with OpenCAP is that calculations can be performed in one-shot within the same python script. Since PySCF allows direct control over data structures such as density matrices, the interface between PySCF and OpenCAP is seamless.

Step 1: Defining the System object

Molden(recommended)

The best way to construct the System object is to import the geometry and basis set from a molden file generated by a PySCF. This ensures proper ordering of the AO basis set.

molden_dict = {"basis_file":"molden_in.molden","molecule": "molden"}
pyscf.tools.molden.from_scf(myhf,"molden_in.molden")
s = pyopencap.System(molden_dict)

Inline

The molecule and basis set can also be specified inline. The “molecule” keyword must be set to “read”, and then an additional keyword “geometry” must be specified, with a string that contains the geometry in xyz format. The “basis_file” keyword must be set to a path to a basis set file formatted in Psi4 style, which can be downloaded from the MolSSI BSE. Other optional keyword for this section include “bohr_coordinates” and “cart_bf”. Please see the keywords section for more details. It is recommended to check the overlap matrix to ensure that the ordering and normalization matches. Up to G-type functions are supported.

pyscf_smat = scf.hf.get_ovlp(mol)
sys_dict = {"geometry":    '''N  0  0   1.039
                          N  0  0   -1.039
                          X   0  0   0.0''',
                    "molecule" : "read",
                    "basis_file":"path/to/basis.bas",
                    "cart_bf":"d",
                    "bohr_coordinates:": "true"}
s.check_overlap_mat(pyscf_smat,"pyscf")

Step 1: Defining 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. An example is provided below. Please see the keywords section for more information on the CAP parameters.

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(my_system,cap_dict,nstates)

Step 2: Passing the density matrices

The simplest interface is with the FCI modules. Transition densities can be obtained using the trans_rdm1() function:

fs = fci.FCI(mol, myhf.mo_coeff)
e, c = fs.kernel()
# tdm between ground and 1st excited states
dm1 = fs.trans_rdm1(fs.ci[0],fs.ci[1],myhf.mo_coeff.shape[1],mol.nelec)

Importantly, trans_rdm1 returns the density matrix in MO basis. Thus before passing it to PyOpenCAP, it must be transformed into AO basis:

dm1_ao = np.einsum('pi,ij,qj->pq', myhf.mo_coeff, dm1, myhf.mo_coeff.conj())

Densities are loaded in one at a time using add_tdm(). Ensure that the indices of each state match those of the zeroth order Hamiltonian.

for i in range(0,len(fs.ci)):
    for j in range(0,len(fs.ci)):
        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")

Note:

The interface with PySCF is not restricted to the FCI module. The add_tdm() function is completely general; it requires only that the densities are in AO basis, and that the basis set ordering matches the system. Examples for ADC, EOM-EA-CCSD, and TDA-TDDFT are provided in the repository.

Step 3: Computing the CAP matrix

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()

Note:

When using cartesian d, f, or g-type basis functions, special care must be taken to ensure that the normalization conventions match what is used by OpenMolcas. In these cases, compute_ao_cap() and then renormalize() or renormalize_cap() should be invoked before calling compute_projected_cap().

pc.compute_ao_cap(cap_dict)
pc.renormalize_cap(pyscf_smat,"pyscf")
pc.compute_projected_cap()

Step 4: Generate and analyze eigenvalue trajectories

H0 and W can be used to construct a CAPHamiltonian object. In many cases, it can be advantageous to use the export() function, which generates an OpenCAP formatted output file, which can be used for later analysis.

from pyopencap.analysis import CAPHamiltonian
CAPH = CAPHamiltonian(H0=H0,W=W_mat)
CAPH.export("output.out")

See the analysis section for more details.

Officially supported methods

  • Full CI

  • ADC (through ADCC)

  • TDA-TDDFT

Untested (use at your own risk!)

Any module which one particle transition densities available can be supported. This includes all methods which can utilize the trans_rdm1 function, including but not limited to:

  • MRPT