samblau
commented
5 months ago FYI @ehermes @mgt16-LANL @Andrew-S-Rosen
Closed samblau closed 5 months ago
samblau
commented
5 months ago FYI @ehermes @mgt16-LANL @Andrew-S-Rosen
juditzador
commented
5 months ago This is not fixing the code, but one way around this, if indeed the problem is the unconnected part, is that you detect that the bond graph is disconnected before running Sella. You can get all internals by internals = Internals(myatoms) and then add missing connections either as bonds or by other types of coordinates - this could be automated then. Did you try something like this?
juditzador
commented
5 months ago Here adding manually a bond to connect a fragment to the central atom (1, 61) - it doesn't throw an error anymore, the optimization starts. Nevertheless, I understand that the code should in principle do this automatically, and it does by the flood algorithm. I'll see if I can find the bug.
from sella import Sella
from sella import Internals
from ase.calculators.lj import LennardJones
from ase.io import read
test_atoms = read("test.xyz")
internals = Internals(test_atoms)
internals.add_bond((1, 61))
internals.find_all_bonds()
internals.find_all_angles()
internals.find_all_dihedrals()
test_atoms.calc = LennardJones()
dyn = Sella(test_atoms, internal=internals, order=0)
dyn.run()I think I know what is going on.
def find_all_bonds(
self,
nbond_cart_thr: int = 6,
max_bonds: int = 20,
) -> None:sets the number of bonds to be max 20. I'd say that's not a very strict limit, but...
If I run the above code with three modifications:
from sella import Sella
from sella import Internals
from ase.calculators.lj import LennardJones
from ase.io import read
test_atoms = read("test.xyz")
internals = Internals(test_atoms, allow_fragments=True)
internals.find_all_bonds()
for bond in internals.internals['bonds']:
print(bond)
internals.find_all_angles()
internals.find_all_dihedrals()
test_atoms.calc = LennardJones()
dyn = Sella(test_atoms, internal=internals, order=0)
dyn.run()The following things happen: (1)+(2) results in a warning only and the optimization runs:
/home/jzador/.local/lib/python3.10/site-packages/sella/internal.py:1559: UserWarning: 291 coords found! Expected 285.(285 is 3*97-6)
And (3) reveals that there are indeed 19 bonds to the central atom already. It is because I guess this structure is not simply fragmented, but with the ballooning radius of the flood algorithm at trial radius R you have still fragments, while at the next trial radius, R * 1.05, you suddenly have more than 20 bonds. This explains why it happens only occasionally. If something is oriented slightly differently, the inclusion of the fragment atoms is more gradual, but sometimes you get unlucky.
OK, what to do. A possible bandaid is to say allow_fragments=True and generate the internals the way I do it. The optimization is still running, and the LJ is not supposed to be good anyway, so you need to see if this bandaid is good for now.
I would imagine that raising the threshold for max_bonds is a bad idea, from a nicely connected easy system now you create a densely connected one, where the internal coordinates lost their chemical meaning. The real solution is to fix this bug, i.e., sella could go more gentle in adding atoms. Of course, if you have a perfect sphere with 20 atoms on its surface, each far from each other, and there is a 21st atom in the middle, it'll still fail, but then maybe you are not doing chemistry - your system is not like that. OK, too many words here, I'll see if there is a good fix and consult with @ehermes if he has ideas.
juditzador
commented
5 months ago Another working version of the script:
from sella import Sella
from sella import Internals
from ase.calculators.lj import LennardJones
from ase.io import read
test_atoms = read("test.xyz")
internals = Internals(test_atoms)
internals.find_all_bonds(max_bonds=40)
for bond in internals.internals['bonds']:
print(bond)
internals.find_all_angles()
internals.find_all_dihedrals()
test_atoms.calc = LennardJones()
dyn = Sella(test_atoms, internal=internals, order=0)
dyn.run()Here increasing max_bonds to 40 (you only need 22 in fact), but not allowing fragments.
I am currently using Sella to run tens of thousands of optimizations of metal-organic coordination complexes. in ~0.8% of cases, we see an error like this during the initialization of the internal coordinates:
This appears to occur for some structures where ligands are not obviously coordinated to the central metal atom. Since I apparently cannot attach XYZ files to this issue, I will provide one example structure here, which is a -2 singlet:
Obviously I do not suggest using LennardJones for this type of system, but here is code that can reproduce the error quickly and easily reading in this structure:
I could easily provide at least fifty other example structures if it would be helpful for debugging.