camm / sassena

Sassena — X-ray and neutron scattering calculated from molecular dynamics trajectories using massively parallel computers
1 stars 4 forks source link

Update input config file documentation #10

Open jmborr opened 7 years ago

jmborr commented 7 years ago

Include in the LYX file the changes to section <background><factor> to allow that scattering length density from an atom selection (see issue #4)

thamnos commented 7 years ago

Hello! Could you please give me a quick pointer on how to use this? I have a box of pure water and calculate the static structure factor -- as well as the coh & inc intermediate scattering functions. In the static structure factor I get this behaviour at low Q, i.e. spurious peaks below Q \sim 0.5\AA^{-1}: fq Could I get these to go away by setting something like <background>all</background>? Would I do that for the intermediate scattering functions, too? Thanks! Sebastian.

jmborr commented 7 years ago

I'm a Quasielastic guy so I don't know much about the devilish details of the coherent scattering. Is 2*pi/0.5 = 12 Angstroms about the size of your simulation box? You are likely seeing the spurious correlations caused by the periodic boundary conditions of the simulation box.

The configuration

selection will remove the structure factor of selection from the total structure factor. If you select all, then you'll end up with no structure factor at all! So as far as I know, there's no way you can remove that artificial peak at low Q. The configuration makes sense when you have a protein surrounded by water, and you want to remove the structure factor due to water in order to find out the structure factor due to the protein only. On Tue, Aug 22, 2017 at 11:03 AM, thamnos wrote: > Hello! > Could you please give me a quick pointer on how to use this? I have a box > of pure water and calculate the static structure factor -- as well as the > coh & inc intermediate scattering functions. In the static structure factor > I get this behaviour at low Q, i.e. spurious peaks below Q \sim 0.5\AA^{-1}: > [image: fq] > > Could I get these to go away by setting something like > all? Would I do that for the intermediate > scattering functions, too? > Thanks! > Sebastian. > > — > You are receiving this because you were assigned. > Reply to this email directly, view it on GitHub > , or mute > the thread > > . >
thamnos commented 6 years ago

Right, that's indeed not the kind of background correction I was looking for... I think you are right that these spurious peaks are due to the finite box. I've played a bit with Fourier transforms. I took the same "signal" (a 1D representation of scattering lengths) and displaced it on the y-axis: oscillating around 0, 1, or 2. As soon as the values don't oscillate around 0, we get the artefacts at low Q; basically the FT of the box. fft_exercise_001 (code below)

So I'm thinking maybe we should subtract the average scattering length from all scattering lengths so that they oscillate around 0? What would that do to the units of the resulting scattering function? How would one do this in the case of x-rays where the scattering lengths are Q-dependent?

from matplotlib.pyplot import *
from scipy import *

a0 = zeros(1024)
a1 = zeros(1024)
a2 = zeros(1024)
num = 128 # length of the box
rnd = random.random(num) # between 0 and 1, mean is 0.5
a0[:num] = rnd - 0.5
b0 = fft(a0)
a1[:num] = rnd + 0.5
b1 = fft(a1)
a2[:num] = rnd + 1.5
b2 = fft(a2)

f, (ax1, ax2) = subplots(2,1)
ax1.plot(a0, 'r-', label='0')
ax1.plot(a1, 'g-', label='1')
ax1.plot(a2, 'b-', label='2')
ax1.legend()
ax1.set_xlim(0,1024)
ax2.plot(abs(b0[:512]), 'r-', label='abs(FFT(0))')
ax2.plot(abs(b1[:512]), 'g-', label='abs(FFT(1))')
ax2.plot(abs(b2[:512]), 'b-', label='abs(FFT(2))')
ax2.set_yscale('log')
ax2.set_xlim(0,512)
ax2.legend()
show()
jmborr commented 6 years ago

I think you found the source of the problem. I looked at the definition of the radial distribution in Wikipedia https://en.wikipedia.org/wiki/Radial_distribution_function#Definition, and the structure factor from a scattering experiment does diverge at the origin (see discussion between equations 5 and 6) because of the average density.

Going back to the definition of S(Q) in neutron scattering,

​this thing will certainly go to zero if we subtract the average scattering length density to each scattering length, while at the same time remaining positive for all values of Q.

In the case of X-rays we have an average density that is Q-dependent, but at sufficiently low Q's it does tend to a constant value, so we end up with the same situation as in the neutron scattering case.

Sassena allows subtraction of a solvent background scattering length density when you have a solute surrounded by a solvent. In this case we don't have any solute and we want to remove the background scattering caused by the whole system. You could try:

system

and see if this works. .Jose

On Fri, Sep 1, 2017 at 11:55 AM, thamnos notifications@github.com wrote:

Right, that's indeed not the kind of background correction I was looking for... I think you are right that these spurious peaks are due to the finite box. I've played a bit with Fourier transforms. I took the same "signal" (a 1D representation of scattering lengths) and displaced it on the y-axis: oscillating around 0, 1, or 2. As soon as the values don't oscillate around 0, we get the artefacts at low Q; basically the FT of the box. [image: fft_exercise_001] https://user-images.githubusercontent.com/7515211/29977281-2df0a536-8f3d-11e7-95cd-92a5aaeca9d9.png (code below)

So I'm thinking maybe we should subtract the average scattering length from all scattering lengths so that they oscillate around 0? What would that do to the units of the resulting scattering function? How would one do this in the case of x-rays where the scattering lengths are Q-dependent?

from matplotlib.pyplot import from scipy import

a0 = zeros(1024) a1 = zeros(1024) a2 = zeros(1024) num = 128 # length of the box rnd = random.random(num) # between 0 and 1, mean is 0.5 a0[:num] = rnd - 0.5 b0 = fft(a0) a1[:num] = rnd + 0.5 b1 = fft(a1) a2[:num] = rnd + 1.5 b2 = fft(a2)

f, (ax1, ax2) = subplots(2,1) ax1.plot(a0, 'r-', label='0') ax1.plot(a1, 'g-', label='1') ax1.plot(a2, 'b-', label='2') ax1.legend() ax1.set_xlim(0,1024) ax2.plot(abs(b0[:512]), 'r-', label='abs(FFT(0))') ax2.plot(abs(b1[:512]), 'g-', label='abs(FFT(1))') ax2.plot(abs(b2[:512]), 'b-', label='abs(FFT(2))') ax2.set_yscale('log') ax2.set_xlim(0,512) ax2.legend() show()

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/camm/sassena/issues/10#issuecomment-326617337, or mute the thread https://github.com/notifications/unsubscribe-auth/ABFVhkeIIzQh-_futH9nSXvhi9DMHnI5ks5seCjZgaJpZM4M2O-8 .