ORCA
Table of Contents

Resources

Online ORCA forum

Papers

Input

Geometry file

System Charge Multiplicity
neutral molecule 0 1
protonated PAHs C16H15+ 1 1
cation H2O+ 1 2
radical cation PAHs C16H14+ 1 2
- Linear molecules have to align with z-axis for vibrational analysis

spin multiplicity

$$\text{multiplicity} = \text{number of } \alpha \text{ electrons} - \text{number of } \beta \text{ electrons} + 1$$ - When #electrons is odd, then multiplicity must be even (2, 4, 6, ...). - When #electrons is even, then multiplicity must be old (1, 3, 5, ...).

SCF convergence

By default, TRAH is launched automatically if standard SCF (DISS/SOSCF) fails. It is possible to turn this option off by adding

! NOTRAH

or

%scf
    AutoTRAH false
end

Adding external electric field

Unit conversion for electric field is 1 a.u. = 51.4220674763 V/Å (1 * ase.units.Ha / ase.units.Bohr)

%scf
  efield X-field-strength-in-a.u. , Y-field-strength-in-a.u. , Z-field-strength-in-a.u.
end

RI

RIJCOSX

XC

meta-GGA

libxc

Using MN15 functional in ORCA 4.2.1 from libxc

! RIJCOSX def2-TZVP def2/J Opt nopop 
%scf maxiter 600 end 
%method
method dft 
exchange hyb_mgga_x_mn15 
correlation mgga_c_mn15 
end
%maxcore 1000 
%pal nprocs 16 end 
*xyz 0 1 
C -0.047131 0.664389 0.0
O -0.047131 -0.758551 0.0
H -1.092995 0.969785 0.0
H 0.878534 -1.048458 0.0
H 0.437145 1.080376 0.891772
H 0.437145 1.080376 -0.891772
*

Broken symmetry

CASSCF

Complete Active Space Self-Consistent Field (CASSCF) - CASSCF(N,M): N-electrons are active in M orbitals.

%casscf nel N
        norb M
        end

Metadynamics

ORCA 5.0.0 and 5.0.1 has a bug in interpreting the second Colvar. This can be solved by adding 3 more zeros in the end.

  Metadynamics Colvar 1 Scale 1.0_A Wall Lower 0.0 50.0 Upper 10.0 50.0 Range 0.0 10.0 100
  Metadynamics Colvar 2 Scale 1.0_A Wall Lower 0.0 50.0 Upper 10.0 50.0 Range 0.0 10.0 100 0 0 0

Alternatively, putting two Colvars in the same line also solves this bug

  Metadynamics Colvar 1 Scale 1.0_A Wall Lower 0.0 50.0 Upper 10.0 50.0 Colvar 2 Scale 1.0_A Wall Lower 0.0 50.0 Upper 10.0 50.0 Range 0.0 10.0 100 0.0 10.0 100

Optimization

For a big system (e.g. around 200 atoms), the local optimization in new redundant internal coordinates can be very slow. For example, updating positions could take several minutes. Diagonalizing the G-matrix is also slow. - %geom MaxIter 200 end sets the maximum optimization cycle to 200.

OptTS

Observation Parameter TS OPTIMIZATION
Energy Change TolE 5.0000e-06 Eh
Max. Gradient TolMAXG 3.0000e-04 Eh/bohr
RMS Gradient TolRMSG 1.0000e-04 Eh/bohr
Max. Displacement TolMAXD 4.0000e-03 bohr
RMS Displacement TolRMSD 2.0000e-03 bohr

In a TS optimization we are looking for a first order saddle point, and thus for a point on the PES where the curvature is negative in the direction of the TS mode (the TS mode is also called transition state vector, the only eigenvector of the Hessian at the TS geometry with a negative eigenvalue). As transition state finder we implemented the quasi-Newton like hessian mode following algorithm.This algorithm maximizes the energy with respect to one (usually the lowest) eigenmode and minimizes with respect to the remaining 3N - 7(6) eigenmodes of the Hessian.

OptTS uses the eigenvector following algorithm to find the nearest stationary point on PES. This can be local minimum or transition state.

NEB

PreOpt_Ends conflicts with NEB_Restart_XYZFile. You cannot use them at the same time.

zoomNEB

zoomNEB is similar to AutoNEB to some extents. Both methods insert images around the saddle point region on the fly. But zoomNEB inserts several images as many as initial images once identifying the saddle point region, while AutoNEB insert images one by one until a maximal pre-defined number.

Vibration

Analytical Frequency

Hessian and IR intensity

In orca.hess - $hessian has the 2nd energy derivatives w.r.t cartesian coordinates (not mass-weighted yet), with the unit of Ha / Bohr**2 and a shape of (3N, 3N) - $dipole_derivatives has the dipole derivative in cartesian coordinates, with the unit of e*Bohr / Bohr and a shape of (3N, 3)

import numpy as np
import ase.units as units
from ase import Atoms

def read_orca_hess(filename='orca.hess'):
    ''' 
        Read Hessian and dipole derevatives from the orca hess output file
    '''
    n_column = 5
    with open(filename, 'r') as f:
        for line in f:
            if '$hessian' in line:
                n_dof = int(next(f).strip().split(' ')[0])
                n_block = np.ceil(n_dof/n_column).astype('int')
                hessian = np.zeros((n_dof, n_dof))
                for idx_block in range(n_block):
                    next(f)
                    for i in range(n_dof):
                        content = next(f).strip().split(" ")[1:]
                        content = [float(item) for item in content if item != ""]
                        hessian[i, idx_block*n_column:(idx_block+1)*n_column] = content
                hessian *= units.Ha / (units.Bohr)**2
            if '$atoms' in line:
                n_atoms = int(next(f).strip().split(' ')[0])
                symbols = []
                masses = []
                positions = []
                for idx_atom in range(n_atoms):
                    content = next(f).strip().split(" ")
                    symbol = content[0]
                    symbols.append(symbol)

                    info = [float(item) for item in content[1:] if item != ""]
                    mass = info[0]
                    xyz = info[1:]
                    masses.append(mass)
                    positions.append(xyz)

                positions = np.array(positions) * units.Bohr
                atoms = Atoms(symbols, positions, masses=masses)
            if '$dipole_derivatives' in line:
                n_dof = int(next(f).strip().split(' ')[0])
                dipole_derivatives = np.zeros((n_dof, 3))
                for idx in range(n_dof):
                    content = next(f).strip().split(" ")
                    content = [float(item) for item in content if item != ""]
                    dipole_derivatives[idx] = content

    n_atoms = len(atoms)
    masses = atoms.get_masses()
    mass_weights = np.repeat(masses**-0.5, 3)
    omega2, vectors = np.linalg.eigh(mass_weights * hessian * mass_weights[:, np.newaxis])

    freq_unit_in_SI = np.sqrt(units._e/units._amu) * units.m
    freq_unit_in_eV = units._hbar * freq_unit_in_SI * units.J
    freq_in_invcm = omega2.astype(complex)**0.5 * freq_unit_in_eV / units.invcm
    freq_in_invcm = freq_in_invcm.real * np.sign(omega2)

    # ASE's version
    #unit_conversion = units._hbar * units.m / np.sqrt(units._e * units._amu)
    #energies = unit_conversion * omega2.astype(complex)**0.5
    #energies /= units.invcm

    modes = vectors.T * mass_weights

    intensities = np.einsum('ab,bg->ag', modes, dipole_derivatives) 
    intensities = np.sum(intensities**2, axis=1).real / units.Debye**2
    intensities *= 42.255
    for i in range(3*n_atoms):
        print(f'{freq_in_invcm[i].real:8.2f} cm^-1     {intensities[i]:6.2f} km/mol')

Thermochemistry

You can get results form a thermochemical analysis (based on ideal gas assumption) once freq is turned on. The most useful quantity for me is the Gibbs free energy G which includes contributions from electronic energies, vibrational frequencies and temperatures.

G = H - T*S

VPT2

VPT2 module is available again since ORCA 5.0.0. It is turned on by the keyword NEARIR. Since the evaluation of third and fourth derivatives are expensive, ORCA uses XTB2 to get high order derivatives by default. The default settings are writen under %FREQ:

%FREQ
XTBVPT2 True
DELQ 0.5
END

DELQ is the step during the numerical differentiation in dimensionless normal mode units. If one wants high order derivatives in the theory level as same as the single point calculation, XTBVPT2 False can be used. The total number of single point calculations in ORCA-VPT2 is (2N_mode+1)(N_mode+1) where N_mode = 3*N_atoms - 6 for non-linear molecules.

orca_pltvib

orca_pltvib is used for generating animations of vibrational modes. It takes orca.out or orca.hess as arguments. When analytical frequency is performed, orca.out is not compatible with orca_pltvib in ORCA 4.2.1.

Tips

One can use the following script to delete all .gbw, .densities and .tmp files in a directory and its subdirectories.

import os
rootdir = 'calculations/orca'

for subdir, dirs, files in os.walk(rootdir):
    for file in files:
        if file.endswith('.gbw') or file.endswith('.densities') or file.endswith('.tmp'):
            full_file_path = os.path.join(subdir, file)
            os.remove(full_file_path)
            #print(full_file_path)

Parallel run on multiple nodes

One has to create JOBNAME.nodes in the running directory to run on multiple nodes. echo $SLURM_NODELIST >> JOBNAME.nodes will create s31n[07-09] which cannot be recognized by ORCA. The following script gives the correct multiple nodes setting.

#SBATCH --nodelist=s31n[07-09]
#SBATCH --ntasks=36
#SBATCH --ntasks-per-node=12
#SBATCH --cpus-per-task=1

# change directory to the local scratch-directory, and run:
cd /scratch/$SLURM_JOB_ID
cat > orca.nodes <<!
s31n07 slots=12
s31n08 slots=12
s31n09 slots=12
!

Warnings

WARNING: For RIJCOSX frequencies GRIDX4 or larger grids are recommended

Errors

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!                        FATAL ERROR ENCOUNTERED               !!!
    !!!                        -----------------------               !!!
    !!!                          I/O OPERATION FAILED                !!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

I/O OPERATION FAILED is printed out when running out of disks. It is better to restart a job at a node with a larger disk.

[file orca_main/mainchk.cpp, line 1630]: ERROR: RI-MP2 needs an AuxC basis but none was defined!

An AuxC basis has to be defined explicitly in the input when using RI-MP2, such as RI-MP2 ANO-pVDZ cc-pVTZ/C TightSCF Engrad.

spin

[file orca_main/mainfcts.cpp, line 625]: Error : multiplicity (2) is even and number of electrons (204) is even -> impossible

Geometry

    ----------------------------------------------------------------------------
                                   ERROR !!!
       The optimization did not converge but reached the maximum 
       number of optimization cycles.
       As a subsequent Frequencies calculation has been requested
       ORCA will abort at this point of the run.
       Please restart the calculation with the lowest energy geometry and/or
       a larger maxiter for the geometry optimization.
    ----------------------------------------------------------------------------
      [file orca_main/run.cpp, line 10680]: ORCA finished with error return - aborting the run