The GöNEB code uses the framework of the Nudged Elastic Band method described by Jonsson et.al.

Note

Only tested with python 3.8, and the exact versions of the libraries given in requirements.txt!

Installation

1. Clone the repository

   git clone https://github.com/cchembio/goeneb.git
   cd goeneb

---

2. Create and activate a virtual environment

   1. Select your Python version
      Ensure you have your desired Python version installed on your system (Python 3.8 is recommended):

      python --version

   2. Create the virtual environment

      python -m venv neb-env

      Replace neb-env with any name you prefer for your environment.

   3. Activate the environment

      source neb-env/bin/activate

   4. Install required Python modules

      Install the dependencies from the requirements.txt file:

      pip install -r requirements.txt

---

3. Prepare the shell environment

   The file scripts/environment.sh is used to generate a temporary directory for the calculations, export all environment variables and source the python environment. The lines

   ENV_PATH="/path/to/neb-env/bin/activate"
   NEB_PATH="/path/to/goeneb"

   need to be changed to the full paths of the python environment and the GöNEB program.

Warning

The environment.sh expects the executable of the quantum chemistry program you are planning to use (ORCA, Molpro, Gaussian) to be in the PATH. Either permanently, or e.g. by loading the module. If this is not the case, the path to the program needs to be specified here.

---

4. Test your installation

   The scripts/test_NEB.sh can now be used to start the test module (testing_module.py). It tests all the energy calculation interfaces (ORCA, Molpro, Gaussian) and checks whether the calculated values for energies and gradients are similar to the values calculated on our system.

   Furthermore, it runs the important functions from the NEB routine together with the dummy potential and tests whether all functions return the correct values. The tested functions are:

   • Alignment
   • Interpolation (Cartesian, Internal (z-matrix), geodesic (non-deterministic) and IDPP)
   • NEB gradients (tangents, springs, complete NEB-gradient)
   • Step predictors (SD, AMGD, RFO, SCT)

   The test returns a file test.log to your current directory. You should get CHECK for all tested functions (except geodesic).

---

The GöNEB program is now ready to use!

Calculating NEBs

1. Set up a job directory.

   OPTION 1: Use the whole directory as input for the GöNEB

   • This directory must include exactly one *.ini file with the options for the NEB
   • It is possible to define start, end and TS structures in the ini-file via their path. This way the xyz-structures may be located elsewhere
   • It is also possible to not specify the input structures. When choosing this option you should name your structures *1*.xyz for the start structure, *2*.xyz for the end structure and *TS*.xyz for the TS if you want to include this structure. If your naming was not recognized by the program, you will get a warning and the structures are sorted alphabetically.

   OPTION 2: Use the *.ini file as the only input for the GöNEB program. You need to specify at least the start and end structure via their path in the ini file. They may be located anywhere on your system.

Important

You can find further information on how to write the ini file later in this ReadMe.

Also: An example for a complete calculation is in the folder goeneb/testjobs/example. There you can see a complete input and output.

---

2. Run the program

   Run the scripts/run_NEB.sh. This sets up the environment, runs the program and then deletes the temporary files.

---

The program will produce:

• batchlog.out: The slurm output, preferably this file should stay almost empty
• output.log: The NEB output, all info you need to follow the NEB calculation
• results:
  • starttraj.xyz: The starting trajectory of the NEB
  • currenttraj.xyz: The current trajectory during the NEB calculation. Continuously updated while the job is running.
  • finaltraj.xyz: The final trajectory (not necessarily converged)
  • optlog.csv: A csv file with different properties for tracking the calculation that can be visualized with the optplot.py program
  • HEI.xyz and HEI_trj.xyz: The highest energy image and its trajectory over the course of the NEB optimization
  • TS_guess.xyz and TS_trj.xyz: A linear interpolation of the images with the highest energy, that should be closest to the actual TS and its trajectory over the course of the NEB optimization

---

Keywords

The neb.ini file should contain all options of the NEB calculation. It starts with the line:

[options]

Each subsequent line should be structured like:

keyword = value

All keywords are explained in the following tables, default values are printed bold.

Initial Set-up and interpolation

Keyword	Options	Explanation
start_structure	path-to/struct.xyz	The start structure, only explicitly needed, if the input is only this file, not the whole directory.
end_structure	path-to/struct.xyz	The start structure, only explicitly needed, if the input is only this file, not the whole directory.
TS_guess	None path-to/struct.xyz	The TS structure, will be included in initial the interpolation.
starttraj	None path-to/traj.xyz	The file with the starting path, from where the NEB should be started. None means, a starting path is generated.
n_images	default: 11 integer $>$ 0	Number of interpolation images. The total number of images (including ends) will be n_images + 2
interp_mode	internal cartesian geodesic	How the initial path is created. Internal uses the z-matrix from the ChemCoord-package, geodesic uses the code by Zhu et.al.
trajtest	False True	Whether to stop after the initial trajectory has been created.
charge	default: 0 any integer	The charge of the system.
spin	default: 1 any integer	The spin of the system.
failed_img_tol_percent	default: 1.0 float between 0 and 1	The maximum fraction of images that can have failed calculations before the NEB aborts.
rot_align_mode	pairwise single_reference	How the rotation is removed from the initial path. Rotation is removed via Kabsch algorithm, either pairwise to the neighboring image, or single reference to the start structure.
remove_gradtrans	True False	Whether or not translations are included in the calculation of the NEB-gradient.
remove_gradrot	False True	Whether or not rotations are included in the calculation of the NEB-gradient.
verbose	debug info warning error critical	The logging level (levels from the logging module). In how much detail should output.log be written.

IDPP keywords

Keyword	Options	Explanation
IDPP	True False	Whether an IDPP pre-optimization is done after the initial interpolation.
SIDPP	False True	Whether a Sequential IDPP should be used for the initial path.
IDPP_maxiter	default: 1000 any integer	Maximum iterations during the IDPP or SIDPP pass.
IDPP_max_RMSF	default: 0.00945 any float	RMSF convergence criteria for the IDPP. In $E_h/\AA$.
IDPP_max_AbsF	default: 0.0189 any float	Maximum absolute value convergence criteria for the IDPP. In $E_h/\AA$.

Interface selection

Keyword	Options	Explanation
interface	orca gaussian molpro dummy	What program should be used for the energy calculation. Dummy is a fast to calculate, non-chemical surface. Only for debugging.
gaussian_path molpro_path orca_path	None path-to/program	This will override the GAUSS_EXE, MOL_EXE or ORCA_EXE environment variables used normally by the program.
gaussian_keywords molpro_keywords orca_keywords orca_keywords2	None Any valid combination of keywords for the chosen interface	Gaussian: 'force' keyword is added automatically Molpro: 'SET,CHARGE=...', 'SET,SPIN=...' and 'force' are added automatically ORCA: 'NoAutoStart EnGrad Angs' is added automatically, second line can be added with orca_keywords2
memory	default: 10000 any integer	Memory keyword, needed for Molpro and Gaussian. MOLPRO: memory in MW, remember: this is assigned to each core. GAUSSIAN: memory in MB, remember: this is the total memory
n_threads	default: 1 any integer	The cores used for parallelizing the energy calculation. Has to be the same as the number of cores specified in the slurm file!

NEB specific options

Tangents and Spring Forces

Keyword	Options	Explanation
maxiter	default: 500 integer $>$ 0	Maximum number of iterations of the NEB.
relaxed_neb	False True	Whether to use relaxed convergence thresholds.
climbing_image	False True	Whether or not climbing image should be activated (after an initial NEB optimization to relaxed convergence criteria).
spring_gradient	difference projected raw	The definition of the spring gradient. The difference-based version is the improved definition by Henkelman and Jonsson, projected is the original implementation of the NEB, and the raw spring force is not projected onto the tangent direction.
k_const	default: 0.003 any float	The spring force constant in $E_h/\AA^2$. When using variable k, this is the maximum value for the springs.
use_vark	False	Whether to use the improved variable k scheme based on the Scheme by Henkelmann and Jonsson.
vark_min_fac	default: 0.1 float between 0 and 1	The minimal spring constant between images is vark_min_fac * k_const
use_analytical_springpos	False True	Analytical spring positioning scheme. Calculates the ideal position of the images in the 1D case. Does not work with SCT step predictor.
tangents	henkjon simple	The type of tangents used. Most importantantly the improved tangent by Henkelman and Jonsson.

Step Predictors

Keyword	Options	Explanation
step_pred_method	AMGD SD SCT RFO NR	Adaptive Momentum Gradient Descent Steepest Descent Self Consistent Tangents Rational Function Optimization with global BFGS Newton Raphson with global BFGS
stepsize_fac	default: 0.2 any float	The step size factor used to scale the step in AMGD, SD, NR. Also in the AMGD steps leading up to any of the other methods (SCT, NR, RF)
max_step	default: 0.05 any float	The maximum step size in $\AA$. The step will be scaled back if it is too long.
AMGD_max_gamma	default: 0.9 float between 0 and 1	The maximum factor for including the last step in the AMGD method.
initial_hessian	diagonal lindh	The initial Hessian is created either as identity or as Lindh hessian. The diagonal identity will be scaled in the first step. The global Lindh hessian is just the diagonal block matrix.
BFGS_start	default: 5 integer $\geq$ 0	The iteration, in which the BFGS update starts.
NR_start	default: 10 integer $\geq$ 0	The iteration, in which the second order method (NR, RF or SCT) starts.
harmonic_stepsize_fac	default: 0.01 any float	The step size used for the AMGD in the sub-cycle of the harmonic NEB (only in the method SCT)
harmonic_conv_fac	default: 0.7 float between 0 and 1	The factor used to lower the convergence thresholds for the harmonic NEB cycle (only in method SCT)

Convergence Thresholds

All convergence Thresholds are the same ones as in ORCA.

Keyword	Options	Explanation
Max_RMSF_tol	default: 0.000945 any float	The RMSF convergence criteria for the normal NEB. In $E_h/\AA$.
Max_AbsF_tol	default: 0.00189 any float	The maximum absolute value convergence criteria for the normal NEB. In $E_h/\AA$.
CI_RMSF_tol	default: 0.000473 any float	The RMSF convergence criteria for the climbing image in the NEB. In $E_h/\AA$.
CI_AbsF_tol	default: 0.000945 any float	The maximum absolute value convergence criteria for the climbing image in the NEB. In $E_h/\AA$.
Relaxed_Max_RMSF_tol	default: 0.00945 any float	Relaxed RMSF convergence criteria for the NEB. This has to be reached for CI to start. In $E_h/\AA$.
Relaxed_Max_AbsF_tol	default: 0.0189 any float	Relaxed maximum absolute value convergence criteria for the NEB. This has to be reached for CI to start. In $E_h/\AA$.

Visualization with optplot

The optplot.py program transforms the optlog.csv file into nice, readable output. To visualize, follow these steps:

1. Open your NEB environment
   

source path/to/neb-env/bin/activate

2. Run the program with
   

path/to/goeneb/optplot.py -f path/to/optlog.csv

The program has the following arguments (can also be accessed by path-to/goeneb/optplot.py --help):

short argument	long argument	Description
-f	--filename	Enter the file that should be visualized. No file opens the user input.
-s	--savefile	Enter the path where the file should be saved. No file means displaying the image directly.

The optplot software gives you 4 different plots to visualize the progress of the NEB:

1. The progression of the convergence measures over the iterations.
2. The energy profile of the reaction, the energy profiles of previous iterations are also shown in lesser occupacity.
3. A 2-D projection of the image coordinates. The initial pathway is also shown and the movement of the images is tracked.
4. The energy gradient for each image in each iteration as a visualized matrix.

The 2D example

The GöNEB comes with a two dimensional NEB model. It uses the functions in the real GöNEB but a model PES which is Himmelblaus function. On this 2D surface the NEB runs instantly and produces results that are easy to analyze since the optimization is in only 2 dimension. This can be used for visualization and testing of various functions.

The code for the 2D example can be run by the following command:

python /path/to/goeneb/2dtest.py -f <filename> -m <step-prediction-method>

The other arguments that can be parsed by the program are:

Short	Long	Description
-s	--stepsize_fac	Step size factor
-m	--step_pred_method	Step predictor method (AMGD, SD, RFO, etc.)
-k	--k_const	Spring constant
-t	--tangents	Tangents calculation method (henkjon, simple)
-i	--maxiter	Maximum number of iterations
-l	--max_step	Maximum step length
-n	--NR_start	Iteration when to start second order method
-b	--BFGS_start	Iteration when to start BFGS updates
-d	--draw_every_iteration	Draw output after every optimization iteration
-h	--show_hessian	Show hessian matrix
-c	--caption	Title of the output graph
-f	--filename	Filename of the produced image
-ff	--fileformat	File format of the produced image (pdf, png, etc.)

Contact

This code is currently maintained and developed by Lynn Meeder (lynn.meeder@uni-goettingen.de).
The base code was written by Björn Hein-Janke (bjoern.hein-janke@uni-goettingen.de).

References

• Jónsson, H.; Mills, G.; Jacobsen, K. W. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations; WORLD SCIENTIFIC: LERICI, Villa Marigola, 1998; pp 385–404.
• Weser, O.; Hein-Janke, B.; Mata, R. A. Automated Handling of Complex Chemical Structures in Z-Matrix Coordinates—the Chemcoord Library. J. Comput. Chem. 2022, 44 (5), 710–726.
• Zhu, X.; Thompson, K. C.; Martínez, T. J. Geodesic Interpolation for Reaction Pathways. J. Chem. Phys. 2019, 150 (16), 164103.
• Smidstrup, S.; Pedersen, A.; Stokbro, K.; Jónsson, H. Improved Initial Guess for Minimum Energy Path Calculations. J. Chem. Phys. 2014, 140 (21), 214106.
• Schmerwitz, Y. L. A.; Ásgeirsson, V.; Jónsson, H. Improved Initialization of Optimal Path Calculations Using Sequential Traversal over the Image-Dependent Pair Potential Surface J. Chem. Theory Comput. 2024, 20 (1), 155–163.
• Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113 (22), 9901–9904.
• Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113 (22), 9978–9985.
• Lindh, R.; Bernhardsson, A.; Karlström, G.; Malmqvist, P.-Å. On the Use of a Hessian Model Function in Molecular Geometry Optimizations. Chem. Phys. Lett. 1995, 241 (4), 423–428.