Evaluating Gradient Matrix Symbolically

Given some cartesian coordinates, I wish to define an Euclidean distance matrix (D), and determine the corrosponding gradient matrix wrt all coordinates.

The following code allows me to define a symbolic matrix B with dimensions 9x9, the analytical expressions for the gradient matrix are as expected, but I

am struggling to evaluate the analytical expressions numerically.

Q = dlmread('geom.xyz')
[natom, nc] = size(Q);  % Q = some molecular geometry in cartesian coordinates: 
                                                        %0         0   -0.5053
                                                        %0    1.0392    0.0947
                                                        %0   -1.0392    0.0947

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); % Numerical version of B to store evaluated gradients
[Nsq, tN] = size(B); 


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

Nabla = func() % some function to generate symb expressions
% Nabla has a value of: [ x1, y1, z1, x2, y2, z2, x3, y3, z3]
R = func() % function to gen symb expressions again
% R is a symb expression for cartesian coordinates, with value: [ x1, y1, z1]
                                                               %[ x2, y2, z2]
                                                               %[ x3, y3, z3]
                                                               
count = 0;
for i=1:natom
    for j=1:natom
        if i == j
            Dij = 0; % avoid dividing by 0
        else
            Dij = 1/norm(R(i,:)-R(j,:)); % symbolic expression for 1/distance generated from R
        end
        count = count + 1;
        B(count,:) = gradient(Dij,Nabla); % take gradient for each element of distance matrix, each element gives 9 derivatives when taken wrt to all cartesian coords in nabla.
        for k=1:tN
            BB(count, k) = subs(B(count,k),Dij,1/(norm(Q(i,:)-Q(j,:)))); % gradient matrix of evaluated analytical expressions with dimension natom^2, 3*natom
                                                                     % each row corrosponds to one Euclidean distance of 2 atoms - i.e. elements of D(:)
        end
    end
end

Taking row B(2,:) (corresponding to gradient of D(1,2), where D(1,2) = 1/(abs(x1 - x2)^2 + abs(y1 - y2)^2 + abs(z1 - z2)^2)^(1/2)), the analytical expressions are correctly defined as,

>> B(2,:)'
 
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  

But when evaluated numerically it throws the following error trying to fill matrix BB:

The following error occurred converting from sym to double: Unable to convert expression into double array.

When inspecting the result outside of the loop, the expression I get is,

>> 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   
                                                                                                                                     

Why can I not seem to substitute the numerator of derivates too? How do I substitute these values correctly? These numerical substitution seems off and I can not figure out how to sub this correctly. If anyone has any insight that would be great. Thanks!

Answer 1

[Kyle]

Solved - the numerical substitution must be split into two steps, so that instead of,

for k=1:tN
BB(count, k) = subs(B(count,k),Dij,1/(norm(Q(i,:)-Q(j,:))));
end

we have,

        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

once the values of 1/Euclidean distance have been substituted, one can sub for the values of (abs(x1 - x2)*conj(sign(x1 - x2))). The index for the axis of each vector R must be reset every three iterations in order to loop through the values of x,y and z three times for each of the nine derivatives in the row of B.