-
Notifications
You must be signed in to change notification settings - Fork 6
Tutorial
Peter Spackman edited this page Nov 14, 2023
·
4 revisions
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
To start, you'll need an input file that describes the geometry of the molecular crystal for which you want to predict the growth energies of the crystal. In this example we'll using 'urea.cif'. Ensure this file is in the working directory and contains a valid crystal structure in CIF format.
The contents of urea.cif are below:
data_urea
_symmetry_cell_setting tetragonal
_symmetry_space_group_name_H-M 'P -4 21 m'
_symmetry_Int_Tables_number 113
_space_group_name_Hall 'P -4 2ab'
loop_
_symmetry_equiv_pos_site_id
_symmetry_equiv_pos_as_xyz
1 x,y,z
2 -y,x,-z
3 -x,-y,z
4 y,-x,-z
5 1/2-x,1/2+y,-z
6 1/2+y,1/2+x,z
7 1/2+x,1/2-y,-z
8 1/2-y,1/2-x,z
_cell_length_a 5.662
_cell_length_b 5.662
_cell_length_c 4.716
_cell_angle_alpha 90
_cell_angle_beta 90
_cell_angle_gamma 90
_cell_volume 151.187
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
C1 C 0.0000 0.5000 0.3300
H1 H 0.2430 0.7430 0.2810
H2 H 0.1420 0.6420 0.0280
N1 N 0.1433 0.6433 0.1847
O1 O 0.0000 0.5000 0.5980
The overall steps that are made in the program are:
- Calculate the gas-phase and solvated wavefunctions for each symmetry unique molecule in the crystal.
- Calculate all the symmetry unique dimer interaction energies, using these to converge the lattice/cohesive energy of the crystal.
- Partition the solvation free energy up in order to obtain a pairwise approximation of the crystallisation free energies for each dimer.
Let's get this started!
An example of how to run this protocol is below, where we limit our lattice energy convergence to molecules within 12 angstroms, and our nearest neighbours will be the default 3.8 angstrom radius. We'll also provide the --threads flag to speed up the program with multiple threads for parallelisation.
occ cg urea.cif --radius=12.0 --threads=8The program output will look something like this:
found 1 symmetry unique molecules:
index name
0 CH4N2O
Molecule (0)
sym x y z
C 0.000000 5.349815 2.940943
H -2.600010 2.749805 2.504258
H 2.600010 7.949825 2.504258
H -1.519347 3.830467 0.249535
H 1.519347 6.869162 0.249535
N -1.533257 3.816558 1.646037
N 1.533257 6.883072 1.646037
O 0.000000 5.349815 5.329345
=============================== Input ================================
Method string b3lyp
Basis name 6-31g**
Shell kind Cartesian
Net charge 0
Multiplicity 1
Geometry '' (au) ------------------------------------------------------
C 0.000000 5.349815 2.940943
H -2.600010 2.749805 2.504258
H 2.600010 7.949825 2.504258
H -1.519347 3.830467 0.249535
H 1.519347 6.869162 0.249535
N -1.533257 3.816558 1.646037
N 1.533257 6.883072 1.646037
O 0.000000 5.349815 5.329345
Inertia tensor (x 10e-46 kg m^2) --------------------------------------
11.022624 -3.912356 -0.000000
-3.912356 11.022624 -0.000000
-0.000000 -0.000000 7.824711
Principal moments of inertia ------------------------------------------
7.110268 7.824711 14.934979
Rotational constants (GHz) --------------------------------------------
11.802673 10.725017 5.619035
Gas-phase properties (at 298.15 K) ------------------------------------
Rotational free energy -24.343770 kJ/mol
Translational free energy -39.009246 kJ/mol
Loaded basis set: cartesian
Number of shells: 36
Number of basis functions: 80
Maximum angular momentum: 2
Functionals:
hyb_gga_xc_b3lyp
0.2 x HF exchange
finished calculating atom grids (51080 points)
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 -224.842295975525 5.79602e-01 1.41842e-01 7.12e-02
2 -224.330620409683 1.44854e-03 3.35193e-01 7.45e-02
3 -223.631479468365 1.98317e-03 4.62802e-01 6.55e-02
4 -225.011799231770 3.90013e-03 2.29495e-02 6.86e-02
5 -225.014985179797 9.00189e-06 1.07745e-02 6.87e-02
6 -225.015996960961 2.85878e-06 1.22960e-03 6.78e-02
7 -225.016017204384 5.71976e-08 2.89652e-04 6.86e-02
8 -225.016018403883 3.38917e-09 8.91529e-05 6.82e-02
9 -225.016018451587 1.34787e-10 1.06870e-05 7.85e-02
10 -225.016018452506 2.59707e-12 4.09214e-06 6.85e-02
restricted spinorbital SCF energy converged after 0.70018 seconds
Component Energy (Hartree)
nuclear ---------------------------------------------------------------
nuclear.repulsion 128.904845051869
nuclear.point_charge 0.000000000000
nuclear.total 128.904845051869
electronic ------------------------------------------------------------
electronic.kinetic 225.098785629196
electronic.nuclear -785.841480821902
electronic.1e -560.742695192706
electronic -353.920863504375
electronic.2e 206.821831688330
electronic.dft_xc -24.087471562040
total -225.016018452506
Gas phase wavefunctions took 0.736804 seconds
Functionals:
hyb_gga_xc_b3lyp
0.2 x HF exchange
finished calculating atom grids (51080 points)
Using SMD solvent 'water'
Parameters:
Dielectric 78.3550
Intrinsic coulomb radii (4 atom types)
C 3.495993
H 2.267671
N 3.571582
O 2.872384
CDS radii
C 3.968425
H 3.023562
N 3.684966
O 3.628274
solvent surface:
total surface area (coulomb) = 82.764 Angstroms^2, 339 finite elements
total surface area (cds) = 105.247 Angstroms^2, 389 finite elements
CDS energy: 2.6461
CDS energy (molecular): 0.0000
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 -224.888431039731 5.78644e-01 1.24435e-01 1.01e-01
2 -224.684339381869 3.57240e-05 2.38763e-01 9.27e-02
3 -224.343488029878 1.61803e-03 3.36156e-01 9.49e-02
4 -225.037871866611 2.43071e-03 7.00708e-03 9.36e-02
5 -225.038108818340 9.70111e-06 1.41884e-03 9.67e-02
6 -225.038119628071 5.00000e-06 4.58227e-04 9.53e-02
7 -225.038122246603 3.02315e-06 1.12642e-04 9.49e-02
8 -225.038122394844 3.19109e-09 4.92348e-05 9.27e-02
9 -225.038122406832 1.00528e-07 2.07232e-06 9.34e-02
restricted spinorbital SCF energy converged after 0.85499 seconds
Component Energy (Hartree)
nuclear ---------------------------------------------------------------
nuclear.repulsion 128.904845051869
electronic ------------------------------------------------------------
electronic.kinetic 225.090705745791
electronic.nuclear -785.964897521006
electronic.1e -561.200339965710
electronic -354.240988763276
electronic.2e 206.959351202434
electronic.dft_xc -24.094280568910
solvation -------------------------------------------------------------
solvation.electronic -0.326148190494
solvation.nuclear 0.261460832693
solvation.surface 0.032343678900
solvation.CDS 0.004216792981
total -225.038122406832
SCF difference (au) -0.022
Solution phase wavefunctions took 0.955395 seconds
Computing monomer energies for gas phase
Computing monomer energies for solution phase
Computing monomer energies took 0.098070 seconds
Computing crystal interactions using gas wavefunctions
Found 31 symmetry unique dimers within max radius 12.000
Lattice convergence settings:
Start radius 3.8000 angstrom
Max radius 12.0000 angstrom
Radius increment 3.8000 angstrom
Energy tolerance 1.0000 kJ/mol
0 (0[x,y,z] - 0[y,-x,-z + (-1,0,0)]), Rc = 4.34
Took 0.595 seconds
1 (0[x,y,z] - 0[x,y,z + (0,0,-1)]), Rc = 4.72
Took 0.579 seconds
2 (0[x,y,z] - 0[y,-x,-z + (-1,0,-1)]), Rc = 5.02
Took 0.551 seconds
Finished calculating 3 unique dimer interaction energies
Cycle 1 lattice energy: -113.2237611912594
Total polarization term: -23.957
Total coulomb term: -85.342
Total dispersion term: -35.628
Total repulsion term: 90.610
Total exchange term: -48.910
Wolf correction: 0.000
3 (0[x,y,z] - 0[x,y,z + (-1,0,0)]), Rc = 5.66
Took 0.456 seconds
4 (0[x,y,z] - 0[x,y,z + (-1,-1,0)]), Rc = 8.01
Took 0.338 seconds
5 (0[x,y,z] - 0[x,y,z + (-1,0,-1)]), Rc = 7.37
Took 0.353 seconds
6 (0[x,y,z] - 0[y,-x,-z + (-1,0,1)]), Rc = 7.55
Took 0.328 seconds
7 (0[x,y,z] - 0[y,-x,-z + (-1,0,-2)]), Rc = 8.72
Took 0.324 seconds
8 (0[x,y,z] - 0[x,y,z + (-1,-1,-1)]), Rc = 9.29
Took 0.319 seconds
9 (0[x,y,z] - 0[y,-x,-z + (-2,0,0)]), Rc = 9.11
Took 0.311 seconds
10 (0[x,y,z] - 0[x,y,z + (0,0,-2)]), Rc = 9.43
Took 0.313 seconds
11 (0[x,y,z] - 0[y,-x,-z + (-2,0,-1)]), Rc = 9.45
Took 0.319 seconds
Finished calculating 9 unique dimer interaction energies
Cycle 2 lattice energy: -109.82984010715977
Total polarization term: -26.066
Total coulomb term: -75.796
Total dispersion term: -41.587
Total repulsion term: 90.736
Total exchange term: -48.954
Wolf correction: 0.000
12 (0[x,y,z] - 0[x,y,z + (-1,1,0)]), Rc = 8.01
Took 0.319 seconds
13 (0[x,y,z] - 0[x,y,z + (-1,1,-1)]), Rc = 9.29
Took 0.322 seconds
14 (0[x,y,z] - 0[x,y,z + (-1,0,-2)]), Rc = 11.00
Took 0.307 seconds
15 (0[x,y,z] - 0[y,-x,-z + (-2,0,-2)]), Rc = 11.84
Took 0.309 seconds
16 (0[x,y,z] - 0[y,-x,-z + (-2,0,1)]), Rc = 11.00
Took 0.318 seconds
17 (0[x,y,z] - 0[x,y,z + (-2,0,0)]), Rc = 11.32
Took 0.318 seconds
18 (0[x,y,z] - 0[x,y,z + (-2,-1,0)]), Rc = 12.66
Took 0.309 seconds
19 (0[x,y,z] - 0[y,-x,-z + (-1,0,2)]), Rc = 11.81
Took 0.309 seconds
20 (0[x,y,z] - 0[x,y,z + (-1,-1,-2)]), Rc = 12.37
Took 0.305 seconds
21 (0[x,y,z] - 0[x,y,z + (-2,0,-1)]), Rc = 12.27
Took 0.310 seconds
22 (0[x,y,z] - 0[y,-x,-z + (-2,-1,0)]), Rc = 12.13
Took 0.308 seconds
23 (0[x,y,z] - 0[y,-x,-z + (-1,0,-3)]), Rc = 13.09
Took 0.320 seconds
24 (0[x,y,z] - 0[x,y,z + (-2,-1,-1)]), Rc = 13.51
Took 0.315 seconds
25 (0[x,y,z] - 0[y,-x,-z + (-2,-1,-1)]), Rc = 12.39
Took 0.308 seconds
26 (0[x,y,z] - 0[x,y,z + (-1,1,-2)]), Rc = 12.37
Took 0.315 seconds
27 (0[x,y,z] - 0[y,-x,-z + (-2,-1,1)]), Rc = 13.61
Took 0.318 seconds
Finished calculating 16 unique dimer interaction energies
Cycle 3 lattice energy: -106.29780794770092
Total polarization term: -26.410
Total coulomb term: -71.618
Total dispersion term: -42.311
Total repulsion term: 90.736
Total exchange term: -48.954
Wolf correction: 0.000
28 (0[x,y,z] - 0[x,y,z + (0,0,-3)]), Rc = 14.15
Took 0.308 seconds
29 (0[x,y,z] - 0[y,-x,-z + (-2,-1,-2)]), Rc = 14.30
Took 0.314 seconds
30 (0[x,y,z] - 0[y,-x,-z + (-3,0,0)]), Rc = 14.53
Took 0.316 seconds
Finished calculating 3 unique dimer interaction energies
Cycle 4 lattice energy: -108.77668469615521
Total polarization term: -26.423
Total coulomb term: -73.928
Total dispersion term: -42.341
Total repulsion term: 90.736
Total exchange term: -48.954
Wolf correction: 0.000
Finished calculating 0 unique dimer interaction energies
Cycle 5 lattice energy: -108.77668469615521
Total polarization term: -26.423
Total coulomb term: -73.928
Total dispersion term: -42.341
Total repulsion term: 90.736
Total exchange term: -48.954
Wolf correction: 0.000
Neighbors for asymmetric molecule urea_0_water
nn Rn Rc Symop E_crys ES_AB ES_BA E_S E_nn E_int
========================================================================================
| 2.14 4.70 y,-x,1-z -31.35 -6.97 -2.40 -11.65 0.71 18.98
| 2.14 4.70 -1+y,-x,1-z -31.35 -2.40 -6.97 -7.09 1.77 22.49
| 2.14 4.70 y,1-x,1-z -31.35 -2.40 -6.97 -7.09 1.77 22.49
| 2.14 4.70 -1+y,1-x,1-z -31.35 -6.97 -2.40 -11.65 0.71 18.98
| 2.32 4.72 x,y,-1+z -40.87 -5.32 -26.28 -21.13 -0.21 19.95
| 2.32 4.72 x,y,1+z -40.87 -26.28 -5.32 -42.09 0.17 -1.40
| 2.77 4.60 -1+y,-x,-z -4.83 -2.59 -1.25 -4.51 1.53 -1.21
| 2.77 4.60 y,-x,-z -4.83 -1.25 -2.59 -3.18 0.46 1.20
| 2.77 4.60 -1+y,1-x,-z -4.83 -1.25 -2.59 -3.18 0.46 1.20
| 2.77 4.60 y,1-x,-z -4.83 -2.59 -1.25 -4.51 1.53 -1.21
Free energy estimates at T = 298 K, P = 1 atm., units: kJ/mol
-------------------------------------------------------
lattice energy (crystal) -108.777 (E_lat)
rotational free energy (molecule) -22.612 (E_rot)
translational free energy (molecule) -38.987 (E_trans)
solvation free energy (molecule) -50.126 (E_solv)
dH sublimation 103.821
dS sublimation -61.599
dG sublimation 42.223
dG solution -7.904
equilibrium_constant 2.43e+01
log S 1.385
solubility (g/L) 1.46e+03
Total E_int 101.486
Writing crystalgrower structure file to 'urea_cg.txt'
Writing crystalgrower net file to 'urea_water_net.txt'
Net dipole for molecule shell 0 = (-0.000 0.000 -18.229)
Energy = -252.430712 kJ/mol (-0.841 per molecule)
Wall clock time by category (s)
4-centre 2-electron integrals 83.154992
3-centre 2-electron integrals 1.712802
file input/output 0.016022
Initial guess 0.028778
MO update 0.018232
DIIS extrapolation 0.000552
DFT XC total 0.445319
GTO evaluation total 0.477069
Fock build 11.320893
Solvation 0.204454
Global (total time) 12.543910
A job well done
This is a lot of information, but the key information is all toward the end. First and most importantly, there's a table of the interactions for each molecule:
Neighbors for asymmetric molecule urea_0_water
nn Rn Rc Symop E_crys ES_AB ES_BA E_S E_nn E_int
========================================================================================
| 2.14 4.70 y,-x,1-z -31.35 -6.97 -2.40 -11.65 0.71 18.98
| 2.14 4.70 -1+y,-x,1-z -31.35 -2.40 -6.97 -7.09 1.77 22.49
| 2.14 4.70 y,1-x,1-z -31.35 -2.40 -6.97 -7.09 1.77 22.49
| 2.14 4.70 -1+y,1-x,1-z -31.35 -6.97 -2.40 -11.65 0.71 18.98
| 2.32 4.72 x,y,-1+z -40.87 -5.32 -26.28 -21.13 -0.21 19.95
| 2.32 4.72 x,y,1+z -40.87 -26.28 -5.32 -42.09 0.17 -1.40
| 2.77 4.60 -1+y,-x,-z -4.83 -2.59 -1.25 -4.51 1.53 -1.21
| 2.77 4.60 y,-x,-z -4.83 -1.25 -2.59 -3.18 0.46 1.20
| 2.77 4.60 -1+y,1-x,-z -4.83 -1.25 -2.59 -3.18 0.46 1.20
| 2.77 4.60 y,1-x,-z -4.83 -2.59 -1.25 -4.51 1.53 -1.21
These are the crystal growth energies, in the last column, with various contributions to them:
-
E_cryscorresponds to the dimer interaction energy (the same quantity calculated by the Total of theocc paircommand) -
E_Sis the portion of the solvation free energy attributed to this interaction -
E_nnis the contribution from theE_crysfor all the interactions that aren't in the nearest neighbors. If the calculation was limited to only nearest neighbours, this term would be 0. See the paper for more details.
Then we have the free energy estimates:
Free energy estimates at T = 298 K, P = 1 atm., units: kJ/mol
-------------------------------------------------------
lattice energy (crystal) -108.777 (E_lat)
rotational free energy (molecule) -22.612 (E_rot)
translational free energy (molecule) -38.987 (E_trans)
solvation free energy (molecule) -50.126 (E_solv)
dH sublimation 103.821
dS sublimation -61.599
dG sublimation 42.223
dG solution -7.904
equilibrium_constant 2.43e+01
log S 1.385
solubility (g/L) 1.46e+03
Total E_int 101.486
Highlights here include a predicted solubility, both in terms of the log S value and in g/L. Urea is highly soluble, but our estimate here isn't too far off the true value. It should be noted that the accuracy of the solubility prediction is largely orthogonal to prediction of the relative rates of growth for different interactions.
A lot of files are written out in this protocol. Many of which are useful for debugging or further calculations.