OpenMolcas

OpenMolcas is an open-source quantum chemistry package which specializes in multiconfigurational approaches to electronic structure. OpenMolcas can be used in tandem with PyOpenCAP to perform complex absorbing potential (extended)multi-state complete active space second order perturbation theory [CAP/MS-CASPT2] calculations, which have been shown to yield accurate energies and lifetimes for metastable electronic states. Here, we outline the steps of performing these calculations using OpenMolcas and PyOpenCAP. Some suggested readings are provided at the bottom of the page.

Step 1: Running OpenMolcas calculation

To generate the one-particle densities required to construct the CAP matrix, the RASSI module must be executed with the TRD1 keyword activated. When using XMS-CASPT2, RMS-CASPT2, or other variants which utilize rotated CASSCF wave functions, the CAP matrix will eventually be rotated into the new basis using the rotation matrix in the output (U^dagger*W*U). RASSI will save transition density matrices between each pair of CASSCF states as well as the one-particle density matrices for each state to a file titled $Jobname.rassi.h5.

Export transition densities with RASSI

&RASSI
TRD1

Generate effective Hamiltonian with MS-CASPT2

The MS-CASPT2 approach is required to generate an appropriate zeroth Hamiltonian for the projected CAP method. To activate MS-CASPT2 in OpenMolcas, use the Multistate keyword in the CASPT2 module.

&CASPT2
Multistate = all

See the OpenMolcas manual for other variants of MS-CASPT2 which can be activated in the &CASPT2 section.

Step 2: Importing the data to PyOpenCAP

System object

To run a PyOpenCAP calculation, the geometry and basis set must be imported into a System object. The constructor takes in a Python dictionary as an argument. The relevant keywords are discussed here, and more information is provided in the keywords page.

Rassi.h5

The rassi.h5 file which contains the one-particle densities also contains the geometry and basis set information. To read in from rassi, “molcas_rassi” must set as the value to the key “molecule”, and the path to the file must be set as the value to the key “basis_file”. Here is an example:

sys_dict = {"molecule": "molcas_rassi","basis_file": "path/to/rassi.h5"}
my_system = pyopencap.System(sys_dict)

Molden

Molden files generated by OpenMolcas contain the geometry and basis set information. To read in from molden, “molden” must be set as the value to the key “molecule”, and the path to the file must be set as the value to the key “basis_file”. Here is an example:

sys_dict = {"molecule": "molden","basis_file": "path/to/file.molden"}
my_system = pyopencap.System(sys_dict)

Inline(not recommended)

The molecule and basis set can also be specified manually. 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. Up to G-type functions are supported.

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"}
my_system = pyopencap.System(sys_dict)

One particle densities/zeroth order Hamiltonian

The CAP matrix is computed by the CAP object. The constructor requires a System, a dictionary containing the CAP parameters, the number of states, and finally the string “openmolcas”, which denotes the ordering of the atomic orbital basis set. An example is provided below. Please see the keywords section for more information on the CAP parameters.

cap_dict = {"cap_type": "box",
            "cap_x":"2.76",
            "cap_y":"2.76",
            "cap_z":"4.88",
            "Radial_precision": "14",
            "angular_points": "110"}
nstates = 10
pc = pyopencap.CAP(my_system,cap_dict,nstates)

Before we can compute the CAP matrix in the state basis, we must load in the density matrices. The best way is to use the read_data() function. As shown below, we define a dictionary which contains the following keys: “package”(openmolcas), “method” (electronic structure method chosen), “rassi_h5”(density matrices), and “molcas_output”(output file containing effective Hamiltonian). The effective Hamiltonian can be retrieved using the get_H() function of the CAP object. Currently, only effective Hamiltonians from MS-CASPT2 calculations can be parsed from an OpenMolcas output file.

es_dict = { "package": "openmolcas",
       "method" : "ms-caspt2",
       "molcas_output":"path/to/output.out",
       "rassi_h5":"path/to/rassi.h5"}
pc.read_data(es_dict)
# save the effective Hamiltonian for later use
h0 = pc.get_H()

Step 3: Computing the CAP matrix

Once all of the densities are loaded, the CAP matrix is computed using compute_projected_cap(). The matrix can be retrieved using get_projected_cap().

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()
pc.compute_projected_cap()

Step 4: Generate and analyze eigenvalue trajectories

H0 and W, or the CAP object can be used to construct a CAPHamiltonian object.

from pyopencap.analysis import CAPHamiltonian
CAPH = CAPHamiltonian(H0=H0,W=W_mat)
# equivalently
CAPH = CAPHamiltonian(pc=pc)

See the analysis section for more details.

Officially supported methods

The following methods have been benchmarked, and the read_data() function is capable of parsing output files to obtain the zeroth order Hamiltonian.

MS-CASPT2, and other variants (e.g. XMS-CASPT2) which utilize unitary rotations of the original CASSCF states. The CAP matrix will be rotated into the new basis using the rotation matrix.

Untested (use at your own risk!)

The following methods are capable of dumping densities using the TRD1 keyword of the RASSI module, but have not been benchmarked for any systems, and the zeroth order Hamiltonian cannot be parsed from the output file using the read_data() function. Use at your own caution, and please contact us if you find success using any of these methods so we can add official support!

(QD)DMRG-(PC/SC)NEVPT2
SS-CASPT2
MC-PDFT