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/(X)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.

Preliminary: Prepare input orbitals¶

As with any multi-reference calculation, the choice of active space is crucial for CAP/(X)MS-CASPT2, and is most often guided by chemical intuition. We refer the reader to the OpenMolcas manual for how to prepare input orbitals for a state-averaged RASSCF calculation.

Step 1: Running the OpenMolcas calculation¶

State-averaged RASSCF

In order to utilize the Perturbative CAP approach, a multi-state excited state calculation must be performed. In the RASSCF module, the keyword ‘CIROOT’ is used to activate state-averaged RASSCF calculations.

&RASSCF
CIROOT = 10 10 1

Export transition densities with RASSI

To generate the one-particle densities required to construct the CAP matrix, the RASSI module must be executed with the TRD1 keyword activated. This keyword saves one-particle transition density matrices between each pair of RASSCF states as well as the one-particle density matrices for each state to a file titled $Jobname.rassi.h5.

&RASSI
TRD1

Generate effective Hamiltonian with (X)MS-CASPT2

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

&CASPT2
Multistate = all
# or
Xmultistate = all

Reference energy

There are multiple strategies for obtaining the reference energy used to define the resonance position. For anionic resonances, one such strategy is to add an additional diffuse orbital to the active space in order to mimic ionization, which obtains the resonance and the ground state of the neutral molecule in a single calculation [Kunitsa2017]. Another strategy (which was used in the tutorial) is to calculate the ground state of the neutral molecule with CASCI/CASPT2 using the optimized orbitals of the anionic state.

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 = pycap.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 = pycap.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 = pycap.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"}
pc = pycap.CAP(my_system,cap_dict,10,"openmolcas")

Before we can compute the CAP matrix in the state basis, we must load in the density matrices. There are two ways of doing this. The first is to use the read_data() function. As shown below, we define a dictionary which contains the following keys: “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 the effective Hamiltonians from (X)MS-CASPT2 calculations can be parsed from an OpenMolcas output file. We note that when read_data() is used, our code symmetrizes the CAP matrix in the state basis.

es_dict = {"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()

Alternatively, one can load in the densities one at a time using add_tdm(). In our examples below, we load in the matrices from rassi.h5 using the h5py package, and then pass them as numpy arrays to the CAP object.

import h5py
f = h5py.File('path/to/rassi.h5', 'r')
dms = f["SFS_TRANSITION_DENSITIES"]
# spin traced
nbasis = int(np.sqrt(dms.shape[2]))
for i in range(0,10):
    for j in range(i,10):
        dm = 0.5*np.reshape(dms[i][j],(nbasis,nbasis))
        pc.add_tdm(dm,i,j,"openmolcas","path/to/rassi.h5")
        # usually a good idea to symmetrize
        if i!=j:
            pc.add_tdm(dm,,j,i,"openmolcas","path/to/rassi.h5")

Step 3: Computing the CAP matrix¶

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

pc.compute_perturb_cap()
W_mat=pc.get_perturb_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_perturb_cap().

pc.compute_ao_cap()
pc.renormalize()
pc.compute_perturb_cap()

Step 4: Generate eigenvalue trajectories¶

Extracting resonance position and width requires analysis of the eigenvalue trajectories. Template scripts are provided in the repository. Development of automated tools for trajectory analysis is a subject of future work.

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
XMS-CASPT2

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/SS)DMRG-(PC/SC)NEVPT2
SS-CASPT2
MC-PDFT