IRC Failing to Converge Silently on Constrained Slab TSs

[mjohnson541]

I'm having trouble with the IRCs in metal slab systems with the metal atoms frozen for many different transition states. All of the log files look something like the below. Only the 0 step is logged. Some warnings show up in the error file. The calculation finishes, but the output irc is not marked as converged.

Out file:

 Step     Time          Energy         fmax

IRC: 0 02:47:38 -238815.926294 0.2857

Error file:

/etc/profile.d/gaussian.sh: line 8: [: too many arguments
/etc/profile.d/gaussian.sh: line 19: [: too many arguments
WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
/home/mjohns9/anaconda3/envs/pynta_env/lib/python3.8/site-packages/sella/peswrapper.py:325: RuntimeWarning: invalid value encountered in double_scalars

The last espresso file seems to have terminated cleanly. And best as I can tell the IRC calculations don't run into any errors within Sella itself. I essentially ran:

opt = IRC(sp,constraints=cons,trajectory=label+"_irc.traj",dx=0.1,eta=1e-4,gamma=0.4)
try:
    run_kwargs = {"direction": "forward"}
    opt.run(**run_kwargs)
    run_kwargs = {"direction": "reverse"}
    opt.run(**run_kwargs)
except Exception as e:
    if not ignore_errors:
        raise e
    else:
        errors.append(e)

Naturally I went back and tried a simple gas phase example C2H5 => C2H4 + H, which largely worked even if it errored explicitly eventually. I then tried to do the above example with GFN1-xTB to make it easier to look at. The IRC errors in that case with IRCInnerLoopConvergenceFailure, but from the DFT guess GFN1-xTB after optimization has two imaginary frequencies so that's probably expected.

Since the gas phase example worked it seems like it might be related to the constraints implementation in Sella.

[ehermes]

Can you provide a minimal example that reproduces this issue? xTB should be fine, but I'd like to avoid having to run PWDFT calculations.

[mjohnson541]

from ase.io import read
from ase import Atoms
from sella import Sella, Constraints, IRC
from xtb.ase.calculator import XTB
from ase.visualize import view
from ase.io.trajectory import Trajectory


import numpy as np
from ase.vibrations import Vibrations

sp = read("opt.xyz")
sp.calc = XTB(method="GFN1-xTB")
cons = Constraints(sp)

for i,atom in enumerate(sp): #freeze the slab
    if i < 3*3*4:
        cons.fix_translation(atom.index)
opt = Sella(sp,constraints=cons,trajectory="xtbtest.traj",order=1)

opt.run(fmax=0.01)

irc = IRC(sp,constraints=cons,trajectory="test_irc.traj",dx=0.1,eta=1e-4,gamma=0.4)

irc.run(direction="forward",fmax=0.1,steps=1000)
irc.run(direction="reverse",fmax=0.1,steps=1000)

This one seems to have only one imaginary frequency after optimizing with xtb.
opt.txt

[ehermes]

The issue is that the constraints aren't getting set properly. I'll look into fixing it, which should be simple, but in the meantime here's a workaround:

from ase.io import read
from ase import Atoms
from sella import Sella, IRC
from xtb.ase.calculator import XTB
from ase.visualize import view
from ase.io.trajectory import Trajectory
from ase.constraints import FixAtoms


import numpy as np
from ase.vibrations import Vibrations

sp = read("opt.xyz")
sp.calc = XTB(method="GFN1-xTB")
sp.set_constraint(FixAtoms([atom.index for atom in sp if atom.symbol == 'Cu']))

opt = Sella(sp,trajectory="xtbtest.traj",order=1)

opt.run(fmax=0.01)

irc = IRC(sp,trajectory="test_irc.traj",dx=0.1,eta=1e-4,gamma=0.4)

irc.run(direction="forward",fmax=0.1,steps=1000)
irc.run(direction="reverse",fmax=0.1,steps=1000)

(You don't have to use my exact construction for setting the constraints, the key is to use ASE's FixAtoms constraint instead of Sella's Constraints object)

[ehermes]

Ok, the constraints are getting set properly, but the problem is that the IRC module is wholly ignorant of the constraints. The reason my workaround works is because ASE constraints project out any forces that would violate the constraints, which is not something we want for general constraints but this works fine for fixed atom constraints in Cartesian coordinates. So the action items to fixing this are the following:

1. Make sure ASE constraints are REMOVED and replaced by Sella Constraints in the IRC class (there should be no difference in behavior between the different ways of defining constraints, so the fact that this workaround works at all is a bug)
2. Verify that the only constraints in use are translation constraints, and raise an error otherwise
3. In IRC.step, zero the components of the force for fixed atoms at every step

That should be sufficient for fixing this bug.