MATLAB Examples

write_atom_cif.m

  • This function writes a basic cif file with the fractional coordinates
  • from the atom struct
  • Tested 15/04/2018
  • Please report bugs to michael.holmboe@umu.se

Contents

Examples

  • write_atom_cif(atom,Box_dim,filename_out)
function write_atom_cif(atom,Box_dim,filename_out)

atom=wrap_atom(atom,Box_dim);

if regexp(filename_out,'.cif') ~= false
    filename_out = filename_out;
else
    filename_out = strcat(filename_out,'.cif');
end

if numel(Box_dim)<9
    if numel(Box_dim)==1
        Box_dim(1)=Box_dim(1);
        Box_dim(2)=Box_dim(1);
        Box_dim(3)=Box_dim(1);
    end
    Box_dim(4:9)=0;
end

disp('Assuming P1 space group. Box can still be triclinic')


Box_dim(Box_dim<0.00001&Box_dim>-0.00001)=0;

lx=Box_dim(1);
ly=Box_dim(2);
lz=Box_dim(3);
xy=Box_dim(6);
xz=Box_dim(8);
yz=Box_dim(9);

a=lx;
b=(ly^2+xy^2)^.5;
c=(lz^2+xz^2+yz^2)^.5;

if numel(find(Box_dim(4:9)))>0
    alfa=rad2deg(acos((ly*yz+xy*xz)/(b*c)));
    beta=rad2deg(acos(xz/c));
    gamma=rad2deg(acos(xy/b));
else
    alfa=90;
    beta=90;
    gamma=90;
end

%     a=Box_dim(1);
%     b=Box_dim(2);
%     c=Box_dim(3);
%     xy=Box_dim(6);
%     xz=Box_dim(8);
%     yz=Box_dim(9);
%     lx = a;
%     ly = (b^2-xy^2)^.5;
%     lz = (c^2 - xz^2 - yz^2)^0.5;
%     alfa=rad2deg(acos((ly*yz+xy*xz)/(b*c)))
%     beta=rad2deg(acos(xz/c));
%     gamma=rad2deg(acos(xy/b));


for i=1:size(atom,2)
    if strncmp(atom(i).type,{'Si'},2);atom(i).element={'Si'};atom(i).formalcharge=4;
    elseif strncmpi(atom(i).type,{'SY'},2);atom(i).element={'Si'};atom(i).formalcharge=4;
    elseif strncmpi(atom(i).type,{'SC'},2);atom(i).element={'Si'};atom(i).formalcharge=4;
    elseif strncmpi(atom(i).type,{'S'},1);atom(i).element={'S'};atom(i).formalcharge=2;
    elseif strncmpi(atom(i).type,{'AC'},2);atom(i).element={'Al'};atom(i).formalcharge=3;
    elseif strncmpi(atom(i).type,{'AY'},2);atom(i).element={'Al'};atom(i).formalcharge=3;
    elseif strncmpi(atom(i).type,{'Al'},2);atom(i).element={'Al'};atom(i).formalcharge=3;
    elseif strncmpi(atom(i).type,{'Mg'},2);atom(i).element={'Mg'};atom(i).formalcharge=2;
    elseif strncmpi(atom(i).type,{'Fe'},2);atom(i).element={'Fe'};atom(i).formalcharge=3;
    elseif strncmpi(atom(i).type,{'O'},1);atom(i).element={'O'};atom(i).formalcharge=-2;
    elseif strncmpi(atom(i).type,{'H'},1);atom(i).element={'H'};atom(i).formalcharge=1;
    elseif strncmpi(atom(i).type,{'K'},1);atom(i).element={'K'};atom(i).formalcharge=1;
    elseif strncmpi(atom(i).type,{'Na'},1);atom(i).element={'Na'};atom(i).formalcharge=0;
    elseif strncmpi(atom(i).type,{'Cl'},2);atom(i).element={'Cl'};atom(i).formalcharge=-1;
    elseif strncmpi(atom(i).type,{'Br'},2);atom(i).element={'Br'};atom(i).formalcharge=-1;
    elseif strncmpi(atom(i).type,{'Ca'},2);atom(i).element={'Ca'};atom(i).formalcharge=2;
    elseif strncmpi(atom(i).type,{'C'},1);atom(i).element={'C'};atom(i).formalcharge=0;
    else
        [atom(i).element]=atom(i).type;atom(i).formalcharge=0;
    end
end


nAtoms=length(atom);
Atom_section=cell(nAtoms,10);

atom=orto_atom(atom,Box_dim);

% Write the file
fid = fopen(filename_out, 'wt');
fprintf(fid, '%s\r\n','data_matlab_gen');
fprintf(fid, '\r\n');
timestamp = datetime;
fprintf(fid, '%s     %s\r\n','_audit_creation_date',timestamp);
fprintf(fid, '%s\r\n','_audit_creation_method   generated by the Matlab ATOM scripts');
fprintf(fid, '\r\n');
fprintf(fid, '\r\n');
fprintf(fid, '%-22s   %8.4f\r\n','_cell_length_a',a);
fprintf(fid, '%-22s   %8.4f\r\n','_cell_length_b',b);
fprintf(fid, '%-22s   %8.4f\r\n','_cell_length_c',c);
fprintf(fid, '%-22s   %8.4f\r\n','_cell_angle_alpha',alfa);
fprintf(fid, '%-22s   %8.4f\r\n','_cell_angle_beta',beta);
fprintf(fid, '%-22s   %8.4f\r\n','_cell_angle_gamma',gamma);
fprintf(fid, '\r\n');
fprintf(fid, '%s\r\n','_symmetry_space_group_name_H-M     P1');
fprintf(fid, '%s\r\n','loop_');
fprintf(fid, '%s\r\n','_space_group_symop_operation_xyz');
fprintf(fid, '  %s\r\n','x,y,z');
fprintf(fid, '\r\n');
fprintf(fid, '%s\r\n','loop_');
fprintf(fid, '%s\r\n','_atom_site_type_symbol');
fprintf(fid, '%s\r\n','_atom_site_label');
fprintf(fid, '%s\r\n','_atom_site_fract_x');
fprintf(fid, '%s\r\n','_atom_site_fract_y');
fprintf(fid, '%s\r\n','_atom_site_fract_z');
for i = 1:nAtoms
    Atom_section(1:5) = [atom(i).element, atom(i).type, atom(i).xfrac, atom(i).yfrac, atom(i).zfrac];
    fprintf(fid,'%-5s%-5s%10.5f%10.5f%10.5f\r\n',Atom_section{1:5});
end
fprintf(fid, '\r\n');

fclose(fid);

% for i = 1:nAtoms
%     Atom_section(1:13) = ['ATOM  ', atom(i).index, atom(i).type, atom(i).resname, 'A',atom(i).molid, atom(i).x, atom(i).y, atom(i).z,1,1,atom(i).element,atom(i).formalcharge];
%     fprintf(fid,'%-6s%5i %-4s %3s %1s%4i    %8.3f%8.3f%8.3f%6.2f%6.2f          %2s%2i\n',Atom_section{1:13});
% end
% assignin('caller','Bond_index',Bond_index);
% assignin('caller','Angle_index',Angle_index);