Issue with constrained optimization in Sella

[csoulie]

Hi,

I am trying to do a constrained optimization, but Sella crashes after a few iterations. Please find attached below an input example for this issue:

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

myatoms = Atoms(numbers=[1,8,1], positions=np.array([[0., -0.5, -0.5], [0., 0., 0.5], [0., 0.5, -0.5]]) ) changes = [[0, 1, 1.0], [1, 2, 1.0], [0, 1, 2, 105]]

calc = EMT() myatoms.calc = calc sella_cons = Constraints(myatoms)

for frozen_coord in changes: match len(frozen_coord[:-1]): case 2: sella_cons.fix_bond(frozen_coord[:-1], target=frozen_coord[-1]) case 3: sella_cons.fix_angle(frozen_coord[:-1], target=frozen_coord[-1]) case 4: sella_cons.fix_dihedral(frozen_coord[:-1], target=frozen_coord[-1])

internals = Internals(myatoms, cons=sella_cons)

dyn = Sella(myatoms, order=0, internal=internals, trajectory="h2o.traj") dyn.run(1e-3, 1000)`

I think the problem is coming from the project() method in ApproximateHessian, where Bproj is returned as an empty (0,0) array. See the screenshot attached.

Am I doing something wrong?

[csoulie]

Building the constrains directly into the Atoms object works:

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

myatoms = Atoms(numbers=[1,8,1],
                        positions=np.array([[0., -0.5, -0.5],
                                            [0., 0., 0.5],
                                            [0., 0.5, -0.5]])
                        )
changes = [[0, 1, 1.0], [1, 2, 1.0]]#, [0, 1, 2, 105]]

calc = EMT()
myatoms.calc = calc
#sella_cons = Constraints(myatoms)

bonds = []
angles_deg = []
dihedrals_deg = []
for frozen_coord in changes:
    match len(frozen_coord[:-1]):
        case 2:
            bonds.append([frozen_coord[-1], frozen_coord[:-1]])
            #sella_cons.fix_bond(frozen_coord[:-1], target=frozen_coord[-1])
        case 3:
            angles_deg.append([frozen_coord[-1], frozen_coord[:-1]])
            #sella_cons.fix_angle(frozen_coord[:-1], target=frozen_coord[-1])
        case 4:
            dihedrals_deg.append([frozen_coord[-1], frozen_coord[:-1]])
            #sella_cons.fix_dihedral(frozen_coord[:-1], target=frozen_coord[-1])
ase_constraints = FixInternals(bonds=bonds, angles_deg=angles_deg, dihedrals_deg=dihedrals_deg)
myatoms.set_constraint(ase_constraints)
#internals = Internals(myatoms, cons=sella_cons)

dyn = Sella(myatoms,
            order=0,
            #internal=internals,
            trajectory="h2o.traj")
dyn.run(1e-3, 1000)

But it only works if not all internal coordinates are specified (doesn't work if the angle constraint is uncommented). Any solutions?

Best,

Clément

[juditzador]

Yes, I was going to post the same solution, here is mine:

from ase import Atoms
from ase.calculators.emt import EMT
from sella import Sella, Constraints, Internals
myatoms = Atoms(numbers=[1,8,1],
          positions=np.array([[0., -0.5, -0.5],
          [0., 0., 0.5],
          [0., 0.5, -0.5]])
          )
print('original internals')
print(myatoms.get_distance(0, 1))
print(myatoms.get_distance(1, 2))
print(myatoms.get_angle(0, 1, 2))
changes = [[1, 2, 1.0], [0, 1, 2, 105]]
calc = EMT()
myatoms.calc = calc
sella_cons = Constraints(myatoms)
for frozen_coord in changes:
    match len(frozen_coord[:-1]):
        case 2:
            sella_cons.fix_bond(frozen_coord[:-1], target=frozen_coord[-1])
        case 3:
            sella_cons.fix_angle(frozen_coord[:-1], target=frozen_coord[-1])
        case 4:
            sella_cons.fix_dihedral(frozen_coord[:-1], target=frozen_coord[-1])
internals = Internals(myatoms, cons=sella_cons)
dyn = Sella(myatoms,
            order=0,
            constraints=sella_cons,
            trajectory="h2o.traj")
dyn.run(1e-3, 1000)
print('final internals')
print(myatoms.get_distance(0, 1))
print(myatoms.get_distance(1, 2))
print(myatoms.get_angle(0, 1, 2))

So yes, you need to either pass the constraints, or not pass anything. The point is not about specifying all internals, they are all specified, but they are all constrained, and that's probably asking too much of the code? You could achieve the desired outcome in two steps perhaps. I think there is nothing to be optimized when all coordinates are constrained and that might be the problem.