MATLAB Examples

# bond_angle_atom.m

• This function tries to find all bonds and angles of the atom struct
• One optional argument like 'more' will give more bond info
• atom is the atom struct
• Box_dim is the box dimension vector
• Tested 15/04/2017
• Please report bugs to michael.holmboe@umu.se

## Examples

• atom=bond_angle_atom(atom,Box_dim,1.25,2.25)
• atom=bond_angle_atom(atom,Box_dim,1.25,2.25,'more')
```function atom=bond_angle_atom(atom,Box_dim,max_short_dist,max_long_dist,varargin) ```

tic

```nAtoms=size(atom,2); % max_short_dist=1.2; % max_long_dist=2.3; XYZ_labels=[atom.type]'; XYZ_data=[[atom.x]' [atom.y]' [atom.z]']; lx=Box_dim(1);ly=Box_dim(2);lz=Box_dim(3); if length(Box_dim) > 3 xy=Box_dim(6);xz=Box_dim(8);yz=Box_dim(9); else xy=0;xz=0;yz=0; end close_count=1; Bond_index=zeros(4*size(XYZ_data,1),3); Angle_index=zeros(4*size(XYZ_data,1),4); dist_matrix=zeros(nAtoms,nAtoms); % dist_matrix = dist_matrix_atom(atom,Box_dim); % Do this if you want to calc dist_matrix_atom first b=1;a=1; overlap_index=[]; for i = 1:size(XYZ_data,1) ```
``` max_neigh_distance = max_long_dist;%2.3 if length(Box_dim)>3 %Calculate Distance Components for triclic cell rx = XYZ_data(i,1) - XYZ_data(:,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(i,2) - XYZ_data(:,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(i,3) - XYZ_data(:,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 ortogonal cell rx = XYZ_data(i,1) - XYZ_data(:,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(i,2) - XYZ_data(:,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(i,3) - XYZ_data(:,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 r = sqrt( rx(:,1).^2 + ry(:,1).^2 + rz(:,1).^2 ); % distance calc. bond_in=intersect(find(r > 0), find(r < max_neigh_distance)); dist_matrix(:,i)=r; % Else if we already called dist_matrix_atom... % r = dist_matrix(:,i); % rx= X_dist(:,i); % ry= Y_dist(:,i); % rz= Z_dist(:,i); % bond_in=intersect(find(r > 0), find(r < max_neigh_distance)); ```
``` n=1; Neigh_ind=zeros(12,1);Neigh_vec=zeros(12,3); for j=1:length(bond_in); if [atom(i).molid]==[atom(bond_in(j)).molid] % Original, uncomment this section % max_distance = max_short_dist; % Only test if strncmp(strtrim(XYZ_labels(i)),{'H'},1) || strncmp(strtrim(XYZ_labels(bond_in(j))),{'H'},1) max_distance = max_short_dist; elseif strcmpi(strtrim(XYZ_labels(i)),{'Oalhh'}) || strcmp(strtrim(XYZ_labels(i)),{'Osih'}) if max_long_dist < 2.25 max_distance = 2.25; else max_distance = max_long_dist; end %disp('Looking for M-O-H bonds') else max_distance = max_long_dist; end if r(bond_in(j)) < 0.6; disp('Atoms to close!!!') r(bond_in(j)) [i,bond_in(j)] XYZ_labels(i) XYZ_labels(bond_in(j)) XYZ_data(i,:) XYZ_data(bond_in(j),:) overlap_index=[overlap_index; {i bond_in(j) r(bond_in(j)) XYZ_labels(i) XYZ_labels(bond_in(j)) XYZ_data(i,:) XYZ_data(bond_in(j),:)}]; elseif r(bond_in(j)) > max_distance && r(bond_in(j)) < 1.25; disp('Atoms pretty close...') r(bond_in(j)) [i,bond_in(j)] XYZ_labels(i) XYZ_labels(bond_in(j)) XYZ_data(i,:) XYZ_data(bond_in(j),:) close_count=close_count+1; overlap_index=[overlap_index; {i bond_in(j) r(bond_in(j)) XYZ_labels(i) XYZ_labels(bond_in(j)) XYZ_data(i,:) XYZ_data(bond_in(j),:)}]; end if r(bond_in(j)) > max_distance/3 && r(bond_in(j)) < max_distance;%strncmp(XYZ_labels(i),strtrim(XYZ_labels(j)),1) == 0; if strncmpi(strtrim(XYZ_labels(i)),{'OW'},2)% || strncmpi(strtrim(XYZ_labels(bond_in(j))),{'OW'},2); % < bond_in(j) && ismember(XYZ_labels(i),{'Ow','Hw','OW','HW','HW1','HW2'}) > 0; if bond_in(j) > i && bond_in(j) < i+3; %strncmpi(strtrim(XYZ_labels(i)),{'Ow'},2) && bond_in(j) > i && bond_in(j) < i+3; Bond_index(b,1)= min([i bond_in(j)]); Bond_index(b,2)= max([i bond_in(j)]); Bond_index(b,3)= r(bond_in(j)); b=b+1; if r(bond_in(j)) < max_short_dist % This should always be true right? Neigh_ind(n,1)= bond_in(j); Neigh_vec(n,1:3) = [rx(bond_in(j)) ry(bond_in(j)) rz(bond_in(j))]; n=n+1; end end elseif strncmpi(strtrim(XYZ_labels(i)),{'HW'},2) % disp('HW') elseif strncmpi(strtrim(XYZ_labels(i)),{'H'},1) && strncmpi(strtrim(XYZ_labels(bond_in(j))),{'H'},1) % disp('disallowed a H-H bond') elseif ismember({'H'},[XYZ_labels(i) XYZ_labels(bond_in(j))]) > 0 || ismember({'Oalhh'},[XYZ_labels(i) XYZ_labels(bond_in(j))]) > 0 || ismember({'Osih'},[XYZ_labels(i) XYZ_labels(j)]) > 0 && ismember(XYZ_labels(bond_in(j)),{'Ow','Hw','OW','HW','HW1','HW2'}) == 0; Bond_index(b,1)= min([i bond_in(j)]); Bond_index(b,2)= max([i bond_in(j)]); Bond_index(b,3)= r(bond_in(j)); b=b+1; if r(bond_in(j)) < max_distance % This should always be true right? Neigh_ind(n,1)= bond_in(j); Neigh_vec(n,1:3) = [rx(bond_in(j)) ry(bond_in(j)) rz(bond_in(j))]; n=n+1; end elseif ~strncmpi(strtrim(XYZ_labels(i)),{'OW'},2) || ~strncmpi(strtrim(XYZ_labels(i)),{'HW'},2) if bond_in(j) < i; Bond_index(b,1)= min([i bond_in(j)]); Bond_index(b,2)= max([i bond_in(j)]); Bond_index(b,3)= r(bond_in(j)); b=b+1; end if r(bond_in(j)) < max_distance % This should always be true right? Neigh_ind(n,1)= bond_in(j); Neigh_vec(n,1:3) = [rx(bond_in(j)) ry(bond_in(j)) rz(bond_in(j))]; n=n+1; end end end end % test Neigh_ind(~any(Neigh_ind,2),:) = []; Neigh_vec(~any(Neigh_vec,2),:) = []; 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 angle > 0 && angle <= 180; % Do we need this?? if v < w; Angle_index(a,1)= Neigh_ind(v,1); Angle_index(a,2)= i; 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)= i; 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 if mod(i,100)==1 i-1 end ```
```end [Y,I]=sort(Bond_index(:,1)); Bond_index=Bond_index(I,:); Bond_index = unique(Bond_index,'rows','stable'); [Y,I]=sort(Angle_index(:,2)); Angle_index=Angle_index(I,:); Angle_index = unique(Angle_index,'rows','stable'); % [Y,I]=sort(overlap_index(:,2)); % overlap_index=overlap_index(I,:); % overlap_index = unique(overlap_index,'rows','stable'); Bond_index(~any(Bond_index,2),:) = []; Angle_index(~any(Angle_index,2),:) = []; % overlap_index(~any(overlap_index,2),:) = []; nBonds=size(Bond_index,1); nAngles=size(Angle_index,1); if nargin > 4; %% This will print a whole lot more info to the calling workspace radius=repmat(2.3,nAtoms,1); % radius_ion([find(strncmpi([atom.type],'Ow',2)) find(strncmpi([atom.type],'H',1))])=max_short_dist; % This is a test % radius_ion(find(strncmpi([atom.type],'H',1)))=max_short_dist; % dist_matrix = dist_matrix_atom(atom,Box_dim); for i=1:size(atom,2) Neigh_ind=intersect(find(dist_matrix(:,i)>0),find(dist_matrix(:,i)<max_neigh_distance));%radius_ion([atom.type]))); Neigh_dist=dist_matrix(Neigh_ind,i); atom(i).neigh.type = [atom(Neigh_ind).type]'; atom(i).neigh.index = Neigh_ind; atom(i).neigh.dist = Neigh_dist; atom(i).neigh.coords = [[atom(Neigh_ind).x]' [atom(Neigh_ind).y]' [atom(Neigh_ind).z]']; % PBC NOT TAKEN INTO ACCOUNT!!! atom(i).bond.type = []; atom(i).bond.index = []; atom(i).bond.dist = []; atom(i).angle.type = []; atom(i).angle.index = []; atom(i).angle.angle = []; atom(i).angle.vec1 = []; atom(i).angle.vec2 = []; if ismember(i,Bond_index(:,1:2)) [A,B]=find(Bond_index(:,1:2)==i); atom(i).bond.type = 1; atom(i).bond.index = Bond_index(A,1:2); atom(i).bond.dist = Bond_index(A,3); if ismember(i,Angle_index(:,1:3)) [C,D]=find(Angle_index(:,1:3)==i); atom(i).angle.type = 1; atom(i).angle.index = Angle_index(C,1:3); atom(i).angle.angle = Angle_index(C,4); atom(i).angle.vec1 = Angle_index(C,5:7); atom(i).angle.vec2 = Angle_index(C,8:10); end end if mod(i,100)==1 i-1 end end %%%%%%%%%% Atom_labels=unique([atom.type]); for i=1:length(Atom_labels) label_ind=find(strcmpi([atom.type],Atom_labels(i))); Tot_dist=[];Tot_type=[];Tot_index=[];Tot_coords=[];Tot_bonds=[];Tot_angles=[]; for j=label_ind; if numel([atom(j).neigh])>0; Tot_dist=[Tot_dist; [atom(j).neigh.dist]]; Tot_type=[Tot_type; [atom(j).neigh.type]]; Tot_index=[Tot_index; [atom(j).neigh.index]]; Tot_coords=[Tot_coords; [atom(j).neigh.coords]]; end if numel([atom(j).bond])>0; Tot_bonds=[Tot_bonds; [atom(j).bond.dist]]; end if numel([atom(j).angle])>0; Tot_angles=[Tot_angles; [atom(j).angle.angle]]; end end assignin('base',strcat(char(Atom_labels(i)),'_dist')',[num2cell(Tot_index) Tot_type num2cell(Tot_dist)]); assignin('base',strcat(char(Atom_labels(i)),'_coords')',[[atom(Tot_index).x]' [atom(Tot_index).y]' [atom(Tot_index).z]']); assignin('base',strcat(char(Atom_labels(i)),'_bonds')',Tot_bonds); assignin('base',strcat(char(Atom_labels(i)),'_angles')',Tot_angles); assignin('base',strcat(char(Atom_labels(i)),'_atom')',atom(ismember([atom.type],Atom_labels(i)))); end assignin('caller','Dist_matrix',dist_matrix); end assignin('caller','overlap_index',overlap_index); assignin('caller','Bond_index',Bond_index); assignin('caller','Angle_index',Angle_index); assignin('caller','nBonds',nBonds); assignin('caller','nAngles',nAngles); % max_short_dist % max_long_dist % toc ```