Molecular representations
This section details some of the main representations used in aqml
.
Many-body representation (MBR) in a nutshell
MBR is perhaps the most popular type of representation, among the many variants in literature.
BAML
[Todo]
Atomic Spectrum of Axilrod-Teller-Muto potential (aSLATM)
Several versions of SLATM are available. The most concise & vectorized
version is obtained through numpy
& itertools
. Take as an
example the generation of aSLATM for an atom with atomic index ia0
.
The prerequisite is to generate all possible types of many-body terms, each being characterized by a tuple of Z’s. I.e., the unique 1-body terms are:
zsu1 = [ (zi,) for zi in np.unique(zs0) ]
where zs0
is the list of nuclear charges of the atoms in the
molecule. Similarly, for 2- & 3-body terms, the associated types are
import itertools as ilt
zsu2 = list(itl.combinations_with_replacement(zsu1, 2))
zsu3 = list(itl.combinations_with_replacement(zsu1, 3))
E.g., for a molecule made up of C
& H
,
zsu1 = [ (1,), (6,) ]
zsu2 = [(1, 1), (1, 6), (6, 6)]
zsu3 = [(1, 1, 1), (1, 1, 6), (1, 6, 6), (6, 6, 6)]
For aSLATM, things are a little different: since we’ve already got the
type of the first (or central) atom, i.e., zs[ia0]
(denoted as
Z_I
), we need only to know the the types of at most 2 remaining
atoms, therefore, for the exemplified molecule cxhy
, 2/3-body terms
(mbs2
/mbs3
) are
mbs2 = [ (1,), (6,) ]
mbs3 = [(1, 1), (1, 6), (6, 6)]
And the complete set of many-body terms associated with atom ia0
are
mbs2 = [ (Z_I, 1), (Z_I, 6) ]
mbs3 = [(Z_I, 1, 1), (Z_I, 1, 6), (Z_I, 6, 6)]
Apparently, the type (Z_I, 1,6)
is equiv. to (Z_J, 6,1)
.
2-body terms
x = np.linspace(r0, rcut, ngx);
xJ = (zs0[ia0]*zs0[jas]) * (const/wdr)*np.exp(-0.5*( x[:, None]-ds_ij[None, :])**2/wdr**2) # of size (ngx, nb)
then sum up the same atomic pair indicated by (Z_I, Z_J):
x1new = np.zeros((len(zsu1), ngx)) # ngx: number of grids along the x-axis
for jz, (zj,) in enumerate(zsu1):
x1new[jz] += np.einsum('xi->x', xJ[:, zj==zs[jas]])
Note that for the case of zs[ia0]=6
, the 2-body term (6,1)
is
present, while the type (1,6)
is not.
3-body terms
For three-body terms, we use as the scaling factor the so-called Axilrod-Teller-Muto potential (resulting from a third-order perturbation correction to the attractive London dispersion interactions, i.e., instantaneous induced dipole-induced dipole):
where \(r_{ij}\) is the distance between atoms i
and j
, and
\gamma_{i}
is the angle between the vectors r_{ij}
and
r_{ik}
. The coefficient E_{0}
is positive and of the order
V\alpha^{3}
, where V
is the ionization energy and \alpha
is
the mean atomic polarizability; the exact value of E_{0}
depends on
the magnitudes of the dipole matrix elements and on the energies of the
p
orbitals ([from wikipedia]
(https://en.wikipedia.org/wiki/Axilrod%E2%80%93Teller_potential)).
To vectorize the 3-body feature generation, we need to define a few arrays first:
iats = np.arange(na) # na: total number of atoms of this molecule
atsr = iats[ np.logical_and(ds0[ia0]<=rcut, iats!=ia0) ]
atsp = np.array(list(itl.combinations(atsr, 2)), dtype=np.int32)
ias = np.array([ia0]*len(atsp), dtype=np.int32)
jas = atsp[:,0]; kas = atsp[:,1]
(i|j|k)as
in particular.
xa = np.linspace(-wa/2.0, np.pi+wa/2., nga)
znew = np.zeros((nab,nab, ngx,ngx,nga))
const = 1./(np.sqrt(2*np.pi))
xJ = (const/wdr) * np.exp(-0.5*( x[:, None]-ds_ij[None, :])**2/wdr**2) # of size (ngx, nb)
yK = (const/wdr) * np.exp(-0.5*( x[:, None]-ds_ik[None, :])**2/wdr**2) # of size (ngx, nb)
aJK = (const/wda) * np.exp(-0.5*(xa[:, None]-a_ijk[None, :])**2/wda**2) # size: (nga, nb)
x3 = * np.einsum('xj,yj,aj->xyaj', xJ, yK, aJK)
where,
a_ijk = ang3[(ias, jas, kas)]
ca_ijk = cos3[(ias, jas, kas)]
ca_jik = cos3[(jas, ias, kas)]
ca_kij = cos3[(kas, ias, jas)]
ds_ij = rs0[(ias, jas)]
ds_ik = rs0[(ias, kas)]
ds_jk = rs0[(jas, kas)]
scale = 1.0 # in the original Axilrod–Teller potential, scale=3.0
prefactor = (zs0[ia0]*zs0[jas]*zs0[kas])
prefactor *= (1. + scale*ca_ijk*ca_jik*ca_kij)/(ds_ij*ds_jk*ds_ik)**3
At last, combine (sum up) the 3-body terms of the same type, as
indicated by (Z_I, Z_J, Z_K)
,
x3new = np.zeros((len(zsu2), ngx,ngx,nga))
for zj,zk in zsu2:
ft = np.logical_or(np.logical_and(zs0[jas]==zj,zs0[kas]==zk), np.logical_and(zs0[jas]==zk,zs0[kas]==zj))
x3new[ zsu2.index((zj,zk)) ] = np.einsum('xyaj->xya', x3[:,:,:, ft])