MATLAB Examples

cell_list_dist_matrix_atom.m

  • This function calculates the distance matrix from the atom struct,
  • using a cell list algorithm adapted from the Matlab MDtoolbox by
  • Yasuhiro Matsunaga, link to github
  • http://github.com/ymatsunaga/mdtoolbox/

Original reference to the cell list algorithm: Heinz, T.N. & Hunenberger, P.H. "A fast pairlist-construction algorithm for molecular simulations under periodic boundary conditions." J Comput Chem 25, Sep;25(12):1474-1486 (2004).

Contents

Examples

  • dist_matrix = cell_list_dist_matrix_atom(atom,Box_dim)
  • dist_matrix = cell_list_dist_matrix_atom(atom,Box_dim,1.25) % Setting the max distance for bonds with H's
  • dist_matrix = cell_list_dist_matrix_atom(atom,Box_dim,1.25,2.25) % Setting the max distance for bonds, otherwise default is 12 Å
  • dist_matrix = cell_list_dist_matrix_atom(atom,Box_dim,1.25,2.25,2.25) % Setting the max distance for angles, default is not to calculate any angles
function dist_matrix = cell_list_dist_matrix_atom(atom,Box_dim,varargin) % ,atom2,Box_dim); % or % ,Box_dim);
disp('Using a cell list algorithm adapted from the Matlab MDtoolbox by');
disp('Yasuhiro Matsunaga, link to github:');
disp('http://github.com/ymatsunaga/mdtoolbox/');

nAtoms=size(atom,2);
dist_matrix=zeros(nAtoms);
X_dist=zeros(nAtoms); % These matrixes are curently not symmetric. A=max(abs(X_dist),abs(X_dist')) makes everything +positive
Y_dist=zeros(nAtoms);
Z_dist=zeros(nAtoms);

H_ind=find(strncmpi([atom.type],'H',1));
if nargin > 2
    rHmax=varargin{1};
else
    rHmax=1.25;
end

if nargin >3
    calc_bond_list=1;
    rmax=varargin{2};
else
    calc_bond_list=0;
    rmax=12;
end

if nargin > 4
    calc_angle_list=1;
    rmax_angle=varargin{3};
else
    calc_angle_list=0;
    %   rmax_angle=2.25;
end

rmax2=rmax.^2;

if numel(Box_dim)==1
    Box_dim(1)=Box_dim(1);
    Box_dim(2)=Box_dim(1);
    Box_dim(3)=Box_dim(1);
    xy=0;xz=0;yz=0;
elseif numel(Box_dim)==3
    xy=0;xz=0;yz=0;
else
    xy=Box_dim(6);xz=Box_dim(8);yz=Box_dim(9);
end

XYZ_data=[[atom.x]' [atom.y]' [atom.z]'];
if numel(Box_dim)==3
    XYZ_data(:,1) = XYZ_data(:,1) - Box_dim(1)*floor(XYZ_data(:,1)./Box_dim(1));
    XYZ_data(:,2) = XYZ_data(:,2) - Box_dim(2)*floor(XYZ_data(:,2)./Box_dim(2));
    XYZ_data(:,3) = XYZ_data(:,3) - Box_dim(3)*floor(XYZ_data(:,3)./Box_dim(3));
elseif numel(Box_dim)==9 && numel(nonzeros(Box_dim(1,4:end))) > 0 %&& sum(abs(Box_dim([6 8 9]))) > 0
    % Untested
    XYZ_data(:,1) = XYZ_data(:,1) - xy*floor(XYZ_data(:,2)./Box_dim(2));
    XYZ_data(:,1) = XYZ_data(:,1) - xz*floor(XYZ_data(:,3)./Box_dim(3));
    XYZ_data(:,2) = XYZ_data(:,2) - yz*floor(XYZ_data(:,3)./Box_dim(3));
    XYZ_data(:,1) = XYZ_data(:,1) - Box_dim(1)*floor(XYZ_data(:,1)./Box_dim(1));
    XYZ_data(:,2) = XYZ_data(:,2) - Box_dim(2)*floor(XYZ_data(:,2)./Box_dim(2));
    XYZ_data(:,3) = XYZ_data(:,3) - Box_dim(3)*floor(XYZ_data(:,3)./Box_dim(3));
end

if any(Box_dim(1:3) < (2*rmax))
    disp('Decrease rmax!!!')
    atom = dist_matrix_atom(atom,Box_dim);
    pause(5)
    return
end

rcell = max([8.0 8.0 8.0], rmax);
ncell = floor(Box_dim(1:3)./rcell); % [Nx Ny Nz] discretized grid cells
rcell = Box_dim(1:3)./ncell; % dimensions of the discretized grid cells

nx = floor(XYZ_data(:,1)./rcell(1))'; % x positions index
ny = floor(XYZ_data(:,2)./rcell(2))'; % y positions index
nz = floor(XYZ_data(:,3)./rcell(3))'; % z positions index

% Number of grid cells
M = prod(ncell);

% Grid cell indexes
m = ncell(1)*ncell(2)*nz + ncell(1)*ny + nx + 1;

mindex = cell(M, 1);
mindex3 = cell(M, 1);
for i = 1:M
    mindex{i} = find(m == i);
    ind3 = zeros(1, numel(mindex{i})*3);
    ind3(1:3:end) = 3.*(mindex{i}-1)+1;
    ind3(2:3:end) = 3.*(mindex{i}-1)+2;
    ind3(3:3:end) = 3.*(mindex{i}-1)+3;
    mindex3{i} = ind3;%to3(mindex{i});
end

calculate cell mask

dm = 1:(M - 1);
dmx = mod(dm, ncell(1));
dmy = floor(mod(dm, ncell(1)*ncell(2)) ./ ncell(1));
dmz = floor(dm ./ (ncell(1)*ncell(2)));

dnx = abs(minimum_image(dmx, ncell(1)));

dny = zeros(size(dmy));
logical_index = (dmx == 0) | (dmy == (ncell(2)-1) & dmz == (ncell(3)-1));
dny(logical_index) = abs(minimum_image(dmy(logical_index), ncell(2)));
dny(~logical_index) = min(abs(minimum_image(dmy(~logical_index), ncell(2))), abs(minimum_image(dmy(~logical_index) + 1, ncell(2))));

dnz = zeros(size(dmz));
logical_index = (dmz == (ncell(3)-1)) | (dmx == 0 & dmy == 0);
dnz(logical_index) = abs(minimum_image(dmz(logical_index), ncell(3)));
dnz(~logical_index) = min(abs(minimum_image(dmz(~logical_index), ncell(3))), abs(minimum_image(dmz(~logical_index) + 1, ncell(3))));

mask = zeros(size(dm));
mask = (max(dnx, 1) - 1).^2 * rcell(1).^2 + (max(dny, 1) - 1).^2 * rcell(2).^2 + (max(dnz, 1) - 1).^2 * rcell(3).^2 <= rmax2;
%mask = true(1, M - 1); % bug check of mask
mask_index = find(mask);

calculate the distances of the atoms in masked cells

pair = zeros(nAtoms*1000, 2);
dist = zeros(nAtoms*1000, 1);
icount = 1;
for m1 = 1:M
    m1index = mindex{m1};
    if numel(m1index) == 0
        continue
    elseif numel(m1index) > 1
        [lpair,ldist,num,rx,ry,rz] = calcpair_atom_pbc(XYZ_data(m1index,:), rmax, Box_dim); % New
        if num > 0
            X_dist(m1index,m1index)=rx;
            Y_dist(m1index,m1index)=ry;
            Z_dist(m1index,m1index)=rz;
            pair(icount:(icount+num-1), :) = [m1index(lpair(:,1)')' m1index(lpair(:,2)')'];
            dist(icount:(icount+num-1)) = ldist;
            icount = icount + num;
        end
    end

    m2 = m1 + mask_index;
    m2 = m2(m2 <= M);
    if isempty(m2)
        continue
    end
    m2index = [mindex{m2}];
    [lpair,ldist,num,rx,ry,rz] = calcpair2_atom_pbc(XYZ_data(m1index,:), XYZ_data(m2index,:), rmax, Box_dim); % New
    if num > 0
        X_dist(m1index,m2index)=rx;
        Y_dist(m1index,m2index)=ry;
        Z_dist(m1index,m2index)=rz;
        X_dist(m2index,m1index)=-rx';
        Y_dist(m2index,m1index)=-ry';
        Z_dist(m2index,m1index)=-rz';
        pair(icount:(icount+num-1),:) = [m1index(lpair(:,1)')' m2index(lpair(:,2)')'];
        dist(icount:(icount+num-1)) = ldist;
        icount = icount + num;
    end
end

pair(icount:end,:) = [];
dist(icount:end) = [];

sort pair index (upper triangular form)

for i = 1:size(pair,1)
    if pair(i,1) > pair(i,2)
        tmp = pair(i,1);
        pair(i,1) = pair(i,2);
        pair(i,2) = tmp;
    end
end

if calc_bond_list==1
    Bond_index = [pair dist];
    Bond_index(Bond_index(:,3)>rmax,:)=[];
    % Special treatment for the short H-distances
    Hcol1=ismember(Bond_index(:,1),H_ind);
    Hcol2=ismember(Bond_index(:,2),H_ind);
    H_ind=Hcol1+Hcol2;
    H_ind=H_ind>0;
    Bond_index(Bond_index(H_ind,3)>rHmax,:)=[];
    [Y,I] = sort(Bond_index(:,2));
    Bond_index = Bond_index(I,:);
    Bond_index = unique(Bond_index,'rows','stable');
    Bond_index(~any(Bond_index,2),:) = [];
    [Y,I] = sort(Bond_index(:,1));
    Bond_index = Bond_index(I,:);

    for i=1:size(Bond_index,1)
        dist_matrix(Bond_index(i,1),Bond_index(i,2))=Bond_index(i,3);
        dist_matrix(Bond_index(i,2),Bond_index(i,1))=Bond_index(i,3); % Do we need this?
    end
    assignin('caller','Bond_index',Bond_index)
else
    for i=1:size(pair,1)
        dist_matrix(pair(i,1),pair(i,2))=dist(i);
        dist_matrix(pair(i,2),pair(i,1))=dist(i); % Do we need this?
    end

end

if calc_angle_list==1
    a=1;
    Angle_dist_index=Bond_index(Bond_index(:,3)<rmax_angle,:);
    Angle_index=zeros(6*size(XYZ_data,1),10);
    mid_angle_ind=unique(Angle_dist_index(:,1));
    for u=1:length(mid_angle_ind)
        angle_ind=find(Angle_dist_index(:,1)==mid_angle_ind(u));
        neigh_ind=Angle_dist_index(angle_ind,2);
        [pair, dist, num,rx,ry,rz] = calcpair2_atom_pbc(XYZ_data(mid_angle_ind(u),:),XYZ_data(neigh_ind,:), 2.25, Box_dim);
        neigh_vec=[rx' ry' rz'];
        for v=1:size(neigh_ind,1)
            for w=1:size(neigh_ind,1)
                angle=rad2deg(atan2(norm(cross(neigh_vec(v,:),neigh_vec(w,:))),dot(neigh_vec(v,:),neigh_vec(w,:))));
                if angle > 0 && angle <= 180
                    if v < w
                        Angle_index(a,1)= neigh_ind(v,1);
                        Angle_index(a,2)= mid_angle_ind(u);
                        Angle_index(a,3)= neigh_ind(w,1);
                        Angle_index(a,4)= angle;
                        Angle_index(a,5:7)= neigh_vec(v,:);
                        Angle_index(a,8:10)= neigh_vec(w,:);
                        a=a+1;
                    else
                        Angle_index(a,1)= neigh_ind(w,1);
                        Angle_index(a,2)= mid_angle_ind(u);
                        Angle_index(a,3)= neigh_ind(v,1);
                        Angle_index(a,4)= angle;
                        Angle_index(a,5:7)= neigh_vec(w,:);
                        Angle_index(a,8:10)= neigh_vec(v,:);
                        a=a+1;
                    end
                end
            end
        end
    end
    [Y,I]=sort(Angle_index(:,2));
    Angle_index=Angle_index(I,:);
    Angle_index = unique(Angle_index,'rows','stable');
    Angle_index(~any(Angle_index,2),:) = [];
    assignin('caller','Angle_index',Angle_index)
end

assignin('caller','dist_matrix',dist_matrix);
assignin('caller','X_dist',(X_dist)');
assignin('caller','Y_dist',(Y_dist)');
assignin('caller','Z_dist',(Z_dist)');

    function mi_n = minimum_image(n, N)
        mi_n = n - N .* sign(n) .* floor((abs(n) + 0.5*N) ./ N);
    end

    function [pair, dist, num, rx, ry, rz] = calcpair_atom_pbc(XYZ_data, rmax, Box_dim)
        if numel(Box_dim)==3
            rx = bsxfun(@minus,XYZ_data(:,1),XYZ_data(:,1)');
            rx = rx - Box_dim(1)*round(rx./Box_dim(1));
            ry = bsxfun(@minus,XYZ_data(:,2),XYZ_data(:,2)');
            ry = ry - Box_dim(2)*round(ry./Box_dim(2));
            rz = bsxfun(@minus,XYZ_data(:,3),XYZ_data(:,3)');
            rz = rz - Box_dim(3)*round(rz./Box_dim(3));
        elseif numel(Box_dim)==9 && numel(nonzeros(Box_dim(1,4:end))) > 0 % Untested
            rx = bsxfun(@minus,XYZ_data(:,1),XYZ_data(:,1)');
            ry = bsxfun(@minus,XYZ_data(:,2),XYZ_data(:,2)');
            rz = bsxfun(@minus,XYZ_data(:,3),XYZ_data(:,3)');

            rx = rx - xz*round(rz./Box_dim(3));
            ry = ry - yz*round(rz./Box_dim(3));
            rz = rz - Box_dim(3)*round(rz./Box_dim(3));

            rx = rx - xy*round(ry./Box_dim(2));
            ry = ry - Box_dim(2)*round(ry./Box_dim(2));

            rx = rx - Box_dim(1)*round(rx./Box_dim(1));
        end
        dist = sqrt(rx.^2 + ry.^2 + rz.^2);
        index = find(tril(dist < rmax, -1));
        if isempty(index)
            pair = [];
            dist = [];
            num = 0;
            rx = [];
            ry = [];
            rz = [];
        else
            [pair1, pair2] = ind2sub(size(dist), index);
            pair = zeros(numel(pair1), 2);
            pair(:,1) = pair1;
            pair(:,2) = pair2;
            dist = dist(index);
            num = numel(dist);
        end
    end

    function [pair, dist, num, rx, ry, rz] = calcpair2_atom_pbc(XYZ_data1, XYZ_data2, rmax, Box_dim)
        if numel(Box_dim)==3
            rx = bsxfun(@minus,XYZ_data1(:,1),XYZ_data2(:,1)');
            rx = rx - Box_dim(1)*round(rx./Box_dim(1));
            ry = bsxfun(@minus,XYZ_data1(:,2),XYZ_data2(:,2)');
            ry = ry - Box_dim(2)*round(ry./Box_dim(2));
            rz = bsxfun(@minus,XYZ_data1(:,3),XYZ_data2(:,3)');
            rz = rz - Box_dim(3)*round(rz./Box_dim(3));
        elseif numel(Box_dim)==9 && numel(nonzeros(Box_dim(1,4:end))) > 0 % Untested
            rx = bsxfun(@minus,XYZ_data1(:,1),XYZ_data2(:,1)');
            ry = bsxfun(@minus,XYZ_data1(:,2),XYZ_data2(:,2)');
            rz = bsxfun(@minus,XYZ_data1(:,3),XYZ_data2(:,3)');

            rx = rx - xz*round(rz./Box_dim(3));
            ry = ry - yz*round(rz./Box_dim(3));
            rz = rz - Box_dim(3)*round(rz./Box_dim(3));

            rx = rx - xy*round(ry./Box_dim(2));
            ry = ry - Box_dim(2)*round(ry./Box_dim(2));

            rx = rx - Box_dim(1)*round(rx./Box_dim(1));
        end
        dist = sqrt(rx.^2 + ry.^2 + rz.^2);
        index = find(dist < rmax);
        if isempty(index)
            pair = [];
            dist = [];
            num = 0;
            rx = [];
            ry = [];
            rz = [];
        else
            [pair1, pair2] = ind2sub(size(dist), index);
            pair = zeros(numel(pair1), 2);
            pair(:,1) = pair1;
            pair(:,2) = pair2;
            dist = dist(index);
            num = numel(dist);
        end
    end
end % end function