replace_atom.m

Contents

Version

2.11

Contact

Please report problems/bugs to michael.holmboe@umu.se

Examples

  1. atom = replace_atom(in_atom,old_atom,3)
function atom=replace_atom(new_atom,prev_atom,MolID)

disp('Assuming all atom structs are unwrapped and whole...')

old_COM=zeros(numel(MolID),3);
for i=1:numel(MolID)
    temp = COM_atom(prev_atom([prev_atom.molid]==MolID(i)));
    position=COM;
    if size(new_atom,2) < 500
        temp_atom=COM_atom(new_atom); % This generates the COM position through assignin. You should use unwrapped molecules
        temp_atom=rmfield(temp_atom,'COM_x');
        temp_atom=rmfield(temp_atom,'COM_y');
        temp_atom=rmfield(temp_atom,'COM_z');
        temp_atom=rmfield(temp_atom,'element');
        temp_atom=rmfield(temp_atom,'Mw');
    else
        COM=[mean([temp_atom.x]) mean([temp_atom.y]) mean([temp_atom.z])]; % Since COM_atom is a bit slow for large molecules, we do this for big molecules
    end

    temp_atom = translate_atom(new_atom,-COM+position,'all');
    [temp_atom.molid]=deal(MolID(i));

    ind_prev=find(ismember([prev_atom.molid],MolID(i)));
    if MolID(i)==min([prev_atom.molid]) % Put the new molid first
        prev_atom=update_atom({temp_atom prev_atom(max(ind_prev)+1:end)});
    elseif MolID(i)==max([prev_atom.molid]) % Put the new molid last
        prev_atom=update_atom({temp_atom prev_atom(1:min(ind_prev)-1)});
    else % Put the new molid in between somewhere
        prev_atom=update_atom({prev_atom(1:min(ind_prev)-1) temp_atom prev_atom(max(ind_prev)+1:end)});
    end

end

% Do we need this?
atom=update_atom(prev_atom);