Example 8 Parallel Electronic Structure Calculation in DVR Mesh PES

NOTES

Electronic structure calculation requires ASE (Atomic Simulation Environment)

Documentation is here

1. Import modules

[1]:

import os

import numpy as np
from ase import Atoms
from ase.calculators.orca import ORCA, OrcaProfile
from ase.units import Bohr, Hartree
from ase.vibrations import Vibrations
from discvar import HarmonicOscillator

from pytdscf import units
from pytdscf.ase_handler import DVR_Mesh
from pytdscf.util.hess_util import get_displce_from_hess

2. Get Initial Coordinate

You can set ase.Atoms mannually. (Arleady optimized)

[2]:

h2co = Atoms(
    "OCHH",
    [
        (0.000000, 0.000000, 0.667507),
        (0.000000, 0.000000, -0.531172),
        (0.000000, 0.938955, -1.119091),
        (0.000000, -0.938955, -1.119091),
    ],
)

3. Select Calculator for opt/hessian

[3]:

nprocs = 4
MyOrcaProfile = OrcaProfile(command=f"{os.environ['ORCA_DIR']}/orca")
calculator_engrad = ORCA(
    profile=MyOrcaProfile,
    directory="calc",
    task="gradient",
    orcasimpleinput="B3LYP def2-TZVP TIGHTSCF Engrad",  # <- to calculate force, Engrad is required.
    orcablocks=f"%scf maxiter 300 end \n %pal nprocs {nprocs} end",
)
h2co.calc = calculator_engrad
h2co.get_forces()

[3]:

array([[ 4.31328299e-07, -2.77679162e-09,  3.50284979e-02],
       [ 2.40706696e-07, -6.77280046e-07, -3.09252283e-02],
       [ 4.06234330e-08,  7.40326482e-04, -1.32232997e-03],
       [ 5.48673456e-08, -7.40052711e-04, -1.32248352e-03]])

4. Run optimization (if needed)

`Berny optimizer were deleted <https://gitlab.com/ase/ase/-/merge_requests/3145>`__ (crying).

[4]:

# BFGS(h2co).run(fmax=0.001)

5. Prepare displacement vectors (normal mode) and frequencies

[5]:

vib = Vibrations(h2co)
vib.run()

hess = vib.get_vibrations().get_hessian_2d()
_freq = vib.get_frequencies()
print(hess, _freq)

vib.summary()
vib.clean()

[[ 1.86153277e+00  2.99790651e-04  5.04179227e-04 -5.57130462e+00
   9.22027231e-04 -6.37208114e-04  1.86644118e+00  9.86876885e-04
   5.61708950e-05  1.86726078e+00 -1.33864496e-05  7.33947164e-05]
 [ 2.99790651e-04  7.14653488e+00 -1.21917736e-03 -3.67693491e-05
  -1.02608808e+01  1.32207163e-03 -1.35799251e-04  1.56815736e+00
   3.56845373e+00  4.88697328e-04  1.56894987e+00 -3.56946525e+00]
 [ 5.04179227e-04 -1.21917736e-03  8.56946629e+01 -9.85218523e-04
   2.00084806e-03 -7.73154928e+01  1.17455715e-04  1.85300888e+00
  -4.19815850e+00 -4.23585421e-04 -1.85354833e+00 -4.19853771e+00]
 [-5.57130462e+00 -3.67693491e-05 -9.85218523e-04  1.68949096e+01
  -9.77686476e-04  8.06353291e-04 -5.67443819e+00 -2.23968813e-04
   9.72147034e-05 -5.67449069e+00  2.29287141e-04  8.88727586e-05]
 [ 9.22027231e-04 -1.02608808e+01  2.00084806e-03 -9.77686476e-04
   5.56456901e+01 -4.77772710e-04  6.37633632e-06 -2.26971598e+01
   7.65508379e+00 -2.08732455e-04 -2.26992016e+01 -7.65359793e+00]
 [-6.37208114e-04  1.32207163e-03 -7.73154928e+01  8.06353291e-04
  -4.77772710e-04  9.58307617e+01 -3.16496395e-04  7.19638854e+00
  -9.24112206e+00  1.99537032e-03 -7.19826156e+00 -9.24207078e+00]
 [ 1.86644118e+00 -1.35799251e-04  1.17455715e-04 -5.67443819e+00
   6.37633632e-06 -3.16496395e-04  1.90673396e+00 -1.59706399e-03
   9.13924398e-05  1.90453621e+00 -1.55207225e-04  9.82791402e-05]
 [ 9.86876885e-04  1.56815736e+00  1.85300888e+00 -2.23968813e-04
  -2.26971598e+01  7.19638854e+00 -1.59706399e-03  2.29844080e+01
  -1.01429642e+01  1.11453345e-03 -1.86153949e+00  1.08973492e+00]
 [ 5.61708950e-05  3.56845373e+00 -4.19815850e+00  9.72147034e-05
   7.65508379e+00 -9.24112206e+00  9.13924398e-05 -1.01429642e+01
   1.25083211e+01 -8.06056328e-04 -1.08852703e+00  9.23814966e-01]
 [ 1.86726078e+00  4.88697328e-04 -4.23585421e-04 -5.67449069e+00
  -2.08732455e-04  1.99537032e-03  1.90453621e+00  1.11453345e-03
  -8.06056328e-04  1.90079164e+00 -7.29364172e-04 -7.85721472e-04]
 [-1.33864496e-05  1.56894987e+00 -1.85354833e+00  2.29287141e-04
  -2.26992016e+01 -7.19826156e+00 -1.55207225e-04 -1.86153949e+00
  -1.08852703e+00 -7.29364172e-04  2.29869210e+01  1.01420041e+01]
 [ 7.33947164e-05 -3.56946525e+00 -4.19853771e+00  8.88727586e-05
  -7.65359793e+00 -9.24207078e+00  9.82791402e-05  1.08973492e+00
   9.23814966e-01 -7.85721472e-04  1.01420041e+01  1.25094645e+01]] [0.00000000e+00+35.59149439j 0.00000000e+00+25.48400563j
 7.12165794e-01 +0.j         9.56649195e+00 +0.j
 1.80922618e+01 +0.j         2.46433264e+01 +0.j
 1.20053422e+03 +0.j         1.26613208e+03 +0.j
 1.53380849e+03 +0.j         1.82104023e+03 +0.j
 2.87825020e+03 +0.j         2.93312380e+03 +0.j        ]
---------------------
  #    meV     cm^-1
---------------------
  0    4.4i     35.6i
  1    3.2i     25.5i
  2    0.1       0.7
  3    1.2       9.6
  4    2.2      18.1
  5    3.1      24.6
  6  148.8    1200.5
  7  157.0    1266.1
  8  190.2    1533.8
  9  225.8    1821.0
 10  356.9    2878.3
 11  363.7    2933.1
---------------------
Zero-point energy: 0.724 eV

[5]:

Absorb \(\sqrt{M}_i\) into \(\vec{R}_i\)

[6]:

sqrt_m = np.repeat(np.sqrt(h2co.get_masses()), 3)
mwhess = hess / sqrt_m[np.newaxis, :] / sqrt_m[:, np.newaxis]
mwhess *= (Bohr * Bohr / Hartree) * units.au_in_dalton
disp_vec, freq = get_displce_from_hess(mwhess, h2co.get_masses())
np.testing.assert_allclose(freq.real[6:], _freq.real[6:])
print(f"freq [cm-1] = {freq}")

freq [cm-1] = [           nan            nan 7.12165793e-01 9.56649195e+00
 1.80922618e+01 2.46433264e+01 1.20053422e+03 1.26613208e+03
 1.53380849e+03 1.82104023e+03 2.87825020e+03 2.93312380e+03]

/Users/hinom/GitHub/PyTDSCF/pytdscf/util/hess_util.py:259: RuntimeWarning: invalid value encountered in sqrt
  freq = np.sqrt(E) * units.au_in_cm1

6. Select DVR grid basis set

DVR grid must be consistent with Wave Function DVR basis

[7]:

ngrid = 7
jobname_pes = f"H2CO_{ngrid}grids"
dvr_info = [
    HarmonicOscillator(ngrid=ngrid, omega=f, units="cm1") for f in freq[6:]
]

7. Generate DVR grid Mesh coordiante

Here one constructs 2-mode representation grid-based PES

[8]:

mesh = DVR_Mesh(dvr_info, h2co, disp_vec[6:] * units.au_in_angstrom)
mesh_id = mesh.save_geoms(jobname_pes, overwrite=True, nMR=2)

START : Displacement Generation
DONE : Displacement Generation

Save geometry in DB of DOFs = (0,)[████████████████████████████████████████] 7/7

Save geometry in DB of DOFs = (1,)[████████████████████████████████████████] 7/7

Save geometry in DB of DOFs = (2,)[████████████████████████████████████████] 7/7

Save geometry in DB of DOFs = (3,)[████████████████████████████████████████] 7/7

Save geometry in DB of DOFs = (4,)[████████████████████████████████████████] 7/7

Save geometry in DB of DOFs = (5,)[████████████████████████████████████████] 7/7

Save geometry in DB of DOFs = (0, 1)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (0, 2)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (0, 3)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (0, 4)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (0, 5)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (1, 2)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (1, 3)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (1, 4)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (1, 5)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (2, 3)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (2, 4)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (2, 5)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (3, 4)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (3, 5)[████████████████████████████████████████] 49/49

Save geometry in DB of DOFs = (4, 5)[████████████████████████████████████████] 49/49

SAVE : DVR Mesh as H2CO_7grids.DVR_Mesh

8. Execute electronic structure calculation

You can check your CPU by cpuinfo command.

[9]:

mesh.execute_multiproc?

Signature:
mesh.execute_multiproc(
    calculator: ase.calculators.calculator.Calculator,
    max_workers: int | None = None,
    timeout: float = 60.0,
    jobname: str | None = None,
    reset_calc: bool = False,
    judge_func: Optional[Callable] = None,
)
Docstring:
Execute electronic structure calculation by multiprocessing

Args:
    calculator (Calculator) : calculator for each geomtry
    max_workers (Optional[int]) : maximum workers in multi-processing.
             Defaults to None. If None, use cpu_count - 1.
    timeout (float) : Timeout calculation in second. Defaults to 60.0
    jobname (Optional[str]) : jobname
    reset_calc (Optional[bool]) : set new calculator in any case.
            Defaults to False.
    judge_func (Optional[Callable[[Any],bool]]) : judge function whether re-calculation is needed.
            Defaults to None.
File:      ~/GitHub/PyTDSCF/pytdscf/ase_handler.py
Type:      method

[10]:

calculator = ORCA(
    profile=MyOrcaProfile,
    directory="calc",
    task="force",
    orcasimpleinput="B3LYP def2-TZVP TIGHTSCF Engrad",  # <- to calculate force, Engrad is required.
    orcablocks="%scf maxiter 300 end \n %pal nprocs 1 end",
)

[11]:

_ = mesh.execute_multiproc(calculator, max_workers=7, timeout=300.0)

unique jobs : 577
START : Electronic Structure Calculations
Connected: H2CO_7grids.db

[████████████████████████████████████████] 777/777

WAIT  : Remaining future task
DONE  : Electronic Structure Calculations
Your calculation completely finished!
DONE  : Shutdown process executor

9. Check Databse is exist.

[12]:

!ls *.db

H2CO_7grids.db