MATLAB Examples

# xrd_atom.m

## Examples

• [twotheta,intensity] = xrd_atom(filename)
• [twotheta,intensity] = xrd_atom(atom,Box_dim))
```function [twotheta,intensity] = xrd_atom(varargin) ```

## Various settings

```mode=0; % Activate sigma_star, DIV, surface roughness FWHM=.5; % Specify the full width at half maximum lambda=1.54187; % Ångstrom exp_twotheta=2:.05:60; % The twotheta range of interest Bfactor=2; % Debye-Waller factor Ångstrom Lorentzian_factor=1; %Enter the fraction of calculated pattern you would like to have described by a lorentzian function vs. a gaussian function Sample_length = 4; % cm Gonio_radius = 24; % cm Div_slit = .1; % Divergence slit setting, 0 for automatic? roughness = 0; % Surface roughness % misalignment=0.0; % Not yet implemented % alfa_strain = 0; % Not yet implemented sigma_star=45; % Reynolds 00l mean preferred orientation, 45 [deg] is random, 1 [deg] is the opposite RNDPWD = 0; % Random powder % mu_star=45; % Not yet implemented L_type='normal'; % Lorent polarization type, else 'Reynolds'; S1=2.3;S2=2.3; % Primary and secondary Soller slit , in deg ```

## 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=0; preferred_l=1; ```

## Fetch either a .pdb|.gro file or use an atom struct with its Box_dim

```if nargin==1 filename=varargin{1}; atom=import_atom(filename); else atom=varargin{1}; Box_dim=varargin{2}; end ```

## Import the structure into an atom struct

```atom=wrap_atom(atom,Box_dim); % atom=replicate_atom(atom,Box_dim,[1 1 1]); % If we want to replicate the struct ```
```frac=orto_atom(atom,Box_dim); % frac = round_atom(frac,Box_dim,3,'orto'); frac=element_atom(frac); atom_type=[frac.type];% atom_type(2:length(atom_type)); x=[frac.xfrac]';y=[frac.yfrac]';z=[frac.zfrac]'; occupancy=ones(numel(x),1); ```

## Set the unit cell parameters

```if size(Box_dim,2) == 9 lx=Box_dim(1); ly=Box_dim(2); lz=Box_dim(3); xy=Box_dim(6); xz=Box_dim(8); yz=Box_dim(9); elseif size(Box_dim,2) == 3 lx=Box_dim(1); ly=Box_dim(2); lz=Box_dim(3); xy=0; xz=0; yz=0; end a=lx; b=(ly^2+xy^2)^.5; c=(lz^2+xz^2+yz^2)^.5; alfa=rad2deg(acos((ly*yz+xy*xz)/(b*c))); beta=rad2deg(acos(xz/c)); gamma=rad2deg(acos(xy/b)); ```

```%%%%% From MOF-FIT %%%%%% alfa_rad=alfa*pi/180; beta_rad=beta*pi/180; gamma_rad=gamma*pi/180; ```

## 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; hkl=[h' k' l']; ```

## Determine the two theta values for each h,k,l

```V_cell=a*b*c*(1-cos(alfa_rad)^2-cos(beta_rad)^2-cos(gamma_rad)^2+2*cos(alfa_rad)*cos(beta_rad)*cos(gamma_rad))^0.5; one_over_dhkl=1/V_cell.*... (h.^2*b^2*c^2*sin(alfa_rad)^2+... k.^2*a^2*c^2*sin(beta_rad)^2+... l.^2*a^2*b^2*sin(gamma_rad)^2+... 2.*h.*k*a*b*c^2*(cos(alfa_rad)*cos(beta_rad)-cos(gamma_rad))+... 2*k.*l*a^2*b*c*(cos(beta_rad)*cos(gamma_rad)-cos(alfa_rad))+... 2*h.*l*a*b^2*c*(cos(alfa_rad)*cos(gamma_rad)-cos(beta_rad))).^(0.5); one_over_dhkl=real(one_over_dhkl); %%%%% End MOF-FIT %%%%%% ```

## Order the hkl after decreasing d-spacings and increasing two_theta

```d_hkl=1./one_over_dhkl; twotheta_rad=2.*asin(one_over_dhkl*lambda/2); [d_hkl,hkl_order]=sort(d_hkl,'ascend'); twotheta_rad=twotheta_rad(hkl_order); hkl=hkl(hkl_order); h=h(hkl_order); k=k(hkl_order); l=l(hkl_order); ```

## Calculate the structure factor for each reflection

```two_theta_disc=real(180*twotheta_rad/pi); n=1; structure_factor=0; while(n<=length(atom_type)) scattering_factor = atomic_scattering_factors(atom_type(n),lambda,two_theta_disc,Bfactor); 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)); twotheta=real(180*twotheta_rad/pi); ```

## Construct the calculated pxrd pattern by adding a lorentzian fraction to a gaussian fraction

Gaussian part

```n=1;gauss_component=0;c_=FWHM/(2*(2*log(2))^0.5); while(n<=length(twotheta)) temp_gauss=F_squared(n).*exp(-(exp_twotheta-twotheta(n)).^2/(2*c_^2)); gauss_component=gauss_component+temp_gauss; n=n+1; end gauss_part=gauss_component/max(gauss_component); % Lorentzian part n=1;lorentz_component=0; while(n<length(twotheta)) temp_lorentz=F_squared(n)./((exp_twotheta-twotheta(n)).^2+(0.5*FWHM)^2); lorentz_component=lorentz_component+temp_lorentz; n=n+1; end lorentz_part=lorentz_component/max(lorentz_component); % Mix together the Lorentzian and Gaussian parts in the ratio specified by eta to generate the pseudo-Voigt function intensity=Lorentzian_factor*lorentz_part+(1-Lorentzian_factor)*gauss_part; ```

## Lorentz factor

```S_bar=((S1/2)^2+(S2/2)^2)^.5; Q=S_bar./(2*2^0.5*sin(exp_twotheta/2*pi()/180)*sigma_star); PSI=erf(Q)*(2*pi())^.5/(2*sigma_star*S_bar)-2*sin(exp_twotheta/2*pi()/180)/S_bar^2.*(1-exp(-Q.^2)); ```

## Lorentz * Polarization factors

```Lorentz=(1+cos(exp_twotheta*pi()/180).^2); SingXtalLorentz=sin(exp_twotheta/2*pi()/180); if strcmp(L_type,'Reynolds') LP=Lorentz./ ( sin(exp_twotheta*pi()/180).*(sin(exp_twotheta/2*pi()/180)).^0.8); else LP=Lorentz./SingXtalLorentz.*PSI; end LP_random = Lorentz./(sin(exp_twotheta/2*pi()/180)) * 1./sin(exp_twotheta*pi()/180); if RNDPWD == 1 LP=LP_random; end ```

## Divergence slit

```if Div_slit>0 DIV=Sample_length/(Gonio_radius*Div_slit*pi()/180).*sin(exp_twotheta./2*pi()/180); DIV(DIV>1)=1; else DIV=[sin(exp_twotheta/2*pi()/180)]; end ```

## Surface roughness

```SR=0.5*(1+(sin((exp_twotheta/2-roughness)*pi()/180)./sin((exp_twotheta/2+roughness)*pi()/180))); if mode==1 intensity=SR.*DIV.*LP.*intensity; intensity=intensity-min(intensity); intensity=real(intensity/max(intensity(1:floor(15/((exp_twotheta(end)-exp_twotheta(1))/length(exp_twotheta)))))); else intensity=intensity.*(1+cos(exp_twotheta*pi/180).^2)./(8*sin(exp_twotheta*pi/180/2).^2.*cos(exp_twotheta*pi/180/2)); intensity=real(intensity/max(intensity)); end ```

## Plot the results

```hold on; plot(exp_twotheta,intensity); ```