MATLAB Examples

interface15_param.m

  • This function holds the extended interface ff parameters
  • Call function with find(ismember({forcefield.interface15.atom},'Hw')) to find Hw's index
  • Tested 15/04/2017
  • Please report bugs to michael.holmboe@umu.se

Contents

Examples

  • interface15_param({'AY1'},'tip3p')
function interface15_param(Atom_label,water_model)
%

% Original from charmm27_interface15_v1_5.prm with comparison to clayff (MHolmobes atomnames...)
% !CLAY MINERALS (KAOLINITE, MICA, MONTMORILLONITES, PYROPHYLLITE)
% !
% Heinz, H.; Koerner, H.; Anderson, K. L.; Vaia, R. A. and Farmer, B. L., Force Field for Mica-Type
% Silicates and Dynamics of Octadecylammonium Chains Grafted to Montmorillonite, Chemistry of Materials,
% 2005, 17, 5658-5669; and extensions to kaolinite.															; Clayff?
% K+     Potassium ion in mica, mmt, and other clays 		(+1.0)											; K 			1
% NA+    Sodium ion in silica, clay, and organic salts	    (+1.0)											; Na 			1
% SY1    Silicon atom in tetrahedral silicate sheet 		(+1.1)											; Si 			2.1
% SY2    Silicon atom in tetrahedral silicate sheet 		(+1.1)											; Si 			2.1
% AYT1   Aluminum defect of SY1 							(+0.8)											; Alt 			1.575
% AYT2   Aluminum defect of SY2 							(+0.8)											; Alt 			1.575
% AY1    Aluminum atom in octahedral aluminate sheet		(+1.45) or magnesium defect	(+1.1)		 		; Al or Mgo		1.575 or 1.36
% AY2    Aluminum atom in octahedral aluminate sheet		(+1.45) or magnesium defect	(+1.1)		 		; Al or Mgo		1.575 or 1.36
% OY1    Oxygen atom in silicate sheet, surface 			(-0.55(000), -0.783(33) if next to Al defect)	; O or Omg		-1.05 or -1.1688
% OY2    Oxygen atom in silicate sheet, surface 			(-0.55(000), -0.783(33) if next to Al defect)	; O or Omg		-1.05 or -1.1688
% OY3    Oxygen atom in silicate sheet, surface 			(-0.55(000), -0.783(33) if next to Al defect)	; O or Omg		-1.05 or -1.1688
% OY4    Oxygen atom in octahedral aluminate sheet, apical 	(-0.758(33), -0.866(66) if next to Mg defect)	; O or Omg		-1.05 or -1.1808
% OY5    Oxygen atom in octahedral aluminate sheet, apical 	(-0.758(33), -0.866(66) if next to Mg defect)	; O or Omg		-1.05 or -1.1808
% OY6    Oxygen atom in octahedral aluminate sheet, hydroxl	(-0.683(33), -0.791(66) if next to Mg defect)	; Oh or Ohmg	-0.95 or -1.0808
% OY7    Oxygen atom in octahedral aluminate sheet, apical 	(-0.758(33), -0.866(66) if next to Mg defect)	; O or Omg		-1.05 or -1.1808
% OY8    Oxygen atom in octahedral aluminate sheet, apical 	(-0.758(33), -0.866(66) if next to Mg defect)	; O or Omg		-1.05 or -1.1808
% OY9    Oxygen atom in octahedral aluminate sheet, hydroxl (-0.683(33), -0.791(66) if next to Mg defect)	; Oh or Ohmg	-0.95 or -1.0808
% HOY    Hydrogen atom in (Al,Mg,Si)OH groups in clay 		(+0.2), silica (+0.4), cement (excl. ettr)		; H				0.425
% HOK    Hydrogen atom in AlOH surface group in kaolinite 	(+0.2)											; H				0.425

% New additions by M Holmboe
% MY1   Octahedral Magnesium defect, replacing Aluminum atom in octahedral aluminate sheet (+1.1)		 	; Al or Mgo		1.575 or 1.36
% Edge O atoms starts here... by M Holmboe
% OAH
% OAHH
% OSH
% OAS
% OASH

% from interface15_param.m  1    2    3    4   5    6     7      8      9      10      11      12    13    14   15   16    17    18   19  20   21   22   23   24   25   26  27    28   29   30  31   32
%    Atom_label_Heinz =  {'Hw', 'H','Ow','Oh','Ob','Op','Omg','Oalt','Ohmg','Oalsi','Oalhh','Oalh','Osih','Si','Al','Alt','Mgo','Li','Na','K','Cs','Mg','Ca','Sr','Ba','F','Cl','Cls','Br','I','He','Du'}';

% Hw
forcefield.interface15(1).type =    {'HW'}; % h* water hydrogen
forcefield.interface15(1).mass =    1.00794;
if strcmpi('SPC/E', water_model) > 0 || strcmpi('SPCE', water_model) > 0;
    forcefield.interface15(1).charge =  0.42380;
elseif strcmp('SPC', water_model) > 0;
    forcefield.interface15(1).charge =  0.41000;
elseif strncmpi('TIP3P', water_model,5) > 0;
    forcefield.interface15(1).charge =  0.41700;
end
forcefield.interface15(1).radius =  0.00000; % R0 Å
forcefield.interface15(1).e_kJmol = 0.00000; % kcal mol-1
%
% H
forcefield.interface15(2).type =    {'HOY'}; % ho hydroxyl hydrogen
forcefield.interface15(2).mass =    1.00794;
forcefield.interface15(2).charge =  0.200;
forcefield.interface15(2).radius =  1.13925;
forcefield.interface15(2).e_kJmol = 0.06276;
%
% Ow
forcefield.interface15(3).type =    {'OW'}; % o* water oxygen
forcefield.interface15(3).mass =    15.99410;
if strcmpi('SPC/E', water_model) > 0 || strcmpi('SPCE', water_model) > 0;
    forcefield.interface15(3).charge =  -0.84760;
elseif strcmp('SPC', water_model) > 0;
    forcefield.interface15(3).charge = -0.82000;
elseif strncmpi('TIP3P', water_model,5) > 0;
    forcefield.interface15(3).charge =  -0.83400;
end
forcefield.interface15(3).radius =  3.553200;
forcefield.interface15(3).e_kJmol = 0.155400;
%
% Oh
forcefield.interface15(4).type =    {'OY6'}; % Oxygen atom in octahedral aluminate sheet, hydroxl
forcefield.interface15(4).mass =    15.99410;
forcefield.interface15(4).charge =  -0.68333;
forcefield.interface15(4).radius =  3.67500;
forcefield.interface15(4).e_kJmol = 0.104600;
%
% Ob
forcefield.interface15(5).type =    {'OY1'}; % Oxygen atom in silicate sheet, surface
forcefield.interface15(5).mass =    15.99410;
forcefield.interface15(5).charge =  -0.55000;
forcefield.interface15(5).radius =  3.67500;
forcefield.interface15(5).e_kJmol = 0.104600;
%
% Op
forcefield.interface15(6).type =    {'OY4'}; % Oxygen atom in octahedral aluminate sheet, apical
forcefield.interface15(6).mass =    15.99410;
forcefield.interface15(6).charge =  -0.75833;
forcefield.interface15(6).radius =  3.67500;
forcefield.interface15(6).e_kJmol = 0.104600;
% Omg
forcefield.interface15(7).type =    {'OY5'}; % Oxygen atom in octahedral aluminate sheet, apical, if next to Mg defect
forcefield.interface15(7).mass =    15.99410;
forcefield.interface15(7).charge =  -0.86666;
forcefield.interface15(7).radius =  3.67500;
forcefield.interface15(7).e_kJmol = 0.104600;
%
% Oalt
forcefield.interface15(8).type =    {'OY2'}; % Oxygen atom in silicate sheet, surface, if next to Al defect
forcefield.interface15(8).mass =    15.99410;
forcefield.interface15(8).charge =  -0.78333;
forcefield.interface15(8).radius =  3.67500;
forcefield.interface15(8).e_kJmol = 0.104600;
%
% % Op
% forcefield.interface15(6).type =    {'OY3'}; % Oxygen atom in octahedral aluminate sheet, apical
% forcefield.interface15(6).mass =    15.99410;
% forcefield.interface15(6).charge =  -0.75833;
% forcefield.interface15(6).radius =  3.67500;
% forcefield.interface15(6).e_kJmol = 0.104600;
% %

% % Oh
% forcefield.interface15(4).type =    {'OY7'}; % Oxygen atom in octahedral aluminate sheet, hydroxl
% forcefield.interface15(4).mass =    15.99410;
% forcefield.interface15(4).charge =  -0.68333;
% forcefield.interface15(4).radius =  3.67500;
% forcefield.interface15(4).e_kJmol = 0.104600;
% %
% % Oh
% forcefield.interface15(4).type =    {'OY8'}; % Oxygen atom in octahedral aluminate sheet, hydroxl
% forcefield.interface15(4).mass =    15.99410;
% forcefield.interface15(4).charge =  -0.68333;
% forcefield.interface15(4).radius =  3.67500;
% forcefield.interface15(4).e_kJmol = 0.104600;
% %
% Ohmg
forcefield.interface15(9).type =    {'OY9'}; % Oxygen atom in octahedral aluminate sheet, hydroxl, if next to Mg defect
forcefield.interface15(9).mass =    15.99410;
forcefield.interface15(9).charge =  -0.79166;
forcefield.interface15(9).radius =  3.67500;
forcefield.interface15(9).e_kJmol = 0.104600;
%
% Oalh
forcefield.interface15(10).type =    {'OAH'}; %
forcefield.interface15(10).mass =    15.99410;
forcefield.interface15(10).charge =  -0.75833; % Dummy charge
forcefield.interface15(10).radius =  3.67500;
forcefield.interface15(10).e_kJmol = 0.104600;
%
% Oalhh
forcefield.interface15(11).type =    {'OAHH'}; % obts bridging oxygen w. tetrahedral substitution
forcefield.interface15(11).mass =    15.99410;
forcefield.interface15(11).charge =  -0.4000; %??? Where does it come from? -0.07500;
forcefield.interface15(11).radius =  3.67500;
forcefield.interface15(11).e_kJmol = 0.104600;
%
% Oalsi
forcefield.interface15(12).type =    {'OAS'}; % obss bridging oxygen w. double substitution
forcefield.interface15(12).mass =    15.99410;
forcefield.interface15(12).charge =  -0.75833; %
forcefield.interface15(12).radius =  3.67500;
forcefield.interface15(12).e_kJmol = 0.104600;
%
% Osih
forcefield.interface15(13).type =    {'OSH'}; % obss bridging oxygen w. double substitution
forcefield.interface15(13).mass =    15.99410;
forcefield.interface15(13).charge =  -0.675;
forcefield.interface15(13).radius =  3.67500;
forcefield.interface15(13).e_kJmol = 0.104600;
%
% Si
forcefield.interface15(14).type =    {'SY1'}; % st tetrahedral silicon
forcefield.interface15(14).mass =    28.08538;
forcefield.interface15(14).charge =  1.1000;
forcefield.interface15(14).radius =  4.20000;
forcefield.interface15(14).e_kJmol = 0.20920;
%
%
% % Si
% forcefield.interface15(14).type =    {'SY2'}; % st tetrahedral silicon
% forcefield.interface15(14).mass =    28.08538;
% forcefield.interface15(14).charge =  1.1000;
% forcefield.interface15(14).radius =  4.20000;
% forcefield.interface15(14).e_kJmol = 0.20920;
%
% Al
forcefield.interface15(15).type =    {'AY1'}; % ao octahedral aluminium
forcefield.interface15(15).mass =    26.98154;
forcefield.interface15(15).charge =  1.45;% 1.283334; % To smooth/average charge Al/Mg charge
forcefield.interface15(15).radius =  4.41000;
forcefield.interface15(15).e_kJmol = 0.20920;
%
% % Al
% forcefield.interface15(15).type =    {'AY2'}; % ao octahedral aluminium
% forcefield.interface15(15).mass =    26.98154;
% forcefield.interface15(15).charge =  1.45;% 1.283334; % To smooth/average charge Al/Mg charge
% forcefield.interface15(15).radius =  4.41000;
% forcefield.interface15(15).e_kJmol = 0.20920;
%
% Alt
forcefield.interface15(16).type =    {'AYT1'}; % at tetrahedral aluminium
forcefield.interface15(16).mass =    26.98154;
forcefield.interface15(16).charge =  1.45;
forcefield.interface15(16).radius =  4.41000;
forcefield.interface15(16).e_kJmol = 0.20920;
% Mgo
forcefield.interface15(17).type =    {'MY1'}; % mgo
forcefield.interface15(17).mass =    24.305;
forcefield.interface15(17).charge =  1.1000;%-0.002574/32;%-8.1250e-05;
forcefield.interface15(17).radius =  4.41000;
forcefield.interface15(17).e_kJmol = 0.20920;

% Ion pair-potentials from Joung-Cheatham 2008 and Babu-Lim, 20??
% Li
forcefield.interface15(18).type =    {'Li'}; %
forcefield.interface15(18).mass =    6.94;
forcefield.interface15(18).charge =  1.000;
forcefield.interface15(18).radius =  1.582;
forcefield.interface15(18).e_kJmol = 1.40890;
%
% Na
forcefield.interface15(19).type =    {'Na'}; %
forcefield.interface15(19).mass =    22.98977;
forcefield.interface15(19).charge =  1.000;
forcefield.interface15(19).radius =  2.424;
forcefield.interface15(19).e_kJmol = 1.47545;
%
% K
forcefield.interface15(20).type =    {'K'}; %
forcefield.interface15(20).mass =    39.09830;
forcefield.interface15(20).charge =  1.000;
forcefield.interface15(20).radius =  3.186;
forcefield.interface15(20).e_kJmol = 1.79789;
%
% Cs
forcefield.interface15(21).type =    {'Cs'}; %
forcefield.interface15(21).mass =    132.90545;
forcefield.interface15(21).charge =  1.000;
forcefield.interface15(21).radius =  4.042;
forcefield.interface15(21).e_kJmol = 0.37596;
%
% Mg
forcefield.interface15(22).type =    {'Mg'}; %
forcefield.interface15(22).mass =    24.305;
forcefield.interface15(22).charge =  2.000;
forcefield.interface15(22).radius =  1.56880;
forcefield.interface15(22).e_kJmol = 3.65874;
%
% Ca
forcefield.interface15(23).type =    {'Ca'}; %
forcefield.interface15(23).mass =    40.078;
forcefield.interface15(23).charge =  2.000;
forcefield.interface15(23).radius =  2.6500;
forcefield.interface15(23).e_kJmol = 1.88011;
%
% Sr
forcefield.interface15(24).type =    {'Sr'}; %
forcefield.interface15(24).mass =    87.620;
forcefield.interface15(24).charge =  2.000;
forcefield.interface15(24).radius =  3.4823;
forcefield.interface15(24).e_kJmol = 0.49433;
%
% Ba
forcefield.interface15(25).type =    {'Ba'}; %
forcefield.interface15(25).mass =    137.330;
forcefield.interface15(25).charge =  2.000;
forcefield.interface15(25).radius =  4.24985;
forcefield.interface15(25).e_kJmol = 0.19692;
%
% F
forcefield.interface15(26).type =    {'F'}; %
forcefield.interface15(26).mass =    18.998403;
forcefield.interface15(26).charge =  -1.000;
forcefield.interface15(26).radius =  4.514;
forcefield.interface15(26).e_kJmol = 0.0074005;
%
% Cl
forcefield.interface15(27).type =    {'Cl'}; %
forcefield.interface15(27).mass =    35.4530;
forcefield.interface15(27).charge =  -1.000;
forcefield.interface15(27).radius =  5.422;
forcefield.interface15(27).e_kJmol = 0.05349;
%
% Cls
forcefield.interface15(28).type =    {'Cls'}; %
forcefield.interface15(28).mass =    35.4530;
forcefield.interface15(28).charge =  -1.000;
forcefield.interface15(28).radius =  5.422;
forcefield.interface15(28).e_kJmol = 0.05349;
%
% Br
forcefield.interface15(29).type =    {'Br'}; %
forcefield.interface15(29).mass =    79.9040;
forcefield.interface15(29).charge =  -1.000;
forcefield.interface15(29).radius =  5.502;
forcefield.interface15(29).e_kJmol = 0.11279;
%
% I
forcefield.interface15(30).type =    {'I'}; %
forcefield.interface15(30).mass =    126.90447;
forcefield.interface15(30).charge =  -1.000;
forcefield.interface15(30).radius =  5.838;
forcefield.interface15(30).e_kJmol = 0.17901;
%
% H
forcefield.interface15(31).type =    {'HE'}; % Edge H
forcefield.interface15(31).mass =    1.00794;
forcefield.interface15(31).charge =  0.400;
forcefield.interface15(31).radius =  1.13925;
forcefield.interface15(31).e_kJmol = 0.06276;
%
% Du
forcefield.interface15(32).type =    {'Du'}; % Edge H
forcefield.interface15(32).mass =    10.0;
forcefield.interface15(32).charge =  0.00;
forcefield.interface15(32).radius =  0.00;
forcefield.interface15(32).e_kJmol = 0.00;
% Oahs
forcefield.interface15(33).type =    {'OAHS'}; %
forcefield.interface15(33).mass =    15.99410;
forcefield.interface15(33).charge =  -0.6625;
forcefield.interface15(33).radius =  3.553200;
forcefield.interface15(33).e_kcalmol = 0.155400;
% Oahs
forcefield.interface15(34).type =    {'OMS'}; %
forcefield.interface15(34).mass =    15.99410;
forcefield.interface15(34).charge =  -0.9755;
forcefield.interface15(34).radius =  3.553200;
forcefield.interface15(34).e_kcalmol = 0.155400;
% Al -edge
forcefield.interface15(35).type =    {'AYE'}; % ao octahedral aluminium
forcefield.interface15(35).mass =    26.98154;
forcefield.interface15(35).charge =  1.8125; %1.5750+0.2375;
forcefield.interface15(35).radius =  4.79430;
forcefield.interface15(35).e_kcalmol = 1.3298E-06;
% Ar
forcefield.interface15(36).type =    {'Ar'}; % Ar with SPC..
forcefield.interface15(36).mass =    39.94800;
forcefield.interface15(36).charge =  0.00000;
forcefield.interface15(36).radius =  3.775333752;
forcefield.interface15(36).e_kcalmol = 0.276720296;
% Co
forcefield.interface15(37).type =    {'Co'}; % C as in CO2
forcefield.interface15(37).mass =    12.01073;
forcefield.interface15(37).charge =  0.6512;
forcefield.interface15(37).radius =  3.094627867;
forcefield.interface15(37).e_kcalmol = 0.05580526;
% Oc
forcefield.interface15(38).type =    {'Oc'}; % O as in CO2
forcefield.interface15(38).mass =    15.99941;
forcefield.interface15(38).charge =  -0.3256	;
forcefield.interface15(38).radius =  3.404427393;
forcefield.interface15(38).e_kcalmol = 0.160036571;


% Forcefield_index=find(ismember([forcefield.interface15.type],Atom_label))';
% forcefield.interface15=forcefield.interface15(Forcefield_index);

for i=1:size(forcefield.interface15,2);
    forcefield.interface15(i).sigma=forcefield.interface15(i).radius/2^(1/6);
    %     forcefield.interface15(i).e_kJmol=forcefield.interface15(i).e_kcalmol*4.184;
    forcefield.interface15(i).e_eV=forcefield.interface15(i).e_kJmol*1000/6.02214149999999E23/1.60217733E-19;
end

for i = 1:length(Atom_label)
    if strcmp(Atom_label(i),{'O'})
        Atom_label(i)={'OY1'};
    end
    if strncmpi(Atom_label(i),'Hw',2)
        Atom_label(i)={'Hw'};
    end

    ind=find(strcmpi([forcefield.interface15.type],Atom_label(i)))';
    if numel(ind)==0;
        ind=find(strcmpi([forcefield.interface15.type],'Du'))';
    end

    ind=find(strcmpi([forcefield.interface15.type],Atom_label(i)))';
    ff.interface15(i).fftype =  [forcefield.interface15(ind).type];
    ff.interface15(i).masses =  [forcefield.interface15(ind).mass];
    ff.interface15(i).charge =  [forcefield.interface15(ind).charge];
    ff.interface15(i).epsilon =  [forcefield.interface15(ind).e_kJmol];
    ff.interface15(i).radius =  [forcefield.interface15(ind).radius];
    ff.interface15(i).sigma =  [forcefield.interface15(ind).sigma];
end

fftype  = [ff.interface15.fftype];
Masses  = [ff.interface15.masses];
Charge  = [ff.interface15.charge];
Epsilon = [ff.interface15.epsilon];
Radius  = [ff.interface15.radius];
Sigma   = [ff.interface15.sigma];

All_Bonds =  [554.1349 1.000;
              553.7600 1.000];

% All_Angles = [45.7697 109.470;
%     30.0000 109.470;
%     30.0000 109.470];

% from Teleman Jönsson, flex SPC model, eV A^2
% 24.029546 (*2) Bond
% 1.984757 % Angle

All_Angles = [45.7696 109.470];

kcalmol_to_eV = 4.184*1000/(6.0221415E+23*1.60217733E-19);

typeID = [];
e_mix = zeros(1,1);
r_mix = zeros(1,1);
k=1;

% if strcmp('SPC/E', water_model,'exact') > 0;
%     ind_Hw=find(ismember({forcefield.interface15.type},'Hw')) =  0.4238;%    h*      Hw %SPCE water
%     ind_Ow=find(ismember({forcefield.interface15.type},'Ow')) = -0.8476; %	o*      Ow       %SPCE water
% end


% Send these variables to the calling script too!
assignin('caller','kcalmol_to_eV',kcalmol_to_eV);
assignin('caller','typeID',typeID);
assignin('caller','Masses',Masses);
assignin('caller','Charge',Charge);
assignin('caller','e_mix',e_mix);
assignin('caller','r_mix',r_mix);
assignin('caller','bondstr',All_Bonds);
assignin('caller','anglebend',All_Angles);
assignin('caller','Radius',Radius);
assignin('caller','ff',ff);
assignin('caller','fftype',fftype);
assignin('caller','forcefield',forcefield);

try
    for i=1:length(Atom_label)
        j = i;
        while j < length(Atom_label)+1
            %         Atom_label(i)
            %         Atom_label(j)
            %         typeID = [typeID;[i j]]
            %         e_mix
            %         Sigma
            e_mix(k,1) = (Epsilon(i)*Epsilon(j))^.5;
            r_mix(k,1) = (Sigma(i)+Sigma(j))/2;
            j = j + 1;
            k = k + 1;
        end
    end
    assignin('caller','Epsilon',Epsilon);
    assignin('caller','Sigma',Sigma);
catch
    disp('Unable to assign all atomtypes')
    disp('We had problems with:')
    ff.interface15(i).fftype
    ff.interface15(i).masses
    ff.interface15(i).charge
    ff.interface15(i).epsilon
    ff.interface15(i).radius
    ff.interface15(i).sigma
    disp('or...')
    ff.interface15(j).fftype
    ff.interface15(j).masses
    ff.interface15(j).charge
    ff.interface15(j).epsilon
    ff.interface15(j).radius
    ff.interface15(j).sigma
end

% XYZ_num=zeros(length(XYZ_labels),1);
% Atom_label=sort(unique([atom.type]));
% for i=1:length(Atom_label)
%     XYZ_num(ismember(XYZ_labels,Atom_label(i)))=i;
% end
%
% XYZ_radius=Sigma(XYZ_num(:))';
% A1=repmat(XYZ_radius,1,length(XYZ_radius));
% A2=repmat(XYZ_radius,1,length(XYZ_radius))';
% B=(A1+A2)*.3;
% D=(Dist_matrix-B)<0;
% [C1,C2]=sort(Dist_matrix(:,1));
% C3=XYZ_labels(C2(:));