MATLAB Examples

bond_angle_type.m

  • This function tries to find all bonds and angles between the atom1 and atom2 structures. One optional argument like 'reverse' reverses the angle_limit from max to min angle
  • Tested 15/04/2017
  • Please report bugs to michael.holmboe@umu.se

Contents

Examples

  • atom1 = bond_angle_type(atom1,atom2,Box_dim,1.25,2.25,150)
  • atom1 = bond_angle_type(atom1,atom2,Box_dim,1.25,2.25,150,1)
  • atom1 = bond_angle_type(atom1,atom2,Box_dim,1.25,2.25,150,1,'reverse')
function atom1 = bond_angle_type(atom1,atom2,Box_dim,rmin,rmax,angle_limit,varargin)
% rmin=0;
% rmax=3.5;
% min_angle=120;
% Box_dim=Box_dim(1,:);

if nargin > 5
    skip_internal=cell2mat(varargin(1));
else
    skip_internal=0;
end

lx=Box_dim(1);ly=Box_dim(2);lz=Box_dim(3);
if numel(Box_dim)==3;
    xy=0;xz=0;yz=0;
elseif numel(Box_dim)==9;
    % Box_dim=[lx ly lz 0 0 xy 0 xz yz];
    xy=Box_dim(6);xz=Box_dim(8);yz=Box_dim(9);
end

for i=1:size(atom1,2)
    XYZ_data=[[atom2.x]' [atom2.y]' [atom2.z]'];

    if sum(strncmpi([atom1.type],'OW',2))>0 && sum(strncmpi([atom2.type],'HW',2))>0 && skip_internal == 1;
        ind_sel=~ismember([atom2.index],[atom1(i).index+1 atom1(i).index+2]);
        XYZ_data=XYZ_data(ind_sel,:);
    end

    solute_index=atom1(i).index;
    XYZ_solute=[[atom1(i).x]' [atom1(i).y]' [atom1(i).z]'];

    rx=zeros(size(XYZ_data,1),1);ry=zeros(size(XYZ_data,1),1);rz=zeros(size(XYZ_data,1),1);

    if size(Box_dim,2)>3;
        %Calculate Distance Components for triclic cell with pbc (untested)
        rx =  XYZ_data(:,1) - XYZ_solute(1);
        x_gt_ind=find(rx > lx/2);x_lt_ind=find(rx < - lx/2);
        rx(x_gt_ind) = rx(x_gt_ind) - lx;
        rx(x_lt_ind) = rx(x_lt_ind) + lx;

        ry = XYZ_data(:,2) - XYZ_solute(2);
        y_gt_ind=find(ry > ly/2);y_lt_ind=find(ry < - ly/2);
        ry(y_gt_ind) = ry(y_gt_ind) - ly;
        ry(y_lt_ind) = ry(y_lt_ind) + ly;
        rx(y_gt_ind) = rx(y_gt_ind) - xy;
        rx(y_lt_ind) = rx(y_lt_ind) + xy;

        rz = XYZ_data(:,3) - XYZ_solute(3);
        z_gt_ind=find(rz > lz/2);z_lt_ind=find(rz < - lz/2);
        rz(z_gt_ind) = rz(z_gt_ind) - lz;
        rz(z_lt_ind) = rz(z_lt_ind) + lz;
        rx(z_gt_ind) = rx(z_gt_ind) - xz;
        rx(z_lt_ind) = rx(z_lt_ind) + xz;
        ry(z_gt_ind) = ry(z_gt_ind) - yz;
        ry(z_lt_ind) = ry(z_lt_ind) + yz;
    else
        %Calculate Distance Components for orthogonal cell with pbc
        rx =  XYZ_data(:,1) - XYZ_solute(1);
        x_gt_ind=find(rx > lx/2);x_lt_ind=find(rx < - lx/2);
        rx(x_gt_ind) = rx(x_gt_ind) - lx;
        rx(x_lt_ind) = rx(x_lt_ind) + lx;

        ry = XYZ_data(:,2) - XYZ_solute(2);
        y_gt_ind=find(ry > ly/2);y_lt_ind=find(ry < - ly/2);
        ry(y_gt_ind) = ry(y_gt_ind) - ly;
        ry(y_lt_ind) = ry(y_lt_ind) + ly;

        rz = XYZ_data(:,3) - XYZ_solute(3);
        z_gt_ind=find(rz > lz/2);z_lt_ind=find(rz < - lz/2);
        rz(z_gt_ind) = rz(z_gt_ind) - lz;
        rz(z_lt_ind) = rz(z_lt_ind) + lz;
    end

    dist = sqrt( rx(:,1).^2 + ry(:,1).^2 + rz(:,1).^2 ); % distance calc.

    % Find points inside circle
    in=intersect(find(dist>rmin),find(dist<rmax));

    %in=sort(unique(in));


    %     in = in(find(in~=i));
    neigh.in = in;
    neigh.dist = [dist(in,1)];
    neigh.coords = [XYZ_data(in,1) XYZ_data(in,2) XYZ_data(in,3)];
    neigh.r_vec = [rx(in) ry(in) rz(in)];

    neigh_dist=[neigh.dist];
    neigh_ind=[neigh.in];
    neigh_vec=[neigh.r_vec];

    neigh_ind(~any(neigh_ind,2),:) = [];
    neigh_vec(~any(neigh_vec,2),:) = [];

    a=1;angle_index=zeros(1,12);
    for v=1:size(neigh_ind,1);
        for w=1:size(neigh_ind,1); % From v or from 1?
            angle=rad2deg(atan2(norm(cross(neigh_vec(v,:),neigh_vec(w,:))),dot(neigh_vec(v,:),neigh_vec(w,:))));
            if nargin > 7;
                angle_lo=angle_limit;
                angle_hi=angle;
            else
                angle_hi=angle_limit;
                angle_lo=angle;
            end
            if angle > 0 && angle_lo < angle_hi;
                if neigh_ind(v,1) < neigh_ind(w,1);
                    angle_index(a,1)= atom2(neigh_ind(v,1)).index;
                    angle_index(a,2)= solute_index;
                    angle_index(a,3)= atom2(neigh_ind(w,1)).index;
                    angle_index(a,4)= angle;
                    angle_index(a,5)= neigh_dist(v);
                    angle_index(a,6)= neigh_dist(w);
                    angle_index(a,7:9)= neigh_vec(v,:);
                    angle_index(a,10:12)= neigh_vec(w,:);
                    a=a+1;
                else
                    angle_index(a,1)= atom2(neigh_ind(w,1)).index;
                    angle_index(a,2)= solute_index;
                    angle_index(a,3)= atom2(neigh_ind(v,1)).index;
                    angle_index(a,4)= angle;
                    angle_index(a,5)= neigh_dist(w);
                    angle_index(a,6)= neigh_dist(v);
                    angle_index(a,7:9)= neigh_vec(w,:);
                    angle_index(a,10:12)= neigh_vec(v,:);
                    a=a+1;
                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),:) = [];

    neigh.angle=angle_index;

    atom1(i).neigh.type=[atom2(neigh.in).type]';
    atom1(i).neigh.index=[atom2(neigh.in).index]';
    atom1(i).neigh.dist=neigh.dist;
    atom1(i).neigh.coords=neigh.coords;
    atom1(i).neigh.vec=neigh.r_vec;
    atom1(i).neigh.angle=neigh.angle;

end

% assignin('caller','Neigh',neigh)
% assignin('caller','Angle_index',angle_index)
end