Quick Start

System : H\({}_2\)CO

1. Install and Satisfy requirements

2. Import modules

Code

import mendeleev
import numpy as np
import scipy

Code

try:
    import bpy
except ImportError:
    print('You cannot use Blender in Jupyter')
    print('You must select jupyter kernel as blender when you use blender')
    import matplotlib
try:
    import ase
    import ase.db
except ImportError:
    print('You cannot use ASE')

Code

try:
    import pyvib
except ImportError:
    print('Try execute in `src` directory or set $PYTHONPATH to pyvib')

3. Prepare Geometry

Define - atom element - coordinate

For example,

geom = [['C', (0.0, 0.0, 0.0)],
        ['H', (1,0, 1.0, 1.0)]]

Here, I use precalculated ase database.

Code

geom, _ = pyvib.read_fchk_g16('./sample/ch2o.fchk')

Code

geom

[['O', (1.1408368, -0.000491967682, -0.000123437584)],
 ['C', (-1.13435014, 8.64849597e-05, -9.5923766e-06)],
 ['H', (-2.27091334, 0.465150475, -1.72441985)],
 ['H', (-2.27091603, -0.464367047, 1.72455288)]]

4. Prepare mass-weighted hessian

Define - mass-weighted hessian (unit is a.u.)

Mass-weighted hessian is matrix \(M_{ij}=\frac{\partial^2 E}{\partial\sqrt{m_i}x_i\partial\sqrt{m_j}x_j}\) of the second derivative of energy \(E\) in terms of mass-weighted coordinates \(\sqrt{m_i}x_i\). You can also obtain mass-weighted hessian from harmonic frequency and displacement vectors.

You can use pyvib.read_fchk_g16() for Gaussian16 or pyvib.read_minfo() for SINDO

Code

_, mw_hess = pyvib.read_fchk_g16('./sample/ch2o.fchk')

Code

mw_hess.shape

(12, 12)

5. Set PyViblocalizer

Set geometry and (hessian or (displacement vector and frequency))
Input units can be specified in options, such as unit_mass='AMU'

In large system, it may takes a few minutes

Code

vib = pyvib.Vibration(geom, mw_hess=np.array(mw_hess), 
                      unit_omega='hartree', unit_mass='a.u.')

/mnt/c/Users/hinom/GitHub/PyVibLocalizer/src/local_cls.py:147: RuntimeWarning: invalid value encountered in sqrt
  np.sqrt(np.diag(self.unitary.T@self.mw_hess@self.unitary).tolist()))

6.1. Visualize in Tkinter+matplotlib

6.2. Visualize in Blender

Code

try:
    import bpy
    vib.visualize(blender=True)
except ImportError:
    '''Not recommended'''
    vib.visualize(blender=False)

7. Localization

Returns displacement vecotors in a.u. and frequency in a.u.

7.1 Group Localization

Code

vib = pyvib.Vibration(geom, mw_hess=np.array(mw_hess), 
                      unit_omega='hartree', unit_mass='a.u.')
disp, freq = vib.group_localize(domain=[[0,1],[2,3]], mw_hess=np.array(mw_hess), 
                                   unit_omega='hartree', unit_mass='a.u.')

Code

try:
    import bpy
    vib.visualize(blender=True)
except ImportError:
    vib.visualize(blender=False)

7.2. Metric Localization

Code

vib = pyvib.Vibration(geom, mw_hess=np.array(mw_hess), 
                      unit_omega='hartree', unit_mass='a.u.')
disp, freq = vib.localize(option='Boys', window= 500)

initial zeta =  0.005301308461217721
0 zeta =  0.006251004066510544 delta = 0.006251004066510544
1 zeta =  0.006463777153999226 delta = 0.0002127730874886824
2 zeta =  0.006469695881229195 delta = 5.918727229968675e-06
Boys metric =  0.006469695881229195

Code

try:
    import bpy
    vib.visualize(blender=True)
except ImportError:
    vib.visualize(blender=False)