function atom_type=breathing(theta) %%%%%breathing(theta) takes as an input the angle theta made between two %%%%%ligands meeting at the metal node. Also required are the name of the %%%%%coordinates file, text file of the powder diffraction pattern, desired %%%%%name of the output dummy.RES file, and unit cell parameters. Other %%%%%parameters that can be altered are the full-width of half maximum, %%%%%wavelength of the x-ray source, isotropic debye-waller parameter, %%%%%degree to which the calculated peak profile is Gaussian or Lorenzian, %%%%%maximum h,k,l values, and direction and degree of preferential %%%%%orientation. Any of the numerical quantities can be varied with %%%%%sliders by making them inputs. %%%%%%%Parameters for the user to change %Specify the file that will supply atom types, x, y, z, and occupancy. This %file should have the format: atom_type x y z occupancy. There should be no %extraneous lines/text. This does not distinguish between tabs and spaces. coordinates='ZnNDIH_coords.txt'; %Specify the file that will supply the experimental powder pattern. The %file should contain two columns - one with the two theta values and the %other with the intensities, in that order. experimental_pattern='pxrd.txt'; %%%%Specify the filename of the output dummy .RES file resfilename='ZnNDIbreathing.res'; %%%%Specify the full width at half maximum fwhm=0.1; %Wavelength of the x-ray source in angstroms lambda=1.54056; %isotropic debye-waller parameter - usually between 3.5 and 6.5 square %angstroms Biso=0; %%%Enter the unit cell parameters here a=24.5606; b=24.5606; c=7.1635; alpha=90; beta=90; gamma=90; %%%%End of unit cell input %Enter % of calculated pattern you would like to have described by a %lorentzian function vs. a gaussian function eta=0.5; %%%%%%Enter the maximum h, k, and l values you would like to calculate. The %%%%%%calculated peaks will be for -hmax <= h <= hmax; -kmax<= k <= kmax; %%%%%%-lmax <= l <= lmax hmax=10; kmax=10; lmax=10; %%%Enter the degree (number greater than or equal to 0) and direction of %preferential orientation pref=0; preferred_h=0; preferred_k=1; preferred_l=0; %%%%No more parameters for the user to change %%%%%%%%Program starts here %%%Converting to radians alpha=alpha*pi/180; beta=beta*pi/180; gamma=gamma*pi/180; theta=theta*pi/180; %Reads in the atom type with number appended, x, y, and z from the %specified file [atom_type_and_number,x2,y2,z2,occupancy]=textread(coordinates,'%s %f %f %f %f'); %%%%Input the symmetry elements present here x=[x2];y=[y2];z=[z2]; atom_type_and_number=[atom_type_and_number.'].'; occupancy=[occupancy]; x=[x2;0.5+x2;0.5-x2;0-x2;0.5+y2;y2;0.5-y2;0-y2;0-x2;0.5-x2;0.5+x2;0+x2;0.5-y2;0-y2;0.5+y2;0+y2]; y=[y2;0-y2;0.5-y2;0.5+y2;0.5+x2;0.5-x2;0+x2;0-x2;0-y2;0+y2;0.5+y2;0.5-y2;0.5-x2;0.5+x2;0-x2;0+x2]; z=[z2;0.5-z2;0+z2;0.5-z2;0-z2;0.5+z2;0.5+z2;0-z2;0-z2;0.5+z2;0-z2;0.5+z2;0+z2;0.5-z2;0.5-z2;0+z2]; %x=[x;0.0-x;0.0+x;0.0+x;0.0-x;0.0+x;0.0-x;0.0-x]; %y=[y;0.0+y;0.0-y;0.0+y;0.0-y;0.0-y;0.0+y;0.0-y]; %z=[z;0.5+z;0.0+z;0.5-z;0.0-z;0.5-z;0.0-z;0.5+z]; atom_type_and_number=[atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.',atom_type_and_number.'].'; occupancy=[occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy;occupancy]; %%%%%%%%%%%%Now trying to strip the number from the atom name atom_type={''}; for n=1:length(atom_type_and_number) current=char(atom_type_and_number(n)); %%%%Change the atom type+number from a cell to text if isempty(regexp(current,'\d')); atomtype_new=current; else atomtype_new=current(1:regexp(current,'\d')-1); end %%%%cut off the part of the name after the first appearance of a number cell={atomtype_new}; atom_type(length(atom_type)+1)=cell; %%%Append this new atom name to the list of atom types end atom_type=atom_type(2:length(atom_type)); %%%%%%%%%%%%%%Done stripping number from atom name %%%Change of unit cell parameters a and b corresponding to the new theta %%%value diagonal=(a^2+b^2)^(0.5); a=diagonal*sin(theta/2); b=diagonal*cos(theta/2); %%%%%%%%%%%%%%Setting up the different h,k,l values l_values_temp=repmat(-1*lmax:lmax,1,(2*kmax+1)*(2*hmax+1)); k_repeat_unit=repmat(-1*kmax:kmax,(2*lmax+1),1); k_values_temp=repmat(reshape(k_repeat_unit,1,(2*lmax+1)*(2*kmax+1)),1,(2*hmax+1)); h_values_temp=reshape(repmat(-1*hmax:hmax,(2*kmax+1)*(2*lmax+1),1),1,(2*hmax+1)*(2*kmax+1)*(2*lmax+1)); h_values_temp(hmax*(2*kmax+1)*(2*lmax+1)+kmax*(2*lmax+1)+lmax+1)=[]; k_values_temp(hmax*(2*kmax+1)*(2*lmax+1)+kmax*(2*lmax+1)+lmax+1)=[]; l_values_temp(hmax*(2*kmax+1)*(2*lmax+1)+kmax*(2*lmax+1)+lmax+1)=[]; h=h_values_temp; k=k_values_temp; l=l_values_temp; %%%%%%%%%%%%%%%Now we have all the h,k,l values with 0,0,0 taken out %%%%%%%%%%%%%%Determine the two theta values for each h,k,l V_cell=a*b*c*(1-cos(alpha)^2-cos(beta)^2-cos(gamma)^2+2*cos(alpha)*cos(beta)*cos(gamma))^0.5; one_over_dhkl=1/V_cell.*(h.^2*b^2*c^2*sin(alpha)^2+k.^2*a^2*c^2*sin(beta)^2+l.^2*a^2*b^2*sin(gamma)^2+2.*h.*k*a*b*c^2*(cos(alpha)*cos(beta)-cos(gamma))+2*k.*l*a^2*b*c*(cos(beta)*cos(gamma)-cos(alpha))+2*h.*l*a*b^2*c*(cos(alpha)*cos(gamma)-cos(beta))).^(0.5); two_theta=2.*asin(one_over_dhkl*lambda/2); %http://classes.uleth.ca/200401/chem4000a/Lecture10.pdf %%%%%%%%%%%%twotheta values now determined %%%%%%%%%%%%%%Calculate the structure factor for each reflection n=1; structure_factor=0; while(n<=length(atom_type)) %%%%read in the parameters for the scattering factor equation and give %%%%them names sf_parameters=scattering_factor_parameters(atom_type(n)); a1=sf_parameters(1); b1=sf_parameters(2); a2=sf_parameters(3); b2=sf_parameters(4); a3=sf_parameters(5); b3=sf_parameters(6); a4=sf_parameters(7); b4=sf_parameters(8); c_=sf_parameters(9); %%%%%Plug the parameters into the scattering factor equation scattering_factor=a1.*exp(-b1.*sin(two_theta/2).^2/lambda^2)+a2.*exp(-b2.*sin(two_theta/2).^2/lambda^2)+a3.*exp(-b3.*sin(two_theta/2).^2/lambda^2)+a4.*exp(-b4.*sin(two_theta/2).^2/lambda^2)+c_;%Check %%%%%%Account for vibration in the form of the debye-waller factor scattering_factor=scattering_factor.*exp(-Biso.*sin(two_theta/2).^2/lambda^2); %%%%%%%Determine the contribution of the current atom to the overall %%%%%%%structure factor structure_factor=structure_factor+scattering_factor.*occupancy(n).*exp(2*pi*1i.*(h*x(n)+k*y(n)+l*z(n))); n=n+1; end %%%%%square each term in the structure factor vector F_squared=structure_factor.*conj(structure_factor); %%%%%Correction for preferred orientation %%Angle between h,k,l and preferrential orientation direction theta_pref_orient=acos((h*preferred_h + k*preferred_k + l*preferred_l)./((h.^2+k.^2+l.^2).^0.5*(preferred_h^2+preferred_k^2+preferred_l^2)^0.5)); if theta_pref_orient>pi/2 theta_pref_orient=pi-theta_pref_orient; end F_squared=F_squared.*exp(pref*cos(2*theta_pref_orient)); %%%%%read in the experimental pattern [experimental_two_theta,unnormalized_experimental_intensity]=textread(experimental_pattern,'%f %f'); %%%%%%convert the calculated two theta values to degrees for easier %%%%%%plotting two_theta=180*two_theta/pi; %%%%%%Construct the calculated pxrd pattern by adding a lorentzian fraction %%%%%%to a gaussian fraction %%Gaussian part c_=fwhm/(2*(2*log(2))^0.5); n=1; unnormalized_gauss=0; while(n<=length(two_theta)) unnormalized_gauss_contribution=F_squared(n).*exp(-(experimental_two_theta-two_theta(n)).^2/(2*c_^2)); unnormalized_gauss=unnormalized_gauss+unnormalized_gauss_contribution; n=n+1; end gauss_part=unnormalized_gauss/max(unnormalized_gauss); %%Lorentzian part n=1; unnormalized_lorentz=0; while(n