MATLAB Examples

# write_atom_pqr.m

• This function writes an pqr file from the atom struct
• Tested 15/04/2017
• Please report bugs to michael.holmboe@umu.se

## Examples

• write_atom_pqr(atom,Box_dim,filename_out)
```function write_atom_pqr(atom,Box_dim,filename,varargin) ```
```% Have you done this? nAtoms=size(atom,2) % if exist('[atom.charge]','var') && exist('[atom.radius]','var') % disp('Charges and radius exist') % else % if nargin > 3 % short_r=cell2mat(varargin(1)); % long_r=cell2mat(varargin(2)); % else % short_r=1.25; % long_r=2.25; % end % % if nargin>5; % ffname=varargin(3); % if nargin>6; % watermodel=varargin(4); % else % disp('Unknown watermodel, will try SPC/E') % watermodel='SPC/E'; % end % if strncmpi(ffname,'clayff',5); % clayff_param(sort(unique([atom.type])),watermodel); % Total_charge = check_clayff_charge(atom) % elseif strcmpi(ffname,'interface'); % clayff_param(sort(unique([atom.type])),watermodel); % Total_charge = check_interface_charge(atom) % else % disp('Unknown forcefield, will try clayff') % clayff_param(sort(unique([atom.type])),watermodel); % Total_charge = check_clayff_charge(atom) % end % end % atom = charge_atom(atom,Box_dim,ffname,watermodel,short_r,long_r); % atom = radius_atom(atom,ffname,watermodel); % end if abs(sum([atom.charge])) > 0.01; disp('Not charge neutral system') sum([atom.charge]) end if abs(sum([atom.charge])) > 0.00001; disp('Tweaking the charge') qtot=sum([atom.charge]); charge=num2cell([atom.charge]-qtot/nAtoms); [atom.charge]=deal(charge{:}); end indH= strncmpi([atom.type],'H',1); [atom(indH).radius]=deal(0.1); % .pqr format usually something like % Field_name Atom_number Atom_name Residue_name Chain_ID Residue_number X Y Z Charge Radius if regexp(filename,'.pqr') ~= false; filename = filename; else filename = strcat(filename,'.pqr'); end Atom_section=cell(nAtoms,10); fid = fopen(filename, 'wt'); fprintf(fid, '%s\n','COMPND UNNAMED'); fprintf(fid, '%s\n','AUTHOR GENERATED BY MATLAB'); % % % 1 - 6 Record name "CRYST1" % % % % % % 7 - 15 Real(9.3) a (Angstroms) % % % % % % 16 - 24 Real(9.3) b (Angstroms) % % % % % % 25 - 33 Real(9.3) c (Angstroms) % % % % % % 34 - 40 Real(7.2) alpha (degrees) % % % % % % 41 - 47 Real(7.2) beta (degrees) % % % % % % 48 - 54 Real(7.2) gamma (degrees) % % % % % % 56 - 66 LString Space group % % % % % % 67 - 70 Integer Z value % % % % % % % % % Example: % % % % % % 1 2 3 4 5 6 7 % % % 1234567890123456789012345678901234567890123456789012345678901234567890 % % % CRYST1 117.000 15.000 39.000 90.00 90.00 90.00 P 21 21 21 8 % % % CRYST1 31.188 54.090 20.000 90.00 90.00 90.00 P1 1 disp('Assuming P1 space group. Box can still be triclinic') if length(Box_dim)==3; fprintf(fid, '%6s%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %11s%4i\n','CRYST1', Box_dim(1:3), 90.00, 90.00, 90.00, 'P1', 1); elseif length(Box_dim)==9 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)); fprintf(fid, '%6s%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %11s%4i\n','CRYST1', a, b, c, alfa, beta, gamma, 'P1', 1); else disp('No proper box_dim information') end % COLUMNS DATA TYPE FIELD DEFINITION % ------------------------------------------------------------------------------------- % 1 - 6 Record name "ATOM " % 7 - 11 Integer serial Atom serial number. % 13 - 16 Atom name Atom name. %SKIPPED % 17 Character altLoc Alternate location indicator. % 18 - 20 Residue name resName Residue name. % 22 Character chainID Chain identifier. % 23 - 26 Integer resSeq Residue sequence number. %SKIPPED % 27 AChar iCode Code for insertion of residues. % 31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms. % 39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms. % 47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms. % 55 - 60 Real(6.2) occupancy Occupancy. % 61 - 66 Real(6.2) tempFactor Temperature factor. % 77 - 78 LString(2) element Element symbol, right-justified. % 79 - 80 LString(2) charge Charge on the atom. % .pqr format usually something like % Field_name Atom_number Atom_name Residue_name Chain_ID Residue_number X Y Z Charge Radius 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,{'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,{'Ow'},2);atom(i).element={'Ow'};atom(i).formalcharge=-2; elseif strncmpi(atom(i).type,{'Hw'},2);atom(i).element={'Hw'};atom(i).formalcharge=1; 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 % [atom.type]=atom.element; % ATOM 1 Na Na A 1 9.160 6.810 1.420 1.00 1.00 Na 0 % ATOM 1 Si MMT A 1 2.140 8.380 2.710 1.00 0.00 S % .pqr format usually something like % Field_name Atom_number Atom_name Residue_name Chain_ID Residue_number X Y Z Charge Radius for i = 1:nAtoms Atom_section = ['ATOM ', atom(i).index, atom(i).type, atom(i).resname, 'A',atom(i).molid, atom(i).x, atom(i).y, atom(i).z,atom(i).charge,atom(i).radius,atom(i).element,atom(i).formalcharge]; %sprintf('%-6s%5i %4s %3s %1s%4i %8.3f%8.3f%8.3f %8.5f%8.5f %2s%2i\n',Atom_section{1:length(Atom_section)}); fprintf(fid,'%-6s%5i %4s %3s %1s%4i %8.3f%8.3f%8.3f %8.5f%8.5f %2s%2i\n',Atom_section{1:length(Atom_section)}); end % Write conect records if nargin>3; if nargin > 4 short_r=cell2mat(varargin(1)); long_r=cell2mat(varargin(2)); else short_r=1.25; long_r=1.25; end atom=bond_angle_atom(atom,Box_dim,short_r,long_r); assignin('caller','Bond_index',Bond_index); assignin('caller','Angle_index',Angle_index); assignin('caller','nBonds',nBonds); assignin('caller','nAngles',nAngles); B=[Bond_index(:,1:2); Bond_index(:,2) Bond_index(:,1)]; b1=sortrows(B); for i=min(b1(:,1)):max(b1(:,1)) ind=find(b1(:,1)==i); b2=b1(ind,2); fprintf(fid,'CONECT%5i%5i%5i%5i%5i%5i%5i',[i;b2]); fprintf(fid,'\n'); if mod(i,100)==1 i-1 end end fprintf(fid,'MASTER %5i%5i%5i%5i%5i%5i%5i%5i%5i%5i%5i%5i\n',[0 0 0 0 0 0 0 0 nAtoms 0 i 0]); fprintf(fid,'END'); fprintf(fid,'\n'); else fprintf(fid, '%s\n','TER'); fprintf(fid, '%s\n','ENDMDL'); end fclose(fid); disp('.pqr structure file written') % assignin('caller','atom',atom); % Total_charge=sum([atom.charge]) ```