Linear bond angle in internal coordinates

[simonaxelrod]

Hi xTB team,

Thanks so much for the great method and the wonderful open-source code. I really enjoy using it, and my research has benefitted greatly from your tools and documentation.

One thing I'm interested in is how a normal optimization deals with linear bond angles. I know that the derivative of the angle with respect to Cartesian coordinates diverges as the angle approaches 180 degrees. Yet I've never run into this issue in xTB, even when adding a constraining potential to force one angle to be 180 degrees. I'd love to know how this issue is dealt with. The reason is that, when doing optimizations in ASE [using a non-xTB model for the forces], the only option is to use Cartesian coordinates, which requires a huge number of cycles compared to internals. It would be great if I could write my own optimizer in internal coordinates, without the issue of diverging derivatives for 180 degree angles.

Thanks so much!

Simon

[awvwgk]

Has been a while since I looked into the constraining potentials implementation. Not sure whether we actually do anything special. You can checkout the implementation details here:

xtb/src/constrain_pot.f90

Lines 238 to 260 in 92252bb

	vab = va-vb
	vcb = vc-vb
	rab2 = vab(1)*vab(1) + vab(2)*vab(2) + vab(3)*vab(3)
	rcb2 = vcb(1)*vcb(1) + vcb(2)*vcb(2) + vcb(3)*vcb(3)
	call crprod(vcb,vab,vp)
	rp = norm2(vp)+1.d-14
	call impsc(vab,vcb,cosa)
	cosa = dble(min(1.0d0,max(-1.0d0,cosa)))
	theta= dacos(cosa)
	dt = theta - c0
	ea = fix%fc* dt**2
	deddt = 2.d0 * fix%fc * dt
	e = e + ea
	call crprod(vab,vp,deda)
	rmul1 = -deddt / (rab2*rp)
	deda = deda*rmul1
	call crprod(vcb,vp,dedc)
	rmul2 = deddt / (rcb2*rp)
	dedc = dedc*rmul2
	dedb = deda+dedc
	g(:,i) = g(:,i) + deda
	g(:,j) = g(:,j) - dedb
	g(:,k) = g(:,k) + dedc

[simonaxelrod]

Thanks for pointing me to this code! I whipped up a quick version in Python and tried it out for an angle of 180 degrees and a target angle of 179. I also copied over the main code used in ASE for getting derivatives of angles. I get small numbers for the xTB approach (order 1e-2), but huge numbers for ASE (order 1e6, approaching the infinity they try to catch when an angle is 180). I'll definitely raise it as a question/issue in ASE and see what they say. I also wonder if you have any other insights into this?

Here's the code I used:

import numpy as np

def xtb_angle_grad(va, 
                   vb, 
                   vc,
                   k,
                   c0):
    """
    Get the gradient of an angle restoring force 1/2 * k * (angle - c0) ** 2 with
    respect to the atoms A, B, and C that form the angle. Uses the approach in 
    xTB.
    
    va: position of atom A
    vb: position of atom B
    vc: position of atom C
    k: force constant
    c0: reference angle (radians)
    """
    
    vab = va - vb
    vcb = vc - vb

    rab2 = (vab * vab).sum()
    rcb2 = (vcb * vcb).sum()

    vp = np.cross(vcb, vab)
    rp = np.linalg.norm(vp) + 1e-14

    cosa = np.arccos((vab * vcb).sum() / (rab2 * rcb2) ** 0.5)
    cosa = min([1, max(-1, cosa)])
    theta = np.arccos(cosa)

    dt = theta - c0

    ea = 1/2 * k * dt **2
    dedt = k * dt


    deda = np.cross(vab, vp)
    rmul1 = -dedt / (rab2 * rp)
    deda = deda * rmul1


    dedc = np.cross(vcb, vp)
    rmul2 = dedt / (rcb2 * rp)
    dedc = dedc * rmul2

    dedb = deda + dedc

    xtb_de = np.stack([deda, dedb, dedc])
    
    return xtb_de

def ase_angle_grad(va, 
                   vb, 
                   vc,
                   k,
                   c0):

    """
    Same as `xtb_angle_grad`, but using the approach in ASE.
    """
    
    v0 = (va - vb).reshape(1, -1)
    v1 = (vc - vb).reshape(1, -1)

    
    nv0 = np.linalg.norm(v0, axis=-1)
    nv1 = np.linalg.norm(v1, axis=-1)

    v0n = v0 / nv0[:, np.newaxis]
    v1n = v1 / nv1[:, np.newaxis]
    
    # We just normalized the vectors, but in some cases we can get
    # bad things like 1+2e-16.  These we clip away:
    angles = np.arccos(np.einsum('ij,ij->i', v0n, v1n).clip(-1.0, 1.0))
        
    sin_angles = np.sin(angles)
    cos_angles = np.cos(angles)
    
    
    if (sin_angles == 0.).any():  # identify singularities
        raise ZeroDivisionError('Singularity for derivative of a planar angle')

    product = nv0 * nv1
    deriv_d0 = (-(v1 / product[:, np.newaxis]  # derivatives by atom 0
                  - np.einsum('ij,i->ij', v0, cos_angles / nv0**2))
                / sin_angles[:, np.newaxis])
    deriv_d2 = (-(v0 / product[:, np.newaxis]  # derivatives by atom 2
                  - np.einsum('ij,i->ij', v1, cos_angles / nv1**2))
                / sin_angles[:, np.newaxis])
    deriv_d1 = -(deriv_d0 + deriv_d2)  # derivatives by atom 1
    derivs = np.stack((deriv_d0, deriv_d1, deriv_d2), axis=1)
    
    dt = angles - c0
    ase_de = k * dt * (derivs).reshape(-1, 3)
    
    return ase_de

# Three vectors that form a 180 degree angle

va = np.array([-0.40743, -0.18792,  0.30083])
vb = np.array([ 0.38103,  0.36053, -0.45036])
vc = np.array([ 1.21928497,  0.94361721, -1.24899119])

k = 3
c0 = np.radians(179)


xtb_de = xtb_angle_grad(va, vb, vc, k, c0)
ase_de = ase_angle_grad(va, vb, vc, k, c0)

print(xtb_de) 
# [[-0.00172038 -0.03377378 -0.02646425]
#  [-0.00333856 -0.06554128 -0.05135643]
#  [-0.00161818 -0.03176751 -0.02489219]]

print(ase_de)
# [[  -59787.678243   -1173727.51733191  -919702.09500646]
#  [  116023.80097985  2277731.94391443  1784771.02707537]
#  [  -56236.12273686 -1104004.42658252  -865068.93206891]]

(Just to make sure I was copying over the ASE approach properly, I confirmed that using their code directly gives the same results)

# for comparison to make sure we got the ase derivatives the same way they did
from ase.geometry.geometry import get_angles_derivatives, get_angles

v0 = (va - vb).reshape(1, -1)
v1 = (vc - vb).reshape(1, -1)
angle = get_angles(v0, v1)

print((k * np.radians(180 - 179) * np.radians(get_angles_derivatives(v0, v1))).reshape(-1, 3))
# [[  -59787.678243   -1173727.51733192  -919702.09500647]
#  [  116023.80097985  2277731.94391446  1784771.02707539]
#  [  -56236.12273686 -1104004.42658254  -865068.93206892]]