Hello!
I have xyz data for a polymer that I simulated in lammps. I need to generate the x-ray scattering intensity profile for my system. Does anyone know how to do this? Here is a matlab script that ChatGPT kindly provided, but I don’t know if it’s right. I’m running another simulation that will allow me to compare with experimental data.
% Define atomic positions
positions = positions4; %(this is the dump file from lammps, brought into Matlab as a table and then table2array for just the xyz columns)
lambda = 0.15; %X-ray wavelength in nm
L = 50; % Largest dimension of the polymer system in nm
theta_max = 30; %Maximum half-scattering angle in degrees
% Define scattering vectors
q_min = 2 * pi / L;
q_max = (4 * pi / lambda) * sind(theta_max);
num_points = 10000; % number of q points
q_values = linspace(q_min, q_max, num_points); % 1D array of q values
% Compute the structure factor
S_q = structure_factor(positions, q_values);
% Calculate the intensity
I_q = abs(S_q).^2;
% Convert q to 2*theta
two_theta = 4 * asind((q_values * lambda) / (4 * pi));
% Plot the intensity vs q
figure;
plot(two_theta, I_q);
xlabel('2theta');
ylabel('Intensity');
title('X-ray Scattering Intensity');
grid on;
% Function to compute the structure factor
function S_q = structure_factor(positions, q_values)
N = size(positions, 1);
S_q = zeros(size(q_values), 'like', 1 + 1i); % Initialize complex array
for q_index = 1:length(q_values)
q = q_values(q_index);
sum_real = 0;
sum_imag = 0;
for j = 1:N
pos = positions(j, :);
qr = dot([q, 0, 0], pos); % q·r_j
sum_real = sum_real + cos(qr);
sum_imag = sum_imag + sin(qr);
end
S_q(q_index) = sum_real + 1i * sum_imag;
end
end
For reference, what I need is found in this paper, figure 9: (DOI 10.1021/ma9810565)
A sample dump file is attached.
Thanks for your time.
AAPoly366dump.20500000.dump (429.3 KB)
Kind regards,
Sean