It is not obvious how to increase the distance included in get_bonds. Help would be appreciated! I have below, but want to increase the range included as bonds.
def printbonds(A,B):
ABBonds = ana.get_bonds(A, B, unique=True)
print(“There are {} {}-{} bonds.”.format(len(ABBonds[0]),A,B))
ABBondValues = ana.get_values(ABBonds)
print(*ABBondValues[0])
print(“The average {}-{} bond length is {}.”.format(A,B,np.average(ABBondValues)))
Sincerely,
Ron
help(Analysis) gives
Help on class Analysis in module ase.geometry.analysis:
class Analysis(builtins.object)
| Analysis(images, nl=None, **kwargs)
|
| Analysis class
|
| Parameters for initialization:
|
| images: :class:`~ase.Atoms` object or list of such
| Images to analyze.
| nl: None, :class:`~ase.neighborlist.NeighborList` object or list of such
| Neighborlist(s) for the given images. One or nImages, depending if bonding
| pattern changes or is constant. Using one Neigborlist greatly improves speed.
| kwargs: options, dict
| Arguments for constructing :class:`~ase.neighborlist.NeighborList` object if :data:`nl` is None.
|
This means it will accept arguments to NeighborList and pass them on. Looking at Neighborlist:
Help on class NeighborList in module ase.neighborlist:
class NeighborList(builtins.object)
| NeighborList(cutoffs, skin=0.3, sorted=False, self_interaction=True, bothways=False, primitive=<class 'ase.neighborlist.PrimitiveNeighborList'>)
|
| Neighbor list object.
|
| cutoffs: list of float
| List of cutoff radii - one for each atom. If the spheres
| (defined by their cutoff radii) of two atoms overlap, they
| will be counted as neighbors. See
| :func:`~ase.neighborlist.natural_cutoffs` for an example on
| how to get such a list.
...
So I think you want to use ana=Analysis(atoms, cutoffs=[...]) to set the range of each atom.
Thank you so much! For me this works: ana=Analysis(crys, cutoffs=[5.]*10)
I am not sure what you mean by
| See :func:~ase.neighborlist.natural_cutoffs for an example on
| how to get such a list.
How do I see that?
Ron
I spoke too soon! I wonder if this is a bug. For a simple perovskite structure it will not give the A-A distance. Perhaps there is a hard maximum? I tried smaller values than 10 and they also do the same.
Pm3m_80GPa.xsf (476 Bytes)
xsf2bonds_ase.py (695 Bytes)
Sincerely,
Ron
I also tried:
nl = NeighborList([8.]*len(crys))
nl.update(crys)
ana=Analysis(crys, cutoffs=[8.]*len(crys))
but still doesn’t find Pb-Pb bonds.
Thanks for any help,
Sincerely,
Ron
This gives the info I want, but how do I instead do this through bonds? Thank you!
i = ase.neighborlist.neighbor_list(‘i’, crys,
{(‘Pb’, ‘Pb’): 5})
coord = np.bincount(i)
print(coord)
d = ase.neighborlist.neighbor_list(‘d’, crys,
{(‘Pb’, ‘Pb’): 5})
print(d)
I think the problem is I don’t understand what makes a “bond” in ase, and I don’t see how to change it. Any help is appreciated. Ron
OK, this is really wierd. It seems the bond lengths from ase are incorrect! Attached is a simple code that computes the atomic distances two ways, and there are two different answers. get_values seems to just be incorrect. Using ase.neighborlist.neighbor_list is close but still slightly off. I include a pdf of the code output, the code, the input, and also the output from Mercury from CCDC. What gives?
Sincerely,
Ron
findsym-output-Bonds.csv (1.9 KB)
bonds_ase_test.pdf (23.7 KB)
bonds_ase_test.py (917 Bytes)
R3c_80GPa_findsymR.cif (1.1 KB)
https://wiki.fysik.dtu.dk/ase/ase/neighborlist.html#ase.neighborlist.natural_cutoffs
I think this function can be used to get a cutoff list for use with Analysis by assigning a covalent radius per element, e.g.
cutoffs = natural_cutoffs(atoms, C=5.0, S=4.0, H=2.1)
Thank you so much, but that does not seem to change or fix anything.
See attached.
bonds_ase_test.pdf (23.8 KB)
Ron
You can set the cutoff range of each atom individually in NeighborList(). (ase.neighborlist.py starting at line 1076)
cutoffs: list of float
List of cutoff radii - one for each atom. If the spheres
(defined by their cutoff radii) of two atoms overlap, they
will be counted as neighbors. See
:func:`~ase.neighborlist.natural_cutoffs` for an example on
how to get such a list.
skin: float
If no atom has moved more than the skin-distance since the
last call to the
:meth:`~ase.neighborlist.NeighborList.update()` method, then
the neighbor list can be reused. This will save some
expensive rebuilds of the list, but extra neighbors outside
the cutoff will be returned.
self_interaction: bool
Should an atom return itself as a neighbor?
bothways: bool
Return all neighbors. Default is to return only "half" of
the neighbors.
I don’t use Analysis.get_bonds, but you can look here for an example that gets bonds without the Analysis module. (get_bondpairs() in ase.io.pov.py line 35)
With that, you can compare the Analysis module results to a simplified piece of code.
I hope that helps.
-Mike
Thank you! But it doesn’t help. Still same wrong results. Something is not right with ase.
Ron
bonds_ase_test.pdf (23.8 KB)
I have a bit of trouble interpreting the outputs. Would you be able to make a completely specific statement of the form: “The distance between atom [A1] and atom [A2] computed by [FEATURE] is [BAD_DISTANCE], but their actual distance is [GOOD_DISTANCE].” Then I have something to work with.
Thanks
Ask
Let me have a crack at this.
The distance between Pb (atom 0) and Pb (atom 0 but in some periodic shifted cell (nonzero S vector using neighbor_list convention in ASE/matscipy) should be 3.61331247. This is the lattice constant of the cell in question.
Analysis.get_bonds(‘Pb’, ‘Pb’) is not finding anything, even with a large cutoff. Same goes for a get_bonds call with (‘Ti’, ‘Ti’).
The problem here is that this is a bit counterintuitive (unless bonds are specifically defined in some way that is 1st NN).
If I use ASE’s neighbor_list to call (with ‘ijdS’ as quantities), I get the following:
I’ve highlighted the important parts, but the summary is that we correct detect the Pb-Pb bond (one of many equivalent) shifted by one cell vector. However, it appears that the get_bonds() API doesn’t do this, even with the Atoms object having pbc=True.
So best case, this is quite counterintuitive (for some baseline of intuitive), at worst this is unintended behaviour.
Ther correct Pb-Ti distances for one cell are:
1 4 Pb-Ti 3.02366 A
2 3 Pb-Ti 3.02366 A
2 4 Pb-Ti 3.11663 A
1 3 Pb-Ti 3.11663 A *
1 3 Pb-Ti 3.12881 A *
1 3 Pb-Ti 3.12881 A *
2 4 Pb-Ti 3.12881 A *
2 4 Pb-Ti 3.12881 A *
2 4 Pb-Ti 3.12881 A *
1 3 Pb-Ti 3.12881 A *
1 4 Pb-Ti 3.15907 A *
1 4 Pb-Ti 3.15907 A *
2 3 Pb-Ti 3.15907 A *
2 3 Pb-Ti 3.15907 A *
1 4 Pb-Ti 3.15907 A *
2 3 Pb-Ti 3.15907 A *
There are 16 numbers. The correct average is 3.12549 .
ase.neighborlist.neighbor_list gives 32 which is twice too many. It gives the correct average.
ana.get_bonds and ana.get_values gives the wrong set of numbers , some are wrong like 3.11663475040000 is not a correc value, and the average is incorrect.
R3c80GPaPbTidist.dat (495 Bytes)
bonds_ase_test_P1.pdf (24.2 KB)
bonds_ase_test_P1.ipynb (3.8 KB)
bonds_ase_test_P1.py (970 Bytes)
After some head banging herte is a script that works (but doesn’t use the bonds api). It may not work in all cases or may need adjustment. It is not ideal. But it gives correct answers for me.
Ron
xsf2bonds_ase2.pdf (23.1 KB)
xsf2bonds_ase2.ipynb (4.7 KB)
xsf2bonds_ase2.py (875 Bytes)
P4mm_80GPa_findsym.cif (1.2 KB)
