基于Tensorflow方便地计算Wilson’s B矩阵及其广义逆进行坐标系转换
Calculating Wilson’s B Matrix and Its Generalized Inverse Conveniently Based on Tensorflow for Coordinate Systems Conversion
Jiawei Xu
Released: 2022-06-19 / Updated: 2022-06-19
1. Wilson’s B矩阵的定义
笛卡尔坐标系具有直观、易于数学处理、定义唯一等优点，而更具有化学意义的坐标系是内坐标或冗余内坐标(redundant internal coordinates, RIC)，合理定义的冗余内坐标在几何优化等问题上具有明显的优势。因此，在以上两类坐标系间进行数据转换是量化计算中不可避免的问题。Wilson’s B矩阵给出了笛卡尔坐标系和冗余内坐标系之间的转换形式，其矩阵元有以下形式：
\begin{equation}
B_{ij}=\frac{\partial q_i}{\partial x_j}
\label{eq1}
\end{equation}
其中q_i和x_j分别为冗余内坐标和笛卡尔坐标。当且仅当冗余内坐标的个数为3N时，Wilson’s B矩阵为一个方阵。由定义可知，Wilson’s B矩阵实际上就是由两组坐标变量给定的Jacobi矩阵。一般来说，对于坐标微元有以下转换成立：
\begin{equation}
\mathbf B\delta \mathbf x=\delta \mathbf q
\label{eq2}
\end{equation}
\begin{equation}
\mathbf B^+\delta \mathbf q=\delta \mathbf x
\label{eq3}
\end{equation}
对于坐标系中的矢量：
\begin{equation}
\mathbf g_x=\mathbf B^{\rm T}\mathbf g_q
\label{eq4}
\end{equation}
\begin{equation}
\mathbf g_q=(\mathbf B^{\rm T})^+\mathbf g_x
\label{eq5}
\end{equation}
Wilson’s B矩阵的形式很简洁，但各元素的求解十分麻烦，尤其是二面角对12个笛卡尔坐标的解析偏导数形式极为复杂，具体数学表达式可参考J. Chem. Phys. 2002, 117, 9160。即使是用现成的解析形式编程，也需要花费不少时间写成并检查代码。下面将介绍Tensorflow自动求导功能的基本用法，并用Tensorflow结合Python求解Wilson’s B矩阵，进一步实现结构、矢量等数据在不同坐标系之间的转换。
2. Tensorflow包的安装和自动求导功能的简介
由于下面将在Python中调用Tensorflow，故假定电脑上已经安装了Python 3.0及以上版本。我们使用pip安装Tensorflow，这是最为便捷的方法。首先检查pip是否是最新版本，如不是则更新至最新版本，然后直接安装即可。
pip install --upgrade pip
pip install tensorflow==2.6.0
如果速度过慢，可以尝试清华大学的镜像。
pip install -i https://pypi.tuna.tsinghua.edu.cn/simple tensorflow==2.6.0
下面简单介绍一下Tensorflow的自动求导功能。Tensorflow自动求导的原理很简单，是中学数学中就学过的链式法则。Tensorflow在整个计算过程中构建一张有向图，在每个计算节点上求当前计算式的导数，完成构图后根据链式法则就可以求出整个计算过程的解析导数。虽然不能直接给出解析形式，但其结果必然是和解析导数一致的。Tensorflow基于有向图的框架使得它成为了图神经网络(graph neural network, GNN)算法中最重要的工具，近几年在深度学习等领域应用十分广泛。Tensorflow自动求导的基本语句如下：
import tensorflow as tf
x=tf.Variable(initial_value=3.0)
with tf.GradientTape() as tape:
y=tf.square(x)
grad=tape.gradient(y,x)
print(grad)
# tf.Tensor(6.0, shape=(), dtype=float32)
首先用tf.Variable声明被求导的变量，并对其赋值initial_value=3.0。然后用with tf.GradientTape() as tape:一行引导一个计算模块，在该模块中的所有计算过程会被用于构图，并记录为tape。对于上面的这个例子，模块中计算了y=x^2。在模块外进行求导，tape.gradient(y,x)用于指定求模块中的函数y对变量x的导数。程序输出结果6.0，并给出变量形状和类型。
也可以声明多个变量求多元函数的偏导数，此时一定要写persistent=True：
import tensorflow as tf
x=tf.Variable(initial_value=2.0)
y=tf.Variable(initial_value=3.0)
with tf.GradientTape(persistent=True) as tape:
z=tf.square(x)+tf.square(y)
grad1=tape.gradient(z,x)
grad2=tape.gradient(z,y)
grad3=tape.gradient(x,y)
print(grad1,grad2,grad3,sep='\n')
# tf.Tensor(4.0, shape=(), dtype=float32)
# tf.Tensor(6.0, shape=(), dtype=float32)
# None (No output since dx/dy does not exist!)
3. 代码实现
3.1 产生内坐标系
为了方便比较不同结构，需要有一套唯一确定的内坐标产生方法。对于更一般的冗余内坐标的情况（当冗余内坐标没有冗余自由度时，即退化为内坐标），也一样需要规定唯一确定的产生方法。这里我们规定，对结构中序号相邻的两个原子定义键长，相邻的三个原子定义键角，相邻的四个原子定义二面角。故一共有(k-1)个键长、(k-2)个键角和(k-3)个二面角。代码实现如下：
def genintval(geom,radunit=False):
bond,angle,dihedral=[],[],[]
N=len(geom)/3
for i in range(N):
if i==0:
pass
else:
coord1=[geom[3*i-3],geom[3*i-2],geom[3*i-1]]
coord2=[geom[3*i],geom[3*i+1],geom[3*i+2]]
bond.append(distance(coord1,coord2))
for i in range(N):
if i<=1:
pass
else:
coord1=[geom[3*i-6],geom[3*i-5],geom[3*i-4]]
coord2=[geom[3*i-3],geom[3*i-2],geom[3*i-1]]
coord3=[geom[3*i],geom[3*i+1],geom[3*i+2]]
angle.append(angle(coord1,coord2,coord3,radunit))
for i in range(N):
if i<=2:
pass
else:
coord1=[geom[3*i-9],geom[3*i-8],geom[3*i-7]]
coord2=[geom[3*i-6],geom[3*i-5],geom[3*i-4]]
coord3=[geom[3*i-3],geom[3*i-2],geom[3*i-1]]
coord4=[geom[3*i],geom[3*i+1],geom[3*i+2]]
dihedral.append(dihedral(coord1,coord2,coord3,coord4,radunit))
intcoord=bond+angle+dihedral
return intcoord,bond,angle,dihedral
该函数的输入geom为笛卡尔坐标系下的结构，可以方便地从各种量化程序的输出文件中读取，排成list型。这里涉及到了内坐标的求算。键长键角是比较方便的，二面角的一种可行的求算方法如下，radunit用于控制计算结果使用何种单位制。经典的半平面法向量求二面角的关键问题在于判断是锐角还是钝角，这里程序实现的替代方法是先求半平面上两点对交线垂足到这两个点的向量，再求向量夹角，绕开了这一问题。
def dotprod(vec1,vec2):
dotprod=0.0
for i in range(len(vec1)):
dotprod+=(vec1[i]*vec2[i])
return dotprod
def norm(vec):
norm=0.0
for i in range(len(vec)):
norm+=(vec[i])**2
norm=math.sqrt(norm)
return norm
def vecdiff(vec1,vec2):
vecdiff=[]
for i in range(len(vec1)):
vecdiff.append(vec1[i]-vec2[i])
return vecdiff
def vecangle(vec1,vec2,radunit=False):
dp=dotprod(vec1,vec2)
norm1=norm(vec1)
norm2=norm(vec2)
angle=math.acos(dp/(norm1*norm2))
if radunit:
pass
else:
angle=180.0*angle/math.pi
return angle
def dotpreptoline(coord1,coord2,coord3):
vec12=vecdiff(coord1,coord2)
vec23=vecdiff(coord2,coord3)
coeffvec1=[vec23[1],-vec23[0],0.0]
coeffvec2=[0.0,-vec23[2],vec23[1]]
coeffvec3=[-vec23[0],-vec23[1],-vec23[2]]
A=numpy.array([coeffvec1,coeffvec2,coeffvec3])
b=numpy.array([coord3[0]*vec23[1]-coord3[1]*vec23[0],coord3[2]*vec23[1]-coord3[1]*vec23[2],-coord1[0]*vec23[0]-coord1[1]*vec23[1]-coord1[2]*vec23[2]])
foot1at23=numpy.linalg.solve(A,b)
distance1to23=distance(foot1at23,coord1)
return distance1to23,foot1at23
def dihedral(coord1,coord2,coord3,coord4,radunit=False):
dis1to23,foot1at23=dotpreptoline(coord1,coord2,coord3)
dis4to23,foot4at23=dotpreptoline(coord4,coord2,coord3)
dihedral=vecangle(vecdiff(coord1,foot1at23),vecdiff(coord4,foot4at23),radunit)
return dihedral
3.2 Tensorflow自动求导计算Wilson’s B矩阵元
下面的函数使用Tensorflow跟踪计算内坐标数值时的数据流计算Wilson’s B矩阵元，在同一个函数里计算内坐标和Wilson’s B矩阵。与直接利用解析式求解的传统方法相比，Tensorflow在代码上明显更为简洁和直观，且能够支持GPU加速，对于处理大尺度体系，或者处理复杂数据（如Hessian矩阵等）的坐标系转换时有显著的优势。
def icm(geom,radunit=True):
dim=len(geom)
N=int(dim/3)
bond12,angle123,dihedral1234=[],[],[]
for i in range(len(geom)):
geom[i]=geom[i]/bohrradius
icm=numpy.zeros((dim-6,dim))
for i in range(N-1):
coord1,coord2=[],[]
for j in range(3):
coord1.append(tf.Variable(initial_value=geom[3*i+j]))
coord2.append(tf.Variable(initial_value=geom[3*i+j+3]))
with tf.GradientTape(persistent=True) as bondcalc:
R12=distance(coord1,coord2)
bond12.append(R12)
for j in range(3):
icm[i,3*i+j]=bondcalc.gradient(R12,coord1[j])
icm[i,3*i+j+3]=bondcalc.gradient(R12,coord2[j])
if i<N-2:
coord3=[]
for j in range(3):
coord3.append(tf.Variable(initial_value=geom[3*i+j+6]))
with tf.GradientTape(persistent=True) as anglecalc:
A123=angle(coord1,coord2,coord3,radunit)
angle123.append(A123)
for j in range(3):
icm[N-1+i,3*i+j]=anglecalc.gradient(A123,coord1[j])
icm[N-1+i,3*i+j+3]=anglecalc.gradient(A123,coord2[j])
icm[N-1+i,3*i+j+6]=anglecalc.gradient(A123,coord3[j])
if i<N-3:
coord4=[]
for j in range(3):
coord4.append(tf.Variable(initial_value=geom[3*i+j+9]))
with tf.GradientTape(persistent=True) as dihedralcalc:
D1234=dihedral(coord1,coord2,coord3,coord4,radunit)
dihedral1234.append(D1234)
for j in range(3):
icm[2*N-3+i,3*i+j]=dihedralcalc.gradient(D1234,coord1[j])
icm[2*N-3+i,3*i+j+3]=dihedralcalc.gradient(D1234,coord2[j])
icm[2*N-3+i,3*i+j+6]=dihedralcalc.gradient(D1234,coord3[j])
icm[2*N-3+i,3*i+j+9]=dihedralcalc.gradient(D1234,coord4[j])
intcoord=bond12+angle123+dihedral1234
return icm,intcoord
3.3 奇异值分解法求Wilson’s B矩阵的Penros广义逆
不满秩的矩阵不可逆，长方矩阵也无法使用方阵求逆的方法。由于Wilson’s B矩阵是不满秩的，且不一定是方阵，即冗余内坐标定义的个数不一定等于3N，故此处使用广义逆(generalized inverse)的概念给出广义逆矩阵\mathbf A^+。广义逆的求法有好几种，这里使用奇异值分解法(singular value decomposition, SVD)。奇异值分解是一种因子分解运算, 将一个矩阵分解为3个矩阵的乘积。其中, 奇异值矩阵\mathbf \Sigma是对角阵，\mathbf U和\mathbf V是两个正交阵。
\begin{equation}
\mathbf A=\mathbf U\mathbf \Sigma \mathbf V
\label{eq6}
\end{equation}
这样，\mathbf A的Penros广义逆可以由\eqref{eq7}给出：
\begin{equation}
\mathbf A^+=\mathbf V^{\rm T}\mathbf \Sigma^{-1} \mathbf U^{\rm T}
\label{eq7}
\end{equation}
对于长方矩阵，通过增加元素全为0的行或列转换为方阵后，可以依照上述方法求广义逆，这在Numpy中有现成的函数。
import numpy as np
pinv=np.linalg.pinv(mat)
求得Wilson’s B矩阵的广义逆或转置广义逆后，即可根据\eqref{eq2}~\eqref{eq5}进行坐标变换。
* Xuliang Wang from Institute of Computing Technology, Chinese Academy of Sciences provided valuable advice for this work.