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 9×9, 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 endNabla = 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 endend
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!
Best Answer