Displacement Gradient Function for ASE Implementation

[GabrielBram]

I am unsure if this is the correct place to ask so I will try another channel if this repository is no longer supported!

Currently, I am working on an ASE based implementation of the geodesic interpolation method using this repository as a template (mostly with the aim of interpolating reaction pathways with periodic boundary conditions and constraints). I only plan on implementing the initialization and sweeping algorithms for the moment. However, I am a bit stumped by the functions used to calculate the displacement gradient:

    def compute_disp_grad(self, start, end, friction=1e-3):
        """Compute derivatives of the displacement vectors with respect to the Cartesian coordinates"""
        # Calculate derivatives of displacement vectors with respect to image Cartesians
        l = end - start + 1
        self.grad = np.zeros((l * 2 * self.nrij + 3 * (end - start) * self.natoms, (end - start) * 3 * self.natoms))
        self.grad0 = self.grad[:l * 2 * self.nrij]
        grad_shape = (l, self.nrij, end - start, 3 * self.natoms)
        grad_l = self.grad[:l * self.nrij].reshape(grad_shape)
        grad_r = self.grad[l * self.nrij:l * self.nrij * 2].reshape(grad_shape)
        for i, image in enumerate(range(start, end)):
            dmid1 = self.dwdR_mid[image - 1] / 2
            dmid2 = self.dwdR_mid[image] / 2
            grad_l[i + 1, :, i, :] = dmid2 - self.dwdR[image]
            grad_l[i, :, i, :] = dmid1
            grad_r[i + 1, :, i, :] = -dmid2
            grad_r[i, :, i, :] = self.dwdR[image] - dmid1
        for idx in range((end - start) * 3 * self.natoms):
            self.grad[l * self.nrij * 2 + idx, idx] = friction

It's is not entirely clear to me how self.grad takes on the value of grad_l or grad_r - to my understanding, taking the slice of the array creates a new memory reference, so the gradient would be 0 along the cartesian co-ordinates and only the friction component would be added? I am sure the code works given that paths are interpolated correctly for the example cases I have ran.

[GabrielBram]

Nevermind, figured it out! After printing the output, self.grad was indeed pointing to the same memory address as grad_l and grad_r.

If others happen to be interested in the solution used for a sweeping algorithm only:

        grad = np.zeros((4 * len(self.bond_list) + 3 * (im_idx2 - im_idx1) * len(self.images[im_idx1]), 3 * len(self.images[im_idx1])))

        mid1 = self.mid_dwij_list[im_idx1 - 1] / 2
        mid2 = self.mid_dwij_list[im_idx1] / 2
        grad[(1 * len(self.bond_list)):(2 * len(self.bond_list)), :] = mid2 - self.dwij_list[im_idx1]
        grad[:(1 * len(self.bond_list)), :] = mid1
        grad[(3 * len(self.bond_list)):(4 * len(self.bond_list)), :] = -mid2
        grad[(2 * len(self.bond_list)):(3 * len(self.bond_list)), :] = self.dwij_list[im_idx1] - mid1

Thank you for the repository! It was invaluable for replication.