象征性地评估梯度矩阵

Question

给定 3 个原子的笛卡尔坐标中的 3x3 距离矩阵，我希望采用梯度 wrt \nabla。我的问题是如下评估梯度矩阵 (B) 的符号元素，

clear all

%Q = dlmread('bentgeom.xyz'); % ref geometry 
Q = [ 0 0 -0.5053; 0 1.0392 0.0947; 0 -1.0392 0.0947]; % cartesian coordinates for tri-atomic
[natom, nc] = size(Q);

for i=1:natom
    for j=i+1:natom
        D(i,j) = norm(Q(i,:)-Q(j,:)); % distance matrix - numerical values
        D(j,i) = 0; % upper triangular
    end                                     
end

d = D(:); % flatten distance matrix

stot=[' '];
for i=1:natom
    A = i;
      s1=sprintf(' x%d, ',A);
      s2=sprintf(' y%d, ',A);
      s3=sprintf(' z%d ',A);
      stot=strcat(stot, s1, s2, s3);
end

Fext1=['[' stot ']']; % defining symbolic nabla (i.e. derivative over all cartesian coordinates for all atoms)
Nabla=str2sym(Fext1);  
        
R(1,:) = Nabla(1:3); % define symbolic cartesian vectors for each of the three atoms
R(2,:) = Nabla(4:6);
R(3,:) = Nabla(7:9);

        
B = sym(zeros(natom^2, 3*natom)); % init symbolic B matrix (derivative of flattened distance mat. wrt cartesian coordinates)
BB = zeros(natom^2, 3*natom); % init gradient matrix for storing numerical values of evaluated analytical elements of symbolic B
[Nsq, tN] = size(B);  % rows of B/BB corrospond to elements of flattened distance matrix (natom^2 rows)
                      % columns of B/BB corrospond to each derivative wrt
                      % nabla - therefore (3*natom cols) - 3/9 are 0
                     

count = 0;
for i=1:natom % Loop over all cartesian coordinates
    for j=1:nc
        if i == j % avoids dividing by 0 for diag elements
            Dij = 0;
        else
            Dij = 1/norm(R(i,:)-R(j,:)); % take 1/Euclidean distance for each off-diag
                                         % element - this is symbolic    
        end
        count = count + 1;
        B(count,:) = gradient(Dij,Nabla); % for each symb distance take gradient wrt nabla and store the 9 derivatives in cols
        for k=1:tN
            BB(count, k) = subs(B(count,k),Dij,1/(norm(Q(i,:)-Q(j,:)))); % evaluate symb matrix B and fill BB
                                                                         % with corrosponding numerical values
        end
    end
end

% G = BB*BB'; Not relevent to the problem
% 
% [eigVec, eigVal] = eig(G);  % Diag matrix & extract eigvals and vecs
% diagVals = diag(eigVal);
% 
% count = 1;
% for i=1:tN  % Extract eigenvectors that have a non-zero eigenvalue 
%     if diagVals(i) > 1.d-3 % non-zero defined by some thresh
%         U(:, count) = eigVec(:, i); % U contains vecs of non-zero eigvals
%         fVal(count) = diagVals(i);  % should be natom (3) no-zero eigvals
%         count = count + 1;
%     else
%         count = count;
%     end
% end
% 
% q = U'*d;  % Transformation of distance matrix into new internal coordinate system

为了清楚起见，nabla 和 R 被定义为，

>> Nabla
 
Nabla =
 
[ x1, y1, z1, x2, y2, z2, x3, y3, z3]
 
>> R
 
R =
 
[ x1, y1, z1]
[ x2, y2, z2]
[ x3, y3, z3]
 

这为我提供了梯度矩阵元素的正确解析表达式，但是我似乎无法正确进行数字替换。取行B(2,:)（对应D(1,2)的梯度，其中D(1,2) = 1/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(1/2)），解析表达式正确定义为，

>> B(2,:)' % transpose for ease of reading, rows of output are in fact cols
 
ans =

 -(abs(x1 - x2)*conj(sign(x1 - x2)))/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(3/2)
 -(abs(y1 - y2)*conj(sign(y1 - y2)))/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(3/2)
 -(abs(z1 - z2)*conj(sign(z1 - z2)))/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(3/2)
  (abs(x1 - x2)*conj(sign(x1 - x2)))/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(3/2)
  (abs(y1 - y2)*conj(sign(y1 - y2)))/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(3/2)
  (abs(z1 - z2)*conj(sign(z1 - z2)))/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(3/2)
                                                                                            0
                                                                                            0
                                                                                            0

但是当用数值计算时，它会在尝试填充矩阵 BB 时抛出以下错误：

从 sym 转换为 double 时出现以下错误： 无法将表达式转换为双精度数组。

在循环外检查结果时，我得到的表达式是，

>> subs(B(2,:),Dij,1/(norm(Q(1,:)-Q(2,:))))'
 
ans =

 -(422888205284047902665468512568529212604620021371*abs(x1 - x2)*conj(sign(x1 - x2)))/730750818665451459101842416358141509827966271488
 -(422888205284047902665468512568529212604620021371*abs(y1 - y2)*conj(sign(y1 - y2)))/730750818665451459101842416358141509827966271488
 -(422888205284047902665468512568529212604620021371*abs(z1 - z2)*conj(sign(z1 - z2)))/730750818665451459101842416358141509827966271488
  (422888205284047902665468512568529212604620021371*abs(x1 - x2)*conj(sign(x1 - x2)))/730750818665451459101842416358141509827966271488
  (422888205284047902665468512568529212604620021371*abs(y1 - y2)*conj(sign(y1 - y2)))/730750818665451459101842416358141509827966271488
  (422888205284047902665468512568529212604620021371*abs(z1 - z2)*conj(sign(z1 - z2)))/730750818665451459101842416358141509827966271488
                                                                                                                                     0
                                                                                                                                     0
                                                                                                                                     0

为什么我似乎也不能替换导数的分子？如何正确替换这些值？这些数字替换似乎关闭了，我无法弄清楚如何正确地替换它。如果有人有任何见解，那就太好了。谢谢！

Answer 1

变电站必须分成两个独立的步骤，

        ax = 0;
        for k=1:tN
            if mod(ax,3) == 0
                ax = 0;
            end
            ax = ax + 1;
            b = subs(B(count,k),Dij,1/(norm(Q(i,:)-Q(j,:)))); % evaluate symb matrix B and fill BB
            c = subs(b,(R(i,ax)-R(j,ax)),(Q(i,ax)-Q(j,ax)));
            BB(count,k) = c;
        end

一旦计算出分母，就可以代入每个元素的分子。前三列是第一个原子的导数，在每个 x,y,z 上，第二组三列是导数每个轴上的第二个原子，最后三个是第三个原子。建立了一个符号矩阵，表示每个 at 的笛卡尔坐标，

>> R
 
R =
 
[ x1, y1, z1]
[ x2, y2, z2]
[ x3, y3, z3]

我们必须每 3 次迭代重置计数器 ax，以便将它们循环 3 次，总共 9 次迭代。