ISISNeutronMuon
[ENHANCEMENT] Add support for NPT/NσT calculations
Is your feature request related to a problem? Please describe. NPT simulations will require extra support for variable cell calculations.
NσT will require further consideration given that the algorithm presented in: https://pubs.acs.org/doi/epdf/10.1021/acs.jctc.3c00308?ref=article_openPDF seems to suggest it only applies to variable cell-lengths in potentially orthogonal systems.
Describe the solution you'd like Full NPT support in the first instance.
Current jobs
- [ ] AngularCorrelation
- [ ] AreaPerMolecule
- [ ] AverageStructure
- [ ] CenterOfMassesTrajectory
- [ ] CoordinationNumber
- [ ] CroppedTrajectory
- [ ] CurrentCorrelationFunction
- [x] Density
- [ ] DensityOfStates
- [ ] DensityProfile
- [ ] DipoleAutoCorrelationFunction
- [ ] DistanceHistogram
- [ ] DynamicCoherentStructureFactor
- [ ] DynamicIncoherentStructureFactor
- [ ] Eccentricity
- [ ] ElasticIncoherentStructureFactor
- [ ] GaussianDynamicIncoherentStructureFactor
- [ ] GeneralAutoCorrelationFunction
- [ ] GlobalMotionFilteredTrajectory
- [ ] Infrared
- [ ] JobStatus
- [ ] McStasVirtualInstrument
- [ ] MeanSquareDisplacement
- [ ] MolecularTrace
- [ ] NeutronDynamicTotalStructureFactor
- [ ] OrderParameter
- [ ] PairDistributionFunction
- [ ] PositionAutoCorrelationFunction
- [ ] PositionPowerSpectrum
- [ ] RadiusOfGyration
- [ ] RigidBodyTrajectory
- [ ] RootMeanSquareDeviation
- [ ] RootMeanSquareFluctuation
- [ ] SolventAccessibleSurface
- [ ] StaticStructureFactor
- [ ] StructureFactorFromScatteringFunction
- [ ] Temperature
- [ ] TrajectoryEditor
- [ ] UnfoldedTrajectory
- [ ] VanHoveFunctionDistinct
- [ ] VanHoveFunctionSelf
- [ ] VelocityAutoCorrelationFunction
- [ ] Voronoi
- [ ] XRayStaticStructureFactor
Describe alternatives you've considered Raise warning in GUI about variable cell. Rely on "average cell" for computing parameter.
Additional context N/A
To correct the CCF and DCSF calculations for NPT/NσT calculations we need to update the latticeqvector generator. Here is my reasoning for this which I will go through only for DCSF.
The equation for the intermediate scattering function (ignoring scattering lengths) is
F(q,t) = \frac{1}{N} \sum_{ij} \langle \exp[-i q \cdot r_{i}(0)] \exp[i q \cdot r_{j}(t)] \rangle
the sum over $i$ and $j$ are over all atoms in the system. When we are working with a periodic system them we can rewrite the equation above so that
F(q,t) = \frac{1}{N}\sum_{\alpha\beta} \sum_{kl} \langle \exp[-i q \cdot (r_{k}(0) + n_{\alpha}(0))] \exp[i q \cdot (r_{l}(t) + n_{\beta}(t) )] \rangle
where the sum over $k$ and $l$ go over all atoms in the unit cell and the sum over $\alpha$ and $\beta$ go over all lattice vectors. Lets rewrite the above equation.
F(q,t) = \frac{1}{N} \sum_{\alpha\beta} \sum_{kl} \langle \exp[-i q \cdot r_{k}(0)] \exp[i q \cdot r_{l}(t)] \exp[-i q \cdot (n_{\alpha}(0) - n_{\beta}(t) )] \rangle
Now we can select specific qvectors so that we can avoid running the summation over all lattice vectors. Specifically we want
\exp[-i q \cdot (n_{\alpha}(0) - n_{\beta}(t) )] = 1
for all $\alpha$ and $\beta$. For NVT when the lattice vectors are constant in time we can just select a reciprocal lattice vector since,
\exp[-i q \cdot (n_{\alpha} - n_{\beta} )] = \exp[-2 (N - M) \pi i ] = 1
here $q$ is a reciprocal lattice vector and $N$ and $M$ are some integers. This is what the latticeqvector generator does in MDANSE.
For NPT/NσT this is more complicated, we need to find qvectors such that
\exp[-i q \cdot (n_{\alpha}(0) - n_{\beta}(t) )] = 1
for all $\alpha$ and $\beta$, which doesn't look possible in general.
Here is one proposal to correct the DCSF and CCF calculations for NPT/NσT trajectories. For DCSF we want to calculate the following
F(q,t) = \frac{1}{N} \sum_{\alpha\beta} \sum_{kl} \langle \exp[-i q \cdot (r_{k}(0) + n_{\alpha}(0))] \exp[i q \cdot (r_{l}(t) + n_{\beta}(t) )] \rangle
but we want to avoid running the summations over the lattice vectors since it's an infinite sum. I propose the following scheme for NPT/NσT. We instead calculate the following function.
F(q,q',t) =\frac{1}{N} \sum_{\alpha\beta} \sum_{kl} \langle \exp[-i q \cdot (r_{k}(0) + n_{\alpha}(0))] \exp[i q' \cdot (r_{l}(t) + n_{\beta}(t) )] \rangle
Now we can set $q$ and $q'$ so that they are the reciprocal lattice vectors for the lattices generated by $n_{\alpha}(0)$ and $n_{\beta}(t)$, respectively so that
\exp[-i q \cdot n_{\alpha}(0)] = \exp[-2 \pi N i] = 1
and
\exp[i q' \cdot n_{\beta}(t)] = \exp[2 \pi M i] = 1.
Now, what I think we need to do is calculate $F(q,q',t)$ but set the value of $q'$ so that it is as close to $q$ as possible (determined from distances or something), so that
q \approx q'
and
F(q,q',t) \approx F(q,q,t) = F(q,t).
The nice thing about this method is that if the user increases the simulation box size, the reciprocal lattice becomes increasingly dense, so that it will be easier to find a better approximation of $q = q'$.
For NPT/NσT, all jobs that use velocities obtained by numerical derivatives will have issues. This is because with NPT/NσT, velocities are not simply the time derivatives of the positions but include other things like the volume when, for example, the Anderson Barostat is used. See 5.8 in Statistical Mechanics: Theory and Molecular Simulation.
Here, I run an NPT simulation of liquid Argon (10 fs time step) and run some calculations. I dump velocities and compare results when we use the actual velocities and velocities from numerical derivatives. Here are the jobs that use the interpolation_order setting for velocities: Temperature, DOS, VACF and CCF.
Temperature
DOS
VACF