camm / sassena

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

periodic boundary conditions #16

Open thamnos opened 7 years ago

thamnos commented 7 years ago

Hello,

I am getting worried that the periodic boundary condition / minimum image convention might not be taken into account. I have a box of water (cubic, 66AA side length) and calculate the incoherent intermediate scattering function I(Q,t). For Q below \sim 0.5\AA^{-1}, it does not decay to 0. My trajectory is a LAMMPS-generated DCD file where the atom positions are wrapped to stay within the box. The behaviour of I(Q,t) seems to indicate that particles are trapped in a finite volume, i.e. the box. Looking through the code, I can't see anywhere where sassena unwraps the trajectory. I think this would be necessary, right? Not only are atoms "confined to the simulation box" now, also the particles close to the border of the box perform huge jumps from one end to the other whilst actually only moving a tiny distance across the border. For pair correlations it's probably a bit trickier since the structure changes if the molecules "flow apart", so it's not possible to just unwrap the coordinates (on the contrary, if given unwrapped coordinates, one would have to wrap them first, I guess). But are we taking correlations between atoms at opposite ends of the simulation box (and therefore close to each other across the box boundary) properly into consideration?

Cheers, Sebastian.

jmborr commented 7 years ago

Sassena does not unwrap the trajectory. For incoherent scattering you must unwrap the trajectory first. For coherent scattering, there's no easy way to get rid of the periodic boundary conditions. As this paper suggest http://orproxy.lib.utk.edu:2053/science/article/pii/S037838120400490X, you could compute the structure factor only for reciprocal lattice vectors, and then average for those vectors with the same modulus. You can pass a file to Sassena specifying the wavectors for which you want to calculate the structure factor.

On Tue, Aug 22, 2017 at 11:43 AM, thamnos notifications@github.com wrote:

Hello,

I am getting worried that the periodic boundary condition / minimum image convention might not be taken into account. I have a box of water (cubic, 66AA side length) and calculate the incoherent intermediate scattering function I(Q,t). For Q below \sim 0.5\AA^{-1}, it does not decay to 0. My trajectory is a LAMMPS-generated DCD file where the atom positions are wrapped to stay within the box. The behaviour of I(Q,t) seems to indicate that particles are trapped in a finite volume, i.e. the box. Looking through the code, I can't see anywhere where sassena unwraps the trajectory. I think this would be necessary, right? Not only are atoms "confined to the simulation box" now, also the particles close to the border of the box perform huge jumps from one end to the other whilst actually only moving a tiny distance across the border. For pair correlations it's probably a bit trickier since the structure changes if the molecules "flow apart", so it's not possible to just unwrap the coordinates (on the contrary, if given unwrapped coordinates, one would have to wrap them first, I guess). But are we taking correlations between atoms at opposite ends of the simulation box (and therefore close to each other across the box boundary) properly into consideration?

Cheers, Sebastian.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/camm/sassena/issues/16, or mute the thread https://github.com/notifications/unsubscribe-auth/ABFVhh8Bwsvy7agcGYAt4LvKI80wNLTIks5savc3gaJpZM4O-0Ml .

thamnos commented 7 years ago

When you say "coherent scattering", do you mean the coherent intermediate scattering function I(Q,t) or the coherent static structure factor S(Q)?

jmborr commented 7 years ago

The problem affects both, since S(Q) is just I(Q,t=0), right? codecogseqn The second term involves subtraction of two real-lattice vectors, so it's just zero unless Q is a reciprocal-lattice vector.