-
Notifications
You must be signed in to change notification settings - Fork 6
Tutorial
To start, you'll need an input file that describes the geometry of the molecule you're studying. In this example we're using 'water.xyz'. Ensure this file is in the working directory and contains the atomic coordinates of the water molecule in XYZ format an example is below:
3
# comment line
O -0.7021961 -0.0560603 0.0099423
H -1.0221932 0.8467758 -0.0114887
H 0.2575211 0.0421215 0.0052190
With 'water.xyz' ready, you can run a simple SCF calculation with a minimal basis set. Here’s the basic command structure:
# Assuming occ is available in your $PATH
occ scf water.xyz rhf 3-21gThis command tells the program to perform a Hartree-Fock calculation using the 3-21G basis set on the 'water.xyz' geometry.
Execute the command in the terminal. If you're writing a script, you can save the command as a script file with a .sh extension and run it with sh your_script.sh.
The program should have an output very similar to below
Open
\
Comp
\
Chem
a library and program for quantum chemistry
copyright (C) 2020-2022 Peter Spackman
this version of occ makes use of the following third party libraries:
CLI11 command line argument parser
eigen3 Linear Algebra (v 3.4.0)
fmt String formatting (v 10.1.0)
gau2grid Gaussian basis function evaluation (v 2.0.7)
gemmi CIF parsing & structure refinement (v 0.6.3)
LBFGS++ L-BFGS implementation
libcint Electron integrals using GTOs (v CINT_VERSION)
libxc Density functional implementations (v 6.2.2)
nanoflann KDtree implementation
nlohmann::json JSON parser
phmap Fast hashmap implementation
pocketFFT Standalone implementation of fast Fourier transform
scnlib String parsing
spdlog Logging (v 1.12.0)
subprocess Calling external subprocesses
Parallelization: 1 OpenMP threads, 1 eigen threads
Spinorbital kind: Restricted
=============================== Input ================================
Method string rhf
Basis name 3-21g
Shell kind Cartesian
Net charge 0
Multiplicity 1
Geometry 'water.xyz' (au) ---------------------------------------------
O -1.326958 -0.105939 0.018788
H -1.931665 1.600174 -0.021710
H 0.486644 0.079598 0.009862
Inertia tensor (x 10e-46 kg m^2) --------------------------------------
0.128732 0.038580 -0.000546
0.038580 0.167540 0.003071
-0.000546 0.003071 0.296123
Principal moments of inertia ------------------------------------------
0.104928 0.191270 0.296197
Rotational constants (GHz) --------------------------------------------
799.788374 438.753506 283.325060
Gas-phase properties (at 298.15 K) ------------------------------------
Rotational free energy -9.658515 kJ/mol
Translational free energy -34.532039 kJ/mol
Loaded basis set: cartesian
Number of shells: 9
Number of basis functions: 13
Maximum angular momentum: 1
Including potential from 0 point charges
Computing core hamiltonian
Computing initial guess using SOAD in minimal basis
starting restricted scf iterations
# E (Hartrees) |dE|/E max|FDS-SDF| T (s)
1 -75.570107461651 6.47838e-04 5.31104e-02 1.01e-03
2 -75.584244261472 1.66824e-04 1.96361e-02 9.96e-04
3 -75.585159527555 1.08006e-05 7.07734e-03 9.67e-04
4 -75.585320851767 1.90371e-06 1.23457e-03 9.61e-04
5 -75.585325636630 5.64639e-08 1.10974e-04 9.70e-04
6 -75.585325703594 7.90205e-10 6.89991e-06 9.82e-04
restricted spinorbital SCF energy converged after 0.00589 seconds
Component Energy (Hartree)
nuclear ---------------------------------------------------------------
nuclear.repulsion 9.156713561444
nuclear.point_charge 0.000000000000
nuclear.total 9.156713561444
electronic ------------------------------------------------------------
electronic.kinetic 75.452812174289
electronic.nuclear -198.019164066723
electronic.1e -122.566351892434
electronic -84.742039265038
electronic.2e 37.824312627395
total -75.585325703594
wavefunction stored in water.owf.json
======================== Converged Properties ========================
Center of Mass -1.259322 -0.000102 0.016023
Molecular Multipole Moments (au) --------------------------------------
q -0.000000
Dx 0.515350 Dy 0.813557 Dz -0.021245
Qxx -3.382082 Qxy -0.458593 Qxz 0.007201
Qyy -3.831282 Qyz -0.029514 Qzz -5.075671
Oxxx 1.630097 Oxxy 0.292888 Oxxz -0.011309
Oxyy -0.627176 Oxyz 0.013666 Oxzz -0.038710
Oyyy 1.290941 Oyyz -0.032436 Oyzz -0.043947
Ozzz 0.003559
Hxxxx -13.034491 Hxxxy 0.788456 Hxxxz -0.022120
Hxxyy -4.954950 Hxxyz 0.004995 Hxxzz -4.723439
Hxyyy -0.426280 Hxyyz 0.023912 Hxyzz 0.248129
Hxzzz -0.011014 Hyyyy -12.637748 Hyyyz -0.017099
Hyyzz -4.457029 Hyzzz 0.053453 Hzzzz -11.128986
Mulliken charges (au) -------------------------------------------------
O -0.724463
H 0.363043
H 0.361419
Wall clock time by category (s)
4-centre 2-electron integrals 0.005724
file input/output 0.005336
Initial guess 0.003523
MO update 0.000128
DIIS extrapolation 0.000021
Fock build 0.005735
Global (total time) 0.014192
A job well done
Further, a listing of the contents of the directory should look something like this:
ls -lh
total 36K
-rw-rw-r-- 1 prs prs 29K Nov 9 08:55 water.owf.json
-rw-rw-r-- 1 prs prs 124 Nov 9 08:55 water.xyzThe water.owf.json file contains the wavefunction information.
There's a lot of information in the output, so let's break some of it down.
This section:
starting restricted scf iterations
# E (Hartrees) |dE|/E max|FDS-SDF| T (s)
1 -75.570107461651 6.47838e-04 5.31104e-02 1.01e-03
2 -75.584244261472 1.66824e-04 1.96361e-02 9.96e-04
3 -75.585159527555 1.08006e-05 7.07734e-03 9.67e-04
4 -75.585320851767 1.90371e-06 1.23457e-03 9.61e-04
5 -75.585325636630 5.64639e-08 1.10974e-04 9.70e-04
6 -75.585325703594 7.90205e-10 6.89991e-06 9.82e-04
Shows the energy at each iteration of the SCF, along with the (relative) change in energy, and another convergence criterion (the maximum change in the commutator). We also get a print out of the time.
This part of the program will typically be where 99% of the calculation is spent, so the time print out can be useful if we want an idea of how long it might take.
We also get a print out of some molecular properties, such as the multipole moments:
Molecular Multipole Moments (au) --------------------------------------
q -0.000000
Dx 0.515350 Dy 0.813557 Dz -0.021245
Qxx -3.382082 Qxy -0.458593 Qxz 0.007201
Qyy -3.831282 Qyz -0.029514 Qzz -5.075671
Oxxx 1.630097 Oxxy 0.292888 Oxxz -0.011309
Oxyy -0.627176 Oxyz 0.013666 Oxzz -0.038710
Oyyy 1.290941 Oyyz -0.032436 Oyzz -0.043947
Ozzz 0.003559
Hxxxx -13.034491 Hxxxy 0.788456 Hxxxz -0.022120
Hxxyy -4.954950 Hxxyz 0.004995 Hxxzz -4.723439
Hxyyy -0.426280 Hxyyz 0.023912 Hxyzz 0.248129
Hxzzz -0.011014 Hyyyy -12.637748 Hyyyz -0.017099
Hyyzz -4.457029 Hyzzz 0.053453 Hzzzz -11.128986
These are provided in atomic units.
After running the calculation, the program will generate output files and data that describe the results of the SCF process. The output will likely include the final energy, the number of iterations the SCF took, and other relevant information.
Let's run the calculation again using b3lyp and def2-svp as the basis set, then we can use the wavefunction in a pair calculation.
occ water.xyz b3lyp def2-svpYou should get a similar print out, though the calculation will take longer largely due to a bigger basis set.
Now that we have a wavefunction for water, we can calculate a dimer interaction energy using the pair subcommand.
First, let's prepare an input file:
[pair]
model = "ce-1p"
monomer_a = "water.owf.json"
monomer_b = "water.owf.json"
rotation_a = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
rotation_b = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
translation_a = [0.0, 0.0, 0.0]
translation_b = [3.0, 0.0, 0.0]The output should look something like the following:
Monomer A energies
E_coul 46.862375921777
E_ex -8.941339023333
E_nn 9.156713561444
E_en -198.548826042149
E_kin 75.513398158215
E_1e -123.035427883934
E_ecp 0.000000000000
Monomer B energies
E_coul 46.862375921777
E_ex -8.941339023333
E_nn 9.156713561444
E_en -198.548826042149
E_kin 75.513398158215
E_1e -123.035427883934
E_ecp 0.000000000000
Dimer
Component Energy (kJ/mol)
Coulomb -4.150470
Exchange -18.903424
Repulsion 34.865956
Polarization -2.547138
Dispersion -1.734868
__________________________________
Total 0.371578
Wall clock time by category (s)
4-centre 2-electron integrals 0.070938
3-centre 2-electron integrals 0.000898
file input/output 0.001619
GTO evaluation total 0.044342
Fock build 0.070944
XDM dispersion 0.138848
Global (total time) 0.500709
A job well done
The key energy quantity here is the Total. Which in this case indicates that the two water molecules are repelling each other.
However, the individual energy components can be quite useful in understanding why, for more details see this paper