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.
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