Tutorial 2: Calcium solvation with Orca¶
The goal of this tutorial is to run the pyMERK program on a set of 10 conformers of \(\ce{Ca_2(BH_4)_4(THF)_2}\) in order to identify and optimize the most relevant structures. All quantum chemistry calculations in this workflow are performed using Orca.
You can download the initial structures here: goat.finalensemble.xyz.
These geometries were generated from a GOAT calculation using rigid fragments.
Note
Running the full workflow may take several hours on a standard workstation. The reference calculation presented here was performed on a supercomputer. Precomputed output files are provided if you prefer to follow the tutorial without running the calculations.
Setting up & running¶
For this tutorial, we use the following input file: input.toml.
Apart from the choice of solvent (tetrahydrofuran instead of water), most parameters are left at their default values as described in Installation & usage.
This tutorial assumes that xTB and Orca are properly installed and available in your environment.
Before running the workflow, verify and update the paths to xTB ([paths] xtb) and Orca ([paths] orca):
[paths]
orca = '<PATH_TO_ORCA_BIN>/orca'
orca_nprocs = 32
xtb = 'xtb'
Tip
To enable MPI parallelization in Orca, you must provide the full path to the Orca executable.
See Parallel instructions for ORCA for details.
You can then adjust [paths] orca_nprocs to control the number of processes.
Run the workflow with:
pymerk_run -i input.toml goat.finalensemble.xyz -o final_ensemble.xyz
Once the workflow starts, progress can be monitored from the standard output.
If you prefer not to run the calculation, you can use the provided run.log.
Analyzing the output¶
Although this system is significantly larger than in the previous tutorial, the use of the r²SCAN-3c approach in Orca greatly accelerates the screening and optimization stages.
The input section corresponding to the optimization stage is:
[optimization]
prog = 'orca'
sm = 'smd'
alternate_solvent = 'tetrahydrofuran'
basis = 'def2-mTZVPP'
func = 'r2scan-3c'
This setup uses the r²SCAN-3c/def2-mTZVPP method in tetrahydrofuran, with solvation described via SMD.
The corresponding output is:
> Macrocycle 1
- Otimizing Conformer #1... E=-1928.77302771, grad=0.016784 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.77295768, grad=0.016667 [NOT CONVERGED]
- Otimizing Conformer #3... E=-1928.77330021, grad=0.018177 [NOT CONVERGED]
- Otimizing Conformer #4... E=-1928.77301709, grad=0.017722 [NOT CONVERGED]
- Otimizing Conformer #5... E=-1928.77229602, grad=0.017388 [NOT CONVERGED]
- Otimizing Conformer #6... E=-1928.77361177, grad=0.017566 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.77349481, grad=0.017240 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.77303233, grad=0.017279 [NOT CONVERGED]
- Otimizing Conformer #9... E=-1928.77382920, grad=0.017675 [NOT CONVERGED]
→ retained 9 conformer(s)
> Macrocycle 2
- Otimizing Conformer #1... E=-1928.77914650, grad=0.010644 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.77919673, grad=0.010655 [NOT CONVERGED]
- Otimizing Conformer #3... E=-1928.77999569, grad=0.011913 [NOT CONVERGED]
- Otimizing Conformer #4... E=-1928.77953865, grad=0.010643 [NOT CONVERGED]
- Otimizing Conformer #5... E=-1928.77937474, grad=0.015141 [NOT CONVERGED]
- Otimizing Conformer #6... E=-1928.77954065, grad=0.011201 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.77967820, grad=0.011955 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.77917046, grad=0.013855 [NOT CONVERGED]
- Otimizing Conformer #9... E=-1928.77916968, grad=0.009322 [NOT CONVERGED]
→ retained 9 conformer(s)
> Macrocycle 3
- Otimizing Conformer #1... E=-1928.78151517, grad=0.010718 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.78133010, grad=0.010217 [NOT CONVERGED]
- Otimizing Conformer #3... E=-1928.78188768, grad=0.012678 [NOT CONVERGED]
- Otimizing Conformer #4... E=-1928.78161302, grad=0.010248 [NOT CONVERGED]
- Otimizing Conformer #5... E=-1928.78257420, grad=0.011785 [NOT CONVERGED]
- Otimizing Conformer #6... E=-1928.78220200, grad=0.011244 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.78249045, grad=0.011445 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.78235279, grad=0.012536 [NOT CONVERGED]
- Otimizing Conformer #9... E=-1928.78129662, grad=0.010044 [NOT CONVERGED]
→ retained 9 conformer(s)
> Macrocycle 4
- Otimizing Conformer #1... E=-1928.78490641, grad=0.018756 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.78381092, grad=0.012692 [NOT CONVERGED]
- Otimizing Conformer #3... E=-1928.78627929, grad=0.012779 [NOT CONVERGED]
- Otimizing Conformer #4... E=-1928.78578860, grad=0.020505 [NOT CONVERGED]
- Otimizing Conformer #5... E=-1928.78839902, grad=0.018355 [NOT CONVERGED]
- Otimizing Conformer #6... E=-1928.78597848, grad=0.017054 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.78661962, grad=0.016879 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.78652716, grad=0.015147 [NOT CONVERGED]
- Otimizing Conformer #9... E=-1928.78458063, grad=0.012629 [NOT CONVERGED]
→ retained 9 conformer(s)
> Macrocycle 5
- Otimizing Conformer #1... E=-1928.79369502, grad=0.012539 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.79460004, grad=0.029011 [NOT CONVERGED]
- Otimizing Conformer #3... E=-1928.79406931, grad=0.012735 [NOT CONVERGED]
- Otimizing Conformer #4... E=-1928.79495595, grad=0.012822 [NOT CONVERGED]
- Otimizing Conformer #5... E=-1928.79560609, grad=0.008616 [NOT CONVERGED]
- Otimizing Conformer #6... E=-1928.79260713, grad=0.011855 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.79171977, grad=0.006051 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.79254052, grad=0.025737 [NOT CONVERGED]
- Otimizing Conformer #9... E=-1928.79014324, grad=0.007283 [NOT CONVERGED]
→ retained 8 conformer(s)
> Macrocycle 6
- Otimizing Conformer #1... E=-1928.79632813, grad=0.009871 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.79706226, grad=0.008321 [NOT CONVERGED]
- Otimizing Conformer #3... E=-1928.79554638, grad=0.005453 [NOT CONVERGED]
- Otimizing Conformer #4... E=-1928.79498081, grad=0.003733 [NOT CONVERGED]
- Otimizing Conformer #5... E=-1928.79596289, grad=0.005593 [NOT CONVERGED]
- Otimizing Conformer #6... E=-1928.79477773, grad=0.007267 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.79423398, grad=0.006896 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.79573457, grad=0.007073 [NOT CONVERGED]
→ retained 8 conformer(s)
> Macrocycle 7
- Otimizing Conformer #1... E=-1928.79837042, grad=0.006560 [NOT CONVERGED]
- Otimizing Conformer #2... E=-1928.79715134, grad=0.008322 [CONVERGED]
- Otimizing Conformer #3... E=-1928.79564309, grad=0.005453 [CONVERGED]
- Otimizing Conformer #4... E=-1928.79437346, grad=0.003732 [CONVERGED]
- Otimizing Conformer #5... E=-1928.79602255, grad=0.005598 [CONVERGED]
- Otimizing Conformer #6... E=-1928.79770014, grad=0.010852 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.79681750, grad=0.013492 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.79886170, grad=0.010926 [NOT CONVERGED]
→ retained 8 conformer(s)
> Macrocycle 8
- Otimizing Conformer #1... E=-1928.79841420, grad=0.006561 [CONVERGED]
- Otimizing Conformer #6... E=-1928.80003175, grad=0.006218 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.80041361, grad=0.005542 [NOT CONVERGED]
- Otimizing Conformer #8... E=-1928.80071451, grad=0.005523 [NOT CONVERGED]
→ retained 6 conformer(s)
> Macrocycle 9
- Otimizing Conformer #6... E=-1928.80123923, grad=0.004664 [NOT CONVERGED]
- Otimizing Conformer #7... E=-1928.80048334, grad=0.005540 [CONVERGED]
- Otimizing Conformer #8... E=-1928.80147675, grad=0.006496 [NOT CONVERGED]
→ retained 5 conformer(s)
> Macrocycle 10
- Otimizing Conformer #6... E=-1928.80133117, grad=0.004664 [CONVERGED]
- Otimizing Conformer #8... E=-1928.80136756, grad=0.003163 [NOT CONVERGED]
→ retained 5 conformer(s)
> Macrocycle 11
- Otimizing Conformer #8... E=-1928.80139099, grad=0.003163 [CONVERGED]
→ retained 5 conformer(s)
* Final values:
G* ΔG*
Conformer #1 -1928.79841420 0.00297679 *
Conformer #2 -1928.79715134 0.00423965 *
Conformer #3 -1928.79564309 0.00574790
Conformer #4 -1928.79437346 0.00701754
Conformer #5 -1928.79602255 0.00536844
Conformer #6 -1928.80133117 0.00005982 *
Conformer #7 -1928.80048334 0.00090765 *
Conformer #8 -1928.80139099 0.00000000 *
Conformer #9 -1928.79014324 0.01124775
* Done, retained 5 conformer(s)
A total of 11 macrocycles are required to reduce the set to the 5 lowest-energy conformers out of the 9 retained after the previous stages.
The selection effectively begins at macrocycle 5, when the gradient norms fall below [optimization] gradthr.
This illustrates the efficiency of the macrocycle procedure in discarding high-energy conformers early.
The refinement stage further reduces the ensemble using a higher level of theory (wb97x-d3/def2-TZVPP):
* Stage: 4_refinement *
***********************
* Boltzmann population filtering on ΔG* (threshold: 0.950)
* Setup Summary:
- Components: elec=MAIN, gsolv=MAIN, gtrv=AUX
- Main driver: OrcaDriver[wb97x-d3/def2-TZVPP,smd(tetrahydrofuran)]
- Aux driver: XtbDriver[gfn2,gbsa(thf)]
> Conformer #1: -1929.08595168 a.u.
> Conformer #2: -1929.08479631 a.u.
> Conformer #6: -1929.08917852 a.u.
> Conformer #7: -1929.08823314 a.u.
> Conformer #8: -1929.08901520 a.u.
* Final Boltzmann population:
G* ΔG* Pop %
Conformer #1 -1929.08595168 0.00322684 1.5
Conformer #2 -1929.08479631 0.00438221 0.4
Conformer #6 -1929.08917852 0.00000000 44.4 *
Conformer #7 -1929.08823314 0.00094538 16.3 *
Conformer #8 -1929.08901520 0.00016332 37.4 *
* Done, retained 3 conformer(s)
At this stage, the conformer set is reduced to 3 structures, since the two first ones has negligible Boltzmann weight at 298 K.
The final geometries are available in 4_refinement.selected.xyz.