求助:分子重组能较小但 RMSD 值偏大，这种情况是否合理？ - 计算化学公社

#!/usr/bin/env python3
import os
import sys
import numpy as np
from numpy.typing import NDArray
def K_cross_mat(vec: NDArray[np.double]) -> NDArray[np.double]: # cross matrix of a 3d vector
return np.array([[0., - vec[2], vec[1]], [vec[2], 0., - vec[0]], [- vec[1], vec[0], 0.]], dtype=np.double)
def L_quat_mat(q: NDArray[np.double]) -> NDArray[np.double]: # left-multiply matrix of a quaternion
return np.block([[q[0], - q[1:4][np.newaxis, :]], [q[1:4][:, np.newaxis], q[0] * np.identity(3, dtype=np.double) + \
K_cross_mat(q[1:4])]])
def R_quat_mat(q: NDArray[np.double]) -> NDArray[np.double]: # right-multiply matrix of a quaternion
return np.block([[q[0], - q[1:4][np.newaxis, :]], [q[1:4][:, np.newaxis], q[0] * np.identity(3, dtype=np.double) - \
K_cross_mat(q[1:4])]])
def quat_rot_mat(q: NDArray[np.double]) -> NDArray[np.double]: # 4x4 rotate matrix of a quaternion
return L_quat_mat(q) @ R_quat_mat(q).T
def get_rotation_matrix(coord1: NDArray[np.double], coord2: NDArray[np.double], allow_reflection: bool=False):
assert coord1.shape == coord2.shape and coord1.shape[1] == 3
H = coord1.T @ coord2
U, Sigma, Vt = np.linalg.svd(H)
R = Vt.T @ U.T
if not allow_reflection:
# only half of O(3) is SO(3)
if np.linalg.det(R) < 0.:
Vt[-1, :] = - Vt[-1, :]
R = Vt.T @ U.T
return R
def get_rot_quaternion(coord1: NDArray[np.double], coord2: NDArray[np.double]) -> NDArray[np.double]:
assert coord1.shape == coord2.shape and coord1.shape[1] == 3
H = coord1.T @ coord2
Delta = np.array([H[1, 2] - H[2, 1], H[2, 0] - H[0, 2], H[0, 1] - H[1, 0]], dtype=np.double)
tr_H = np.trace(H)
N = np.block([[tr_H, Delta[np.newaxis, :]], [Delta[:, np.newaxis], H + H.T - tr_H * np.identity(3, dtype=np.double)]])
eig_val, eig_vec = np.linalg.eigh(N)
q = eig_vec[:, -1].copy()
if q[0] < 0.: q = - q # SU(2) has a 2->1 mapping to SO(3)
return q
def rotate_by_rotation_matrix(coord1: NDArray[np.double], coord2: NDArray[np.double], allow_reflection: bool=False):
assert coord1.shape == coord2.shape and coord1.shape[1] == 3
coord1 -= coord1.mean(axis=0)
coord2 -= coord2.mean(axis=0)
R = get_rotation_matrix(coord1, coord2, allow_reflection=allow_reflection)
coord1 @= R.T
def rotate_by_quaternion(coord1: NDArray[np.double], coord2: NDArray[np.double]):
assert coord1.shape == coord2.shape and coord1.shape[1] == 3
coord1 -= coord1.mean(axis=0)
coord2 -= coord2.mean(axis=0)
q = get_rot_quaternion(coord1, coord2)
coord1 = ((quat_rot_mat(q) @ (np.hstack([np.zeros((len(coord1), 1), dtype=np.double), coord1])).T).T)[:, 1:]
def align_molecules(ifilename1: str, ifilename2: str, imethod: int=2, allow_reflection: bool=False):
# 1 for Kabsch, 2 for quaternion
if imethod == 1:
align_molecules_Kabsch(ifilename1, ifilename2, allow_reflection=allow_reflection)
elif imethod == 2:
if allow_reflection:
raise ValueError("Quaternion does not allow reflection.")
align_moleculs_quaternion(ifilename1, ifilename2)
else:
raise ValueError(f"Unsupported imethod: {imethod:d}")
def align_molecules_Kabsch(ifilename1: str, ifilename2: str, allow_reflection: bool=False):
ifilename1_split = os.path.splitext(ifilename1)
ifilename2_split = os.path.splitext(ifilename2)
assert ifilename1_split[1] == ".xyz" and ifilename2_split[1] == ".xyz"
elements = np.loadtxt(ifilename1, unpack=True, usecols=0, skiprows=2, dtype=np.dtype("<U2"))
coord1 = np.loadtxt(ifilename1, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
coord2 = np.loadtxt(ifilename2, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
rotate_by_rotation_matrix(coord1, coord2, allow_reflection=allow_reflection)
ofilename1 = ifilename1_split[0] + "_aligned.xyz"
ofilename2 = ifilename2_split[0] + "_aligned.xyz"
with open(ofilename1, "w") as ofile1, open(ofilename2, "w") as ofile2:
print(len(elements), file=ofile1)
print(len(elements), file=ofile2)
print(f"aligned from {ifilename1:s}", file=ofile1)
print(f"aligned from {ifilename2:s}", file=ofile2)
for iatom in range(len(elements)):
print(f" {elements[iatom]:<2s}" + (" {:13.8f}" * 3).format(* coord1[iatom]), file=ofile1)
print(f" {elements[iatom]:<2s}" + (" {:13.8f}" * 3).format(* coord2[iatom]), file=ofile2)
def align_molecules_quaternion(ifilename1: str, ifilename2: str):
ifilename1_split = os.path.splitext(ifilename1)
ifilename2_split = os.path.splitext(ifilename2)
assert ifilename1_split[1] == ".xyz" and ifilename2_split[1] == ".xyz"
elements = np.loadtxt(ifilename1, unpack=True, usecols=0, skiprows=2, dtype=np.dtype("<U2"))
coord1 = np.loadtxt(ifilename1, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
coord2 = np.loadtxt(ifilename2, usecols=(1, 2, 3), skiprows=2, dtype=np.double)
rotate_by_quaternion(coord1, coord2)
ofilename1 = ifilename1_split[0] + "_aligned.xyz"
ofilename2 = ifilename2_split[0] + "_aligned.xyz"
with open(ofilename1, "w") as ofile1, open(ofilename2, "w") as ofile2:
print(len(elements), file=ofile1)
print(len(elements), file=ofile2)
print(f"aligned from {ifilename1:s}", file=ofile1)
print(f"aligned from {ifilename2:s}", file=ofile2)
for iatom in range(len(elements)):
print(f" {elements[iatom]:<2s}" + (" {:13.8f}" * 3).format(* coord1[iatom]), file=ofile1)
print(f" {elements[iatom]:<2s}" + (" {:13.8f}" * 3).format(* coord2[iatom]), file=ofile2)
if __name__ == "__main__":
ifilename1 = sys.argv[1]
ifilename2 = sys.argv[2]
# ifilename1 = "1.xyz"
# ifilename2 = "2.xyz"
# align_molecules(ifilename1, ifilename2, imethod=1, allow_reflection=True)
align_molecules(ifilename1, ifilename2)

复制代码