Question about relaxed PES scan with Sella

[yangyinuo823]

Does Sella support relaxed PES scan in a molecular system? I suppose I need to do optimization with Constraints, but not sure how to scan the coordinates such as bond distances. Highly appreciated for any advices

[juditzador]

Yes, you need to invoke a constrained optimization for the coordinate you wish to scan, and write an outer loop that changes that coordinate as you like. You need to use target to set the desired value of the coordinate in each step. See https://github.com/zadorlab/sella/wiki/Constraints

[yangyinuo823]

This is very helpful, thank you!

[yangyinuo823]

Hi, I have a followup question on relaxed PES scan with internal coordinates. Below is a mininal test example on scaning the O-O bond in HOOH, following the instrustions here https://github.com/zadorlab/sella/wiki/Constraints:

from ase import Atoms
from ase.calculators.emt import EMT
from sella import Sella, Constraints, Internals

mol = Atoms('HOOH', positions=[(0,0,0),(1,0,0),(1,0,1),(1,1,1)])
mol.calc = EMT()
scan_dis = np.linspace(1, 2, 5)
for i, di in enumerate(scan_dis):
    print('%sth iteration with O-O bond distance %f'%(i, di))
    cons = Constraints(mol)
    cons.fix_bond((1,2), target = di)
    internals = Internals(mol, cons = cons)
    internals.add_bond((1,2))
    dyn = Sella(mol, internal = internals)
    dyn.run(1e-3,1000)

But it always gets stuck at the second iteraction with below error

0th iteration with O-O bond distance 1.000000
No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/internal.py:1559: UserWarning: 1 coords found! Expected 6.
  warnings.warn(
     Step     Time          Energy         fmax         cmax       rtrust          rho
Sella   0 12:37:17        3.542592       0.0000       0.0000       0.1000       1.0000
1th iteration with O-O bond distance 1.250000
     Step     Time          Energy         fmax         cmax       rtrust          rho
Sella   0 12:37:17        3.542592       0.0000       0.2500       0.1000       1.0000
Traceback (most recent call last):
  File "PES_minimal.py", line 19, in <module>
    dyn.run(1e-3,1000)
  File "/home/yinuo/.local/lib/python3.8/site-packages/ase/optimize/optimize.py", line 269, in run
    return Dynamics.run(self)
  File "/home/yinuo/.local/lib/python3.8/site-packages/ase/optimize/optimize.py", line 156, in run
    for converged in Dynamics.irun(self):
  File "/home/yinuo/.local/lib/python3.8/site-packages/ase/optimize/optimize.py", line 135, in irun
    self.step()
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/optimize/optimize.py", line 238, in step
    s, smag = self._predict_step()
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/optimize/optimize.py", line 217, in _predict_step
    self.pes.diag(**self.diagkwargs)
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/peswrapper.py", line 252, in diag
    P = self.get_HL().project(Ufree)
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/linalg.py", line 202, in project
    return ApproximateHessian(n, 0, Bproj, self.update_method,
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/linalg.py", line 155, in __init__
    self.set_B(B0)
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/sella/linalg.py", line 170, in set_B
    self.evals, self.evecs = eigh(self.B)
  File "/home/yinuo/miniforge3/envs/torchani/lib/python3.8/site-packages/scipy/linalg/_decomp.py", line 547, in eigh
    w, v, *other_args, info = drv(a=a1, **drv_args, **lwork_args)
_flapack.error: (il>=1&&il<=n) failed for 6th keyword il: dsyevr:il=1

Any suggestions on how to fix this are high appreciated!

[juditzador]

You don't have the internal coordinates set up correctly, you only added one. Try this:

import numpy as np
from ase.io import Trajectory
from ase import Atoms
from ase.calculators.emt import EMT
from sella import Sella, Constraints, Internals

mol = Atoms('HOOH', positions=[(0,0,0),(1,0,0),(1,0,1),(1,1,1)])
mol.calc = EMT()
scan_dis = np.linspace(1, 2, 5)
for i, di in enumerate(scan_dis):
    print('%sth iteration with O-O bond distance %f'%(i, di))
    cons = Constraints(mol)
    cons.fix_bond((1,2), target=di)
    internals = Internals(mol, cons = cons)
    internals.find_all_bonds()
    internals.find_all_angles()
    internals.find_all_dihedrals()
    internals.add_bond((1,2))
    dyn = Sella(mol, internal=internals, trajectory=f'h2o2_{i}.traj')
    dyn.run(1e-3,1000)

[yangyinuo823]

this is exactly what I need, thank you so much!