External Methods¶
This tutorial demonstrates how to perform an ORCA calculation via OPI with external wrapper scripts to allow the use of methods not natively implemented in ORCA. As an example, we will run an optimization of the phenylalanine–serine dipeptide with PM7. Therefore, the respective wrapper must be downloaded (available via GitHub) and MOPAC has to be installed (also available via GitHub). Besides geometry optimization, also other ORCA functionalities like GOAT or a NEB-TS are available in combination with wrapper scripts.
Note
Besides MOPAC and PM7, you can use any program and method of your liking as long as you have the respective program installed and the wrapper available. If you create or use a custom wrapper script, please consider contributing on GitHub to make it accessible for everyone.
In this notebook we will:
Import the required python dependencies.
Define a working directory.
Set up the input structure.
Define the calculation settings.
Perform the ORCA calculation.
Check for proper termination.
Process the results.
Step 1: Import Dependencies¶
The first thing to do for using OPI is to import the modules needed. We additionally import also the modules for visualization like py3Dmol. For this, it might be necessary to install py3Dmol into your OPI venv (e.g., by activating the .venv and using pip install py3Dmol).
from pathlib import Path
import shutil
# > OPI imports for performing ORCA calculations and reading output
from opi.core import Calculator
from opi.output.core import Output
from opi.input.simple_keywords import Opt
from opi.input.structures.structure import Structure
from opi.input.blocks import BlockMethod
from opi.input.simple_keywords import ExternalTools
# > for visualization of molecules
import py3Dmol
Step 2: Working Directory¶
When working with OPI, it is always a good idea to define a subdirectory in which the actual ORCA calculations will take place. In this case, we set up a RUN directory:
# > Calculation is performed in `RUN`
working_dir = Path("RUN")
# > The `working_dir`is automatically (re-)created
shutil.rmtree(working_dir, ignore_errors=True)
working_dir.mkdir()
Step 3: Setup the Input Structure¶
We define the structure with Cartesian coordinates in angstrom and visualize it:
# > Define cartesian coordinates in angstrom as python string
xyz_data = """\
43
H -4.66498485894935 -0.74241492470540 0.69974381312263
C -4.17671963772691 -0.41946813411513 1.62085663387098
H -4.31556295019038 0.65835557933743 1.74617763893195
H -4.64735205939740 -0.91420097769670 2.47504134616210
C -2.69928997074851 -0.71449720892044 1.64944528147625
O -1.98975107774701 -0.42928398115826 2.62017129759124
N -2.17773611732287 -1.30327445426666 0.54436388078954
H -2.78004021432631 -1.52785357956148 -0.23421300168168
C -0.76472603182702 -1.64082193548350 0.46301524117174
H -0.48934972945319 -2.31114326509938 1.28555265397361
C -0.49680188103055 -2.34091857378744 -0.88594861326624
H -1.12591070939924 -3.23791564349719 -0.92078590384981
H -0.81348050791786 -1.67401763664883 -1.69607155624943
C 0.95062385323134 -2.71378129189350 -1.06213996849627
C 1.47509728895600 -3.83883958097587 -0.41893480660121
C 1.80206493224889 -1.92339666967511 -1.83795922104143
C 2.82283098211029 -4.16358701199868 -0.54351664024624
C 3.15157817659875 -2.24582496944767 -1.96455890072444
C 3.66615600696023 -3.36582314961534 -1.31551415939035
H 0.82015336645497 -4.46361266116421 0.18316223338911
H 1.40383614174953 -1.04827263842458 -2.34534914498881
H 3.21532334050428 -5.04205196812944 -0.04017998407723
H 3.80053282584475 -1.62233049211834 -2.57203001397934
H 4.71717343897748 -3.61922096939729 -1.41415620606549
C 0.13897619924998 -0.40963929818479 0.62323338935093
O 1.18982465311154 -0.47050762380759 1.27051960118465
N -0.26042356408273 0.69305102705598 -0.03880513029755
H -1.17884625655101 0.71224550095703 -0.46453347384447
C 0.46632507370965 1.95671966072207 0.00061241530592
H 1.53922104685905 1.73383868295454 0.04807577903327
C 0.17476690677966 2.75108272784102 -1.26783749196997
H 0.72597721347984 3.69989315677944 -1.22599656255061
H 0.52311351799454 2.18298120102531 -2.13543127758411
O -1.22607123349587 2.96938238237637 -1.44869581637241
H -1.49646216555062 3.50039758318850 -0.68132009170136
C 0.06279389695701 2.75700018989079 1.25788607468049
O -0.64234099967927 3.77463473172127 1.18777548552859
N 0.51212557634628 2.25231891234389 2.41467071805516
H 1.04291833952237 1.39006469799826 2.38265045795352
C 0.12919661680901 2.81402640326745 3.70411596033363
H 0.45784513235661 3.85472478200394 3.78674503246121
H -0.95777954492060 2.78312669629151 3.83756598938075
H 0.60317498350469 2.22685472401805 4.49159704123113\n
"""
# > Visualize the input structure
view = py3Dmol.view(width=400, height=400)
view.addModel(xyz_data, 'xyz')
view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})
view.zoomTo()
view.show()
# > Write the input structure to a file
with open(working_dir / "struc.xyz","w") as f:
f.write(xyz_data)
# > Read structure into object
structure = Structure.from_xyz(working_dir / "struc.xyz")
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Step 4: Define the Calculation Settings¶
Next, we have to define the settings used in the ORCA optimization. This can be done by defining a calculator and adding simple keywords to it:
def setup_calc(basename : str, working_dir: Path, structure: Structure, ncores: int = 1) -> Calculator:
# > Set up a Calculator object, the basename for the ORCA calculation is also set
calc = Calculator(basename=basename, working_dir=working_dir)
# > Assign structure to calculator
calc.structure = structure
# > List of simple keywords for ORCA
sk_list = [
ExternalTools.EXTOPT, # Set the extopt keyword for the ORCA input
Opt.OPT # > Perform an optimization
]
# > Use simple keywords in calculator
calc.input.add_simple_keywords(*sk_list)
# > Set the path to the wrapper script that should be used
# > Without adjusting the path this notebook will not work!
path_to_wrapper = "/home/full/path/to/orca-external-tools/mopac.sh"
# Additional arguments for the wrapper script
arguments = '"--method PM7"'
# Add the block list to calculator
calc.input.add_blocks(BlockMethod(ProgExt=path_to_wrapper,Ext_Params=arguments))
# > Define number of CPUs for the calcualtion
calc.input.ncores = ncores # > CPUs for this ORCA run (default: 1)
return calc
basename = "pm7_opt"
calc = setup_calc(basename=basename, working_dir=working_dir, structure=structure)
Note
When defining Ext_Params, both types of quotes are required in a format like ('"arguments"').
Step 5: Perform the Calculation¶
To perform the single-point calculation, we have to write the input in ORCA format first, and then run the job by using the run() function of the calculator class. The output of ORCA can be written to the calculator with its get_output() function.
def run_calc(calc: Calculator) -> Output:
# > Write the ORCA input file
calc.write_input()
# > Run the ORCA calculation
print("Running ORCA calculation ...", end="")
calc.run()
print(" Done")
# > Get the output object
output = calc.get_output()
return output
output = run_calc(calc)
Running ORCA calculation ... Done
Step 6: Check for Proper Termination¶
When performing automated computations, a check for proper termination should always be included. A simple check could look like:
def check_and_parse_output(output: Output):
# > Check for proper termination of ORCA
status = output.terminated_normally()
if not status:
# > ORCA did not terminate normally
raise RuntimeError(f"ORCA did not terminate normally. Please check RUN/{basename}.out")
else:
# > ORCA did terminate normally so we can parse the output
# > We must switch off the processing of the gbw output file
# as a gbw file is not written for extopt calculations
output.parse(read_gbw_json=False)
check_and_parse_output(output)
Step 7: Access the Results¶
After a successful calculation, we can access the resulting properties. They are generally available via the output object:
# > Get the PM7 energy of the final structure (last element in geometries)
Etot = output.results_properties.geometries[-1].single_point_data.finalenergy
# > Print total energy
print("Final energy in Hartree: ",Etot)
# > The final geometry after optimization is stored as the last element in geometries:
coords = output.results_properties.geometries[-1].geometry.coordinates.cartesians
# > Convert the coordinates from bohr to angstrom
bohr_to_angstrom = 0.529177
coords = [
(atom, x * bohr_to_angstrom, y * bohr_to_angstrom, z * bohr_to_angstrom)
for atom, x, y, z in coords
]
# > Make a simple string to be processed by py3Dmol
xyz_str = f"{len(coords)}\n\n"
for atom, x, y, z in coords:
xyz_str += f"{atom} {x:.6f} {y:.6f} {z:.6f}\n"
# Visualize with py3Dmol
view = py3Dmol.view(width=400, height=400)
view.addModel(xyz_str, 'xyz')
view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})
view.zoomTo()
view.show()
Final energy in Hartree: -0.279731653465
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Summary¶
In this notebook we demonstrated how to perform an ORCA optimization with the ORCA Python Interface (OPI) and external methods via wrapper scripts.