我正在尝试编写一个小代码，给定分子的一组正常振动模式，将它们转换为局部振动模式。为此，我遵循J. Chem 的程序。物理。130, 084106 (2009)他们引入作为正交变换 ( ) 从正常模式 ( ) 到局部模式 ( )。是一个矩阵，Q 是一个矩阵，其中是分子的数量，k 是正常模式的数量。来自于每个分子都有 x,y,z 坐标（$\tilde{Q}=QU$$U$$Q$$\tilde{Q}$$U$$k\times k$$3n\times k$$n$$3n$$Q$矩阵有时表示为三个索引量，因此可以对笛卡尔坐标求和）。

为了找到进行这种转换的，他们定义了函数：并说转换的正确是最大化这个表达式的那个。 $U$

我已经在 Python 中编写了这个函数及其梯度，并试图scipy.optimize.minimize针对一个小案例对其进行优化，但所有方法都失败，并显示消息“由于精度损失而不一定实现所需的错误。”。

关键问题似乎是保持的正交性，因为如果没有这个约束，我怀疑问题是不确定的。我已尝试按照Cross Validated 上的问题的建议将其包含在渐变中来包含此约束，但它似乎不起作用。$U$

我的问题基本上可以归结为：是否有明显的我遗漏的东西可以改善与scipy.optimize.minimizescipy 使用的方法的收敛性，或者问题只是难以解决，我需要尝试不同的方法？我希望不必重新发明轮子，因此首选已经实施的建议。

如果有帮助，我在下面包含了我正在使用的代码。

import numpy as np
from scipy.optimize import minimize

#Q is the matrix of normal modes (3N x k). 

#U is the desired unitary transform matrix (k x k). Optimize only accepts functions of vectors, so have to flatten U. 

#Pipek-Mezey function/criterion
def PM(U,Q,NAtoms,NModes):
    U=np.reshape(U,(NModes,NModes))
    Qtilde=np.einsum('ij,jk->ik',Q,U)

    #Make a 3-index array of cartesian, number of atoms, number of modes for ease of summing
    Qtilde=np.reshape(Qtilde,(3,NAtoms,NModes))
    ##A**2 squares all elements of array, not the square of a matrix
    ##i...->... sums over first index
    C=np.einsum('i...->...',Qtilde**2)

    result=-np.einsum('ij->',C**2)

    return result  

#Gradient of PM condition, shares args with PM. Gradient includes additional factors to 
#ensure U remains orthogonal 
def PM_grad(U,Q,NAtoms,NModes):
    gradient=np.zeros((NModes,NModes))
    U=np.reshape(U,(NModes,NModes))

    #Taking derivatives with respect to U_rs. Sum over s first to reuse intermediates
    for s in range(NModes):
        #Each term of the gradient only contains a column of Qtilde
        Qtildecol=np.einsum('ij,j->i',Q,U[:,s])
        Qtildecol=np.reshape(Qtildecol,(3,NAtoms))
        #Similarly, only a column of C
        Ccol= np.einsum('i...->...',Qtildecol**2)    
        for r in range(NModes):
            Qcol=np.reshape(Q[:,r],(3,NAtoms))
            Intermediate=np.einsum('ij,ij->j',Qcol,Qtildecol)
            gradient[r,s]=-4*np.einsum('i,i->',Intermediate,Ccol)    

    #Gradient correction to ensure orthogonality
    gradient=gradient-(U@np.transpose(gradient)@U) 

    return np.reshape(gradient,(NModes**2))

NModes=3
NAtoms=3
#U is the desired unitary transform matrix (k x k). Optimize only accepts functions of vectors, so have to flatten U. 
U=[0.0,-.800,-.6,.8,-.36,.48,.6,.48,-.64]
#Q is the matrix of normal modes (3N x k).
Q=[0.00000,0.00000,0.00000,0.00000,0.00000,-0.06812,-0.06867,0.05247,0.00000,0.00000,0.00000,0.00000,0.44800,
0.57030,0.54056,0.54492,-0.41639,-0.45330,0.00000,0.00000,0.00000,-0.44800,-0.57030,0.54056,0.54492,-0.41639,0.45330]
Q=np.transpose(np.reshape(Q,(3,9)))


result=minimize(PM,U,args=(Q,NAtoms,NModes),method="BFGS",jac=PM_grad,options={'disp': True})
print(result.x)