MATLAB Examples

charge_atom.m

  • This function tries to charge the atom according to clayff or interface ff
  • If more than 4 arguments, it tries to smear out the excess charge
  • atom is the atom struct
  • Box_dim is the box dimension vector
  • ffname is 'clayff' or 'interface'
  • watermodel is not always used but should be 'spc' 'spc/e' 'tip3p'
  • Tested 15/04/2017
  • Please report bugs to michael.holmboe@umu.se

Contents

Examples

  • atom = charge_atom(atom,Box_dim,'clayff','spc')
  • atom = charge_atom(atom,Box_dim,'interface','tip3p','more')
function atom = charge_atom(atom,Box_dim,ffname,watermodel,varargin)
if strcmpi(ffname,'clayff')
    clayff_param(sort(unique([atom.type])),watermodel);
    for i=1:length(atom)
        if strncmpi([atom(i).type],{'Hw'},2)
            ind=strncmpi({'Hw'},[forcefield.clayff.type],2);
        else
            ind=strcmpi([atom(i).type],[forcefield.clayff.type]);
        end
        atom(i).charge=[forcefield.clayff(ind).charge];
    end
    if nargin > 4
        %         Atom_label=varargin{1}(:);
        %         Charge=cell2mat(varargin(2));
        % Total_charge = charge_clayff_atom(atom,Box_dim,{'Al' 'Mgo' 'Si' 'H'},[1.575 1.36 2.1 0.425])
        Atom_label=sort(unique([atom.type]));
        clayff_param(sort(Atom_label),watermodel);
        no_adjust_labels=[Atom_label(~strncmp(Atom_label,'O',1))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Ow',2))];
        no_adjust_ind=ismember(Atom_label,no_adjust_labels);
        no_adjust_ind
        Atom_label
        Charge
        atom = charge_clayff_atom(atom,Box_dim,Atom_label(no_adjust_ind),Charge(no_adjust_ind));
    else
        atom=check_clayff_charge(atom);
    end

elseif strcmpi(ffname,'interface')
    interface_param(unique([atom.type]),watermodel);
    for i=1:length(atom)
        if strncmpi([atom(i).type],{'Hw'},2)
            ind=strncmpi({'Hw'},[forcefield.interface.type],2);
        else
            ind=strcmp([atom(i).type],[forcefield.interface.type]);
        end
        atom(i).charge=[forcefield.interface(ind).charge];
    end
    if nargin > 4
        %         Atom_label=varargin{1}(:);
        %         Charge=cell2mat(varargin(2));
        % Total_charge = charge_interface_atom(atom,Box_dim,{'Al' 'Mgo' 'Si' 'H'},[1.575 1.36 2.1 0.425])
        Atom_label=sort(unique([atom.type]));
        interface_param(sort(Atom_label),watermodel);
        no_adjust_labels=[Atom_label(~strncmp(Atom_label,'O',1))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Ow',2))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Ob',2))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Op',2))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Oh',2))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Omg',3))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Ohmg',4))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Osih',4))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Oalsi',5))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Oalhh',5))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Oalh',4))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Oalt',4))];
        no_adjust_ind=ismember(Atom_label,no_adjust_labels);
        atom = charge_interface_atom(atom,Box_dim,Atom_label(no_adjust_ind),Charge(no_adjust_ind));
    else
        atom=check_interface_charge(atom);
    end
elseif strcmpi(ffname,'interface15')
    interface15_param(unique([atom.fftype]),watermodel);
    for i=1:length(atom)
        if strncmpi([atom(i).fftype],{'Hw'},2)
            ind=strncmpi({'Hw'},[forcefield.interface15.type],2);
        else
            ind=strcmp([atom(i).type],[forcefield.interface15.type]);
        end
        atom(i).charge=[forcefield.interface15(ind).charge];
    end
    if nargin > 4
        %         Atom_label=varargin{1}(:);
        %         Charge=cell2mat(varargin(2));
        % Total_charge = charge_interface_atom(atom,Box_dim,{'Al' 'Mgo' 'Si' 'H'},[1.575 1.36 2.1 0.425])
        Atom_label=sort(unique([atom.fftype]));
        interface15_param(sort(Atom_label),watermodel);
        no_adjust_labels=[Atom_label(~strncmp(Atom_label,'O',1))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'Ow',2))];
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OY1',3))]; % Ob
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OY4',3))]; % Op
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OY5',3))]; % Omg
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OY6',3))]; % Oh
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OY9',3))]; % Ohmg
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OSH',3))]; % Osih
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OAS',3))]; % Oalsi
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OAHH',4))];% Oalhh
        no_adjust_labels=[no_adjust_labels Atom_label(strncmp(Atom_label,'OAH',3))]; % Oalh
        no_adjust_ind=ismember(Atom_label,no_adjust_labels);
        Atom_label(no_adjust_ind)
        Charge(no_adjust_ind)
        atom = charge_interface15_atom(atom,Box_dim,Atom_label(no_adjust_ind),Charge(no_adjust_ind));
    else
        atom=check_interface15_charge(atom);
    end
end

disp('Total charge')
Total_charge=sum([atom.charge])
% atom=tweak_charge_atom(atom);
assignin('caller','Total_charge',Total_charge);