iterativeOptimization.cmp seems to hang using GFN2-xTB

[camattelaer]

Hello

I am trying to use the iterativeOptimization.cmp script to optimize until no negative frequencies are observed.

the command ${ORCAHOME}/orca iterative_opt.inp > iterative_opt.out seems to stall . The last thing that is being executed is writing iterative_opt.citations.tmp.

the iterative_opt.inp simply contains:

$ cat iterative_opt.inp 
%Compound "iterativeOptimization.cmp"

my (modified) iterativeOptimization.cmp (method and system info - xyz input file, charge - are modified):

# Author: Dimitrios G. Liakos, Zikuan Wang and Frank Neese
# Date  : May/June of 2024
# Edited by James O'Brien, November 2024
#
# *************************************** DESCRIPTION ***********************************************
# iterative Optimization protocol to find structure with no negative
# frequencies (i.e. real minima), or in case of transition state search,
# find structure with exactly one imaginary frequency (i.e. saddle point)
#
# Step 1. Run a single point calculation (we need it for the first property file)
#
# Step 2. Loop and perform calculations with (optimization and frequencies)
#
# Step 3. Check the frequencies. If there are negative ones use the hessian
#         of the appropriate normal mode to adjust the geometry
#
# ------ Variables to adjust (e.g. using 'with') -------------------
Variable molecule = "/home/camattelaer/mount_data/data02_18TB/website/01_HNAduplex/01_hhA/xyzInput.xyz";
Variable charge = -4;
Variable mult   = 1;
Variable method          = "GFN2-XTB ALPB(water)";
Variable MaxNTries       =  25;   # Number of maximum tries
Variable CutOff          = -10.0; # CutOff for a negative frequency
Variable scaling         = 0.25;   # Scaling factor for normal mode (WARNING: large values can cause molecules to "blow up")
Variable NNegativeTarget = 0;     # Number of negative frequencies wanted. For TS optimization, change to 1
Variable myFilename      = "xyzInput.xyz";
# ------------------------------------------------------------------
# ------            Rest of variables            -------------------
Geometry myGeom;
Variable freqs, modes;
Variable res = -1;
Variable NNegative =   0;
Variable OptDone;

# -----------------------------------------------------------
# Perform a single point calculation. We need it for
# the initial geometry from the property file
# -----------------------------------------------------------
NewStep
  !&{method}
  *xyzfile &{charge} &{mult} &{molecule}
StepEnd
myGeom.Read();
myGeom.WriteXYZFile(filename=myFilename);

# -----------------------------------------------------------
# Start a for loop over number of tries
# ----------------------------------------------------------
For itry From 1 To maxNTries Do
  # --------------------------------------------
  # Perform a geometry optimization/Frequency calculation
  # --------------------------------------------
  New_Step
    ! &{method} freq Opt
    *xyzfile &{charge} &{mult} &{myFilename}
  Step_End
  res = freqs.readProperty(propertyName = "THERMO_FREQS");
  res = modes.readProperty(propertyName = "HESSIAN_MODES");
  myGeom.Read();

  # ---------------------------------------------
  #  check for sufficiently negative frequencies
  # ---------------------------------------------
  NNegative = 0;
  For ifreq From 0 to freqs.GetSize()-1 Do
    if ( freqs[ifreq] < CutOff )  then
       myGeom.FollowNormalMode(vibrationSN=ifreq+1, scalingFactor=scaling);
       NNegative = NNegative + 1;
    endif
  endfor
  myGeom.WriteXYZFile(filename=myFilename);
  If ( NNegative <= NNegativeTarget ) then
    goto OptDone;
  endif
endfor

# -----------------------------------------------------------------
# Either found correct geometry or reached maximum number of tries.
# -----------------------------------------------------------------
OptDone :
if (NNegative > NNegativeTarget) then
  print("ERROR The program did not find a structure with the desired\n number of imaginary frequencies.\n There are %d negative frequencies after %3d steps,\n which is larger than the desired number %d.", NNegative,itry,NNegativeTarget);
else if (NNegative < NNegativeTarget) then
  print("ERROR The program did not find a structure with the desired\n number of imaginary frequencies.\n There are %d negative frequencies after %3d steps,\n which is smaller than the desired number %d.", NNegative,itry,NNegativeTarget);
else
  print("\nSUCCESS optimized structure with (%d) negative\n frequencies found after %3d steps", NNegative, itry);
endif

End

The last thing written to the orca job output file are timings:

...
For completeness - the Gibbs free energy minus the electronic energy
G-E(el)                           ...      1.42329356 Eh    893.13 kcal/mol


Timings for individual modules:

Sum of individual times          ...     1322.421 sec (=  22.040 min)
Geometry relaxation              ...       12.162 sec (=   0.203 min)   0.9 %
XTB module                       ...     1310.260 sec (=  21.838 min)  99.1 %

I suspect this might be related to this bug (it would be still present in orca v6...) : https://orcaforum.kofo.mpg.de/viewtopic.php?f=55&t=9490 which causes the reading of properties/Hessian to be extremly slow/bugged ?

Kind regards and thanks in advance

Charles-Alexandre