Implement Linear bends

[jhrmnn]

• [ ] Implement ghost atoms.

This might fix #23. It might also make #9 obsolete.

[mjw99]

Indeed; and may also replace #21

[jhrmnn]

I'll add the representation of ghost atoms, which is trivial, then feel free to implement the linear bends if you're up to it. That will be quite a bit of work though.

[jhrmnn]

Internally, I would represent the ghost atoms as regular atoms of species "_". The only unclear thing to me is whether the current API, such as the coords property or iterating over the Geometry object should also include the ghost atoms, or whether they should be retrievable (together with the regular atoms) via some new API, such as a centers property. @mjw99 Your opinion?

[mjw99]

I am not sure at the moment.

Initially, I was going to extend Angle, i.e. class LinearBend(Angle) and then choose to populate these as a function of angle in InternalCoords.__init__(). I was thinking something (this is a crude mock-up off the top of my head), like this:

        for j in range(len(geom)):
            for i, k in combinations(np.flatnonzero(bondmatrix[j, :]), 2):
                ang = Angle(i, j, k, C=C)

                if ang.eval(geom.coords) > (165 * (2*pi/360)):
                    # place ghost atom 
                    #np.append(geom.coords, [0,0,0])
                    #np.append(geom.species, "TODO")

                    # Place first LinearBend with ghost as k
                    lbend = LinearBend(i, j, d, C=C)
                    self.append(lbend)

                    # Place second LinearBend with ghost as i
                    lbend = LinearBend(k, j, d, C=C)
                    self.append(lbend)

                    # Update position of ghost to satify constraints
                    # ijd == kjd, r2d = 2 bohr
                else:
                    if ang.eval(geom.coords) > pi/4:
                       self.append(ang)

However, I am not sure how to update geom.coords whilst in this loop, or even if I should be doing it like that. Also, apologies for my non pythonic style. (edit, spelling and apology)

[mjw99]

Yeah... this was quite a naive first attempt. By attempting to add ghost atoms, I am essentially changing the state of Geometry as I iterate over it in InternalCoords.__init__().

[jhrmnn]

Right. The dummy atoms do not fit well in the current factoring. On one hand, they are not really part of the structure, so they shouldn't be in Geometry, on the other, they do carry a Cartesian coordinate, so they should be there. I think the best resolution is to keep them in InternalCoords. The only disadvantage of this is that the Cartesian information will be in two different places, so one won't be able to, for example, rotate the molecule in the middle of an optimization. But I guess that's ok. The API will require that the Geometry object yielded by Berny cannot be modified in-place.

So I suggest that you add InternalCoords.dummy_atoms which will be an Nda-by-3 Numpy array of the coordinates of the dummy atoms. I don't think that a new class for the linear bend is necessary—the set of internal coordinates will still consist only of bonds, angles, and dihedrals. InternalCoords will just need to remember which coordinates belong to a given linear bend to do all the machinery described in the paper.

[askhl]

Hi @jhrmnn. Unfortunately, at least as far as I am aware, we can't very well use the optimizer for production calculations in ASE as it will normally crash due to this issue.

I'd therefore be inclined to remove the feature from ASE so we don't need to keep infrastructure and tests passing. But should the issue be solved, I'd be very happy to add the interface again. If there are nevertheless good reasons to keep it in ASE, please let me know.

[jhrmnn]

Hi Ask, understood. I unfortunately don’t have the bandwidth in the foreseeable future to finish this off, so someone else would have to. I understand the incentive to remove the interface from ASE.Though I will point out that other than this issue, Berny seems to be doing well: https://pubs.acs.org/doi/10.1021/acs.jctc.3c00188                                              

[askhl]

Thanks @jhrmnn. Those numbers are indeed impressive. However I don't quite understand: When I run calculations, they often crash. I therefore see the problem as a blocker, at least for serious production calculations. This was also the take-home message I got from conversations here. But apparently that didn't stop the optimizer from performing exceedingly well in the publication you posted.

So why is it that I get trouble and others do not? I think one contribution was linear bends, and another was that maybe the forces are 'unclean' (e.g. affected by things like grids). If I only see problems because of a kind of bias, then it would probably be wrong to remove the optimizer.

[jhrmnn]

Let me reach out to the authors, and perhaps even loop them into the discussion here. I haven’t had any contact with them so far.