MATLAB Examples

merge_atom.m

  • This function returns the atom2 struct with atoms in the atom2 struct
  • with a distance r [1x1 or 1x2] away from the atoms in the atom1 struct. There is
  • also a possibility to use a twin-range cutoff approach (suitable for
  • OH2), by setting r(2) to a smaller value than r(1).
  • You must decide if atom1w should be wrapped or not before running this function
  • Let's not wrap atomw1 here, do it before calling this function since Box1
  • does not always equal the full Box_dim, but ratherr a region in the full Box_dim
  • Tested 15/04/2017
  • Please report bugs to michael.holmboe@umu.se

Contents

Examples

  • atom2w = merge_atom(SOLUTE,Box_dim,SOLVENT,'index','C',1.4)
  • atom2w = merge_atom(SOLUTE,Box_dim,SOLVENT,'molid','Hw',[1.6 1.0])
function atom2w = merge_atom(atom1,Box1,atom2,type,Atom_label,r)
% % % %% For testing the script
% % % clear all;
% % % format compact;
% % % atom1=import_atom('newTCH.gro');Box1=Box_dim;% atom1=translate_atom(atom1,[0 0 Box1(3)/2],'all');
% % % atom2=import_atom('spc216_1.00.gro');Box2=Box_dim;
% % % atom2 = replicate_atom(atom2,Box2,[5 5 3]);Box2=Box_dim;
% % % atom1=center_atom(atom1,Box1,'all','xyz');
% % % atom2=center_atom(atom2,Box2,'all','xyz');
% % % r=[2.2 3.2];
% % % limits=[0 0 0 Box1]
% % % Atom_label={'H' 'O'};
% % % type='molid';
% % %%

% if max([atom1.x])>Box1(1)
%     disp('Box1 smaller that Box_dim in x direction, perhaps using unwrapped atom1?')
%     disp('Will set Box1(1) to max([atom1.x])+.001')
%     Box1(1)=max([atom1.x])+.001;
% end
%
% if max([atom1.y])>Box1(2)
%     disp('Box1 smaller that Box_dim in y direction, perhaps using unwrapped atom1?')
%     disp('Will set Box1(2) to max([atom1.y])+.001')
%     Box1(2)=max([atom1.y])+.001;
% end
%
% if max([atom1.z])>Box1(3)
%     disp('Box1 smaller that Box_dim in z direction, perhaps using unwrapped atom1?')
%     disp('Will set Box1(3) to max([atom1.z])+.001')
%     Box1(3)=max([atom1.z])+.001;
% end

minAtom=Atom_label(1);
atom1w=wrap_atom(atom1,Box1);
% atom1w=atom1;% atom1 should hence be wrapped beforehand
atom2w=slice_molid(atom2,[0 0 0 Box1]);
nAtoms1=size(atom1w,2);
nAtoms2=size(atom2w,2);
indvec=zeros(1,nAtoms2);

if size(atom2w,2)>6000
    % To do it stepwise to save memory...
    disp('Splitting up the distance matrix')
    size(atom1w);
    size(atom2w);
    X1_dist=pdist2([atom1w(1:floor(nAtoms1/2)).x]',[atom2w(1:floor(nAtoms2/2)).x]');X1_dist(X1_dist>Box1(1)/2)=X1_dist(X1_dist>Box1(1)/2)-Box1(1);
    Y1_dist=pdist2([atom1w(1:floor(nAtoms1/2)).y]',[atom2w(1:floor(nAtoms2/2)).y]');Y1_dist(Y1_dist>Box1(2)/2)=Y1_dist(Y1_dist>Box1(2)/2)-Box1(2);
    Z1_dist=pdist2([atom1w(1:floor(nAtoms1/2)).z]',[atom2w(1:floor(nAtoms2/2)).z]');Z1_dist(Z1_dist>Box1(3)/2)=Z1_dist(Z1_dist>Box1(3)/2)-Box1(3);

    X2_dist=pdist2([atom1w(1:floor(nAtoms1/2)).x]',[atom2w(floor((nAtoms2)/2)+1:end).x]');X2_dist(X2_dist>Box1(1)/2)=X2_dist(X2_dist>Box1(1)/2)-Box1(1);
    Y2_dist=pdist2([atom1w(1:floor(nAtoms1/2)).y]',[atom2w(floor((nAtoms2)/2)+1:end).y]');Y2_dist(Y2_dist>Box1(2)/2)=Y2_dist(Y2_dist>Box1(2)/2)-Box1(2);
    Z2_dist=pdist2([atom1w(1:floor(nAtoms1/2)).z]',[atom2w(floor((nAtoms2)/2)+1:end).z]');Z2_dist(Z2_dist>Box1(3)/2)=Z2_dist(Z2_dist>Box1(3)/2)-Box1(3);

    X3_dist=pdist2([atom1w(floor((nAtoms1)/2)+1:end).x]',[atom2w(1:floor(nAtoms2/2)).x]');X3_dist(X3_dist>Box1(1)/2)=X3_dist(X3_dist>Box1(1)/2)-Box1(1);
    Y3_dist=pdist2([atom1w(floor((nAtoms1)/2)+1:end).y]',[atom2w(1:floor(nAtoms2/2)).y]');Y3_dist(Y3_dist>Box1(2)/2)=Y3_dist(Y3_dist>Box1(2)/2)-Box1(2);
    Z3_dist=pdist2([atom1w(floor((nAtoms1)/2)+1:end).z]',[atom2w(1:floor(nAtoms2/2)).z]');Z3_dist(Z3_dist>Box1(3)/2)=Z3_dist(Z3_dist>Box1(3)/2)-Box1(3);

    X4_dist=pdist2([atom1w(floor((nAtoms1)/2)+1:end).x]',[atom2w(floor((nAtoms2)/2)+1:end).x]');X4_dist(X4_dist>Box1(1)/2)=X4_dist(X4_dist>Box1(1)/2)-Box1(1);
    Y4_dist=pdist2([atom1w(floor((nAtoms1)/2)+1:end).y]',[atom2w(floor((nAtoms2)/2)+1:end).y]');Y4_dist(Y4_dist>Box1(2)/2)=Y4_dist(Y4_dist>Box1(2)/2)-Box1(2);
    Z4_dist=pdist2([atom1w(floor((nAtoms1)/2)+1:end).z]',[atom2w(floor((nAtoms2)/2)+1:end).z]');Z4_dist(Z4_dist>Box1(3)/2)=Z4_dist(Z4_dist>Box1(3)/2)-Box1(3);

    X_dist=[X1_dist X2_dist;X3_dist X4_dist];
    Y_dist=[Y1_dist Y2_dist;Y3_dist Y4_dist];
    Z_dist=[Z1_dist Z2_dist;Z3_dist Z4_dist];

    dist_matrix=(X_dist.^2+Y_dist.^2+Z_dist.^2).^.5;
    clear X*dist Y*dist Z*dist;
elseif size(atom2w,2)>0
    % To do them all at once
    %disp('Single distance matrix')
    size([atom1w.x]');
    size([atom2w.x]');
    X_dist=pdist2([atom1w.x]',[atom2w.x]');
    X_dist(X_dist>Box1(1)/2)=X_dist(X_dist>Box1(1)/2)-Box1(1);
    Y_dist=pdist2([atom1w.y]',[atom2w.y]');
    Y_dist(Y_dist>Box1(2)/2)=Y_dist(Y_dist>Box1(2)/2)-Box1(2);
    Z_dist=pdist2([atom1w.z]',[atom2w.z]');
    Z_dist(Z_dist>Box1(3)/2)=Z_dist(Z_dist>Box1(3)/2)-Box1(3);
    dist_matrix=(X_dist.^2+Y_dist.^2+Z_dist.^2).^.5;
    clear X*dist Y*dist Z*dist;
else
    disp('atom2w is empty')
end

% vmd([atom1w atom2w],Box1)
% disp('sizes X,Y,Z and distance matrix')
% size(X_dist)
%
% size(Y_dist)
%
% size(Z_dist)
%
% size(dist_matrix)



% This section loops though the atom1 atoms
% indvec=zeros(1,size(dist_matrix,2));
% for i=1:size(dist_matrix,1);
%     indvec(dist_matrix(i,:)<r)=1;
% end
if exist('dist_matrix','var')
    % This section loops though the atom2 atoms, using rsmall and rlarge, where rsmall operates on minAtom=Atom_label(1)
    if numel(r)==1
        rsmall=r;
        rlarge=r;
    else
        rsmall=r(1);
        rlarge=r(2);
    end

    for i=1:size(dist_matrix,2)
        if strncmpi([atom2w(i).type],minAtom,1)
            if any((dist_matrix(:,i)-rsmall)<0)
                indvec(1,i)=1;
            end
        else
            if any((dist_matrix(:,i)-rlarge)<0)
                indvec(1,i)=1;
            end
        end
    end

    clear dist_matrix;

    removed_molid=unique([atom2w(find(indvec)).molid]);
    removed_molid_index=[atom2w(ismember([atom2w.molid],removed_molid)).index];
    removed_index=find(indvec);

    pre_size=size(atom2w,2);
    if strcmpi(type,'molid')
        atom2w(ismember([atom2w.molid],removed_molid))=[];
    else
        atom2w(removed_index)=[];
    end
    post_size=size(atom2w,2);

%     vmd([atom1w atom2w],Box1)

    % disp('Removed this many atoms')
    % pre_size-post_size

    assignin('caller','removed_molid',removed_molid);
    assignin('caller','removed_molid_index',removed_molid_index);
    assignin('caller','removed_index',removed_index);

end
% vmd([atom1w atom2w],Box1)