Running vibronic calculations from a Molcas RASSI set of calculations

Requirements:

• Spin-free hamiltonian for every structure displaced along the normal modes for the numerical first derivative.

• Zero-order file with the following spin-free properties for each of the vibronic coupling calculations:

• Electric dipole moments ('electric_dipole')

• Magnetic dipole moments ('magnetic_dipole')

• Electric quadrupole moments ('electric_quadrupole')

• Zero-order file with the spin-free and spin-orbit energies.

The Zero-order file referenced above is the calculation at the equilibrium geometry with the EJOB option in the molcas input file.

Some imports for the resource to deal with the data files for the purpose of this example.

from vibrav.base import resource
import tarfile
import os
import glob
import shutil

tarball = resource('molcas-ucl6-2minus-vibronic-coupling.tar.xz')
with tarfile.open(tarball, 'r:xz') as tar:
    tar.extractall()
parent = os.getcwd()
os.chdir('molcas-ucl6-2minus-vibronic-coupling')

Begin the vibronic coupling calculations

The configuration file must be give as it will tell the program where to look for all of the things it needs. Refer to the vibrav.vibronic.vibronic.Vibronic documentation for more information about the required and default arguments.

from vibrav.vibronic import Vibronic

vib = Vibronic(config_file='va.conf')

Brief rundown of the most important input arguments

• property: Tell the program which vibronic property to calculate.

• write_oscil: If set to True the program will write the ocillator strengths for the x, y and z directions along with the averaged values in the vibronic-outputs directory under the filenames oscillators-1.txt, oscillators-2.txt, oscillators-3.txt and oscillators-0.txt, respectively.

• write_property: If set to True the program will write the spin orbit vibronic properties that are calculated into the directories under the names vib??? for the respective vibrational mode. It will then create two directories plus and minus referring to the transitions of dE + E_vib and dE - E_vib, respectively. Inside those directories it will write the spin-orbit property values for the three catrtesian directions.

Refer to the vibrav.vibronic.vibronic.Vibronic.vibronic_coupling documentation for more information on the input parameters.

vib.vibronic_coupling(property='electric-dipole', print_stdout=False, temp=298,
                      write_property=True, write_oscil=True, boltz_states=2,
                      verbose=False, eq_cont=False, select_fdx=-1)

/home/herbertl/github/vibrav/vibrav/numerical/boltzmann.py:114: Warning: Calculating only the first 2 states for the Boltzmann distribution.
  warnings.warn("Calculating only the first {} ".format(states) \

Directory listing of the vibronic-outputs directory.

os.listdir('vibronic-outputs')

['oscillators-0.txt',
 'oscillators-1.txt',
 'oscillators-3.txt',
 'boltzmann-populations.csv',
 'oscillators-2.txt',
 'alpha.txt']

The different directories created for each vibrational mode. Note, this uses a one based index as we typically reserve 0 for the normal mode coordinate at the equilibrium structure.

sorted(glob.glob('vib???'))

['vib001',
 'vib002',
 'vib003',
 'vib004',
 'vib005',
 'vib006',
 'vib007',
 'vib008',
 'vib009',
 'vib010',
 'vib011',
 'vib012',
 'vib013',
 'vib014',
 'vib015']

Directory listing of one of the vibronic coupling property output directories.

os.listdir(os.path.join('vib001', 'plus'))

['dipole-1.txt', 'dipole-2.txt', 'dipole-3.txt', 'energies.txt']

Reading the .txt files that are generated

There is a utility in the vibrav code that allows the user to easily open and read the .txt files that are generated by both molcas and vibrav. You cna find more information in the documentation for vibrav.util.io.open_txt.

from vibrav.util.io import open_txt

Reading all of the oscillator files that are generated.

iso_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-0.txt'),
                     rearrange=False)
iso_oscil.head()

	nrow	ncol	oscil	energy	freqdx	sign
0	0	1	2.931277e-11	0.002609	0	minus
1	0	2	1.970756e-11	0.002609	0	minus
2	0	3	3.586465e-11	0.002609	0	minus
3	0	4	1.083229e-10	0.004795	0	minus
4	1	4	2.710447e-11	0.001807	0	minus

x_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-1.txt'),
                   rearrange=False)
x_oscil.head()

	nrow	ncol	oscil	energy	freqdx	sign
0	0	1	5.378963e-12	0.002609	0	minus
1	0	2	1.240604e-11	0.002609	0	minus
2	0	3	4.269799e-13	0.002609	0	minus
3	0	4	9.580531e-11	0.004795	0	minus
4	1	4	3.322260e-11	0.001807	0	minus

y_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-2.txt'),
                   rearrange=False)
y_oscil.head()

	nrow	ncol	oscil	energy	freqdx	sign
0	0	1	6.395906e-12	0.002609	0	minus
1	0	2	9.883249e-15	0.002609	0	minus
2	0	3	6.191645e-11	0.002609	0	minus
3	0	4	7.761667e-11	0.004795	0	minus
4	1	4	1.647976e-12	0.001807	0	minus

z_oscil = open_txt(os.path.join('vibronic-outputs', 'oscillators-3.txt'),
                   rearrange=False)
z_oscil.head()

	nrow	ncol	oscil	energy	freqdx	sign
0	0	1	7.616344e-11	0.002609	0	minus
1	0	2	4.670676e-11	0.002609	0	minus
2	0	3	4.525051e-11	0.002609	0	minus
3	0	4	1.515468e-10	0.004795	0	minus
4	1	4	4.644282e-11	0.001807	0	minus

Cleaning up for the purpose of the example

os.chdir(parent)
shutil.rmtree('molcas-ucl6-2minus-vibronic-coupling')