ISISNeutronMuon / MDANSE

MDANSE: Molecular Dynamics Analysis for Neutron Scattering Experiments
https://www.isis.stfc.ac.uk/Pages/MDANSEproject.aspx
GNU General Public License v3.0
18 stars 4 forks source link

[ENHANCEMENT] Q vector generation: path in reciprocal space #446

Open MBartkowiakSTFC opened 1 month ago

MBartkowiakSTFC commented 1 month ago

Is your feature request related to a problem? Please describe. As far as I can tell, none of the q vector generators in MDANSE can create vectors along a path. (E.g. (0,0,0) -> (0.5,0,0) -> (0.5,0.5,0) ). Also, while I could use the DispersionLatticeQVectors to generate each of these segments separately, it seems that the HKL indices of the points are not written to the output.

Describe the solution you'd like A q vector generator should be added which creates vectors along an arbitrary path. The way vectors are written to the output file needs to be changed so that the HKL indices of the vectors are saved. Plotter has to be extended to plot the results as a function of q. (At the moment it only does |q|).

Describe alternatives you've considered Suggestions welcome.

Additional context N/A

MBartkowiakSTFC commented 3 weeks ago

In an MD simulation of a crystalline solid, the simulation box is most likely going to be a supercell. The HKL values in this situation should correspond to the unit cell, and not the supercell.

gonzalezma commented 3 weeks ago

Are you sure? I do not see how you could handle this easily (in most cases when you read a trajectory you will not know if it is a single unit cell or a nmp supercell) and nothing ensures that the MD will maintain the symmetry of the unit cell. I would imagine that if you have e.g. a 4x4x4 supercell, the user will interpret the vector (4,0,0) as the vector (1,0,0) of the unit cell. But I don't work with crystals, so you decide.

MBartkowiakSTFC commented 3 weeks ago

You are right that this will require some effort to get it right. So, the more problems with this idea you point out now, the better we can prepare our implementation.

In my own simulations of crystals, I typically specify the input structure as a single cell, and then use the MD engine commands to replicate it, as in: https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/CELL.html#CP2K_INPUT.FORCE_EVAL.SUBSYS.CELL.MULTIPLE_UNIT_CELL or https://docs.lammps.org/replicate.html and then the three numbers that follow the command specify the multiplicity of the supercell. (The random initialisation of velocities will ensure that the atoms move in different directions, and each copy of the unit cell is in a different configuration starting from frame 2.) The simplest approach from the MDANSE side would be to ask the user to input those three numbers again. Or leave them as 1,1,1 if they prefer the HKL value to be expressed relative to the total cell size. (Which could be valid in case of symmetry breaking which leads to the unit cell size doubling, tripling, etc. - but I would expect the finite-temperature effects of MD simulations to make the structure more symmetrical when integrating over long enough time.)

General thoughts:

  1. Technically, we should only consider something a crystal if the average positions of atoms remain fixed (within tolerance) at one position. The mean square displacement calculation could be one way of confirming it.
  2. At each specific time frame, we should expect no symmetry in the system. That is, a single frame should normally be triclinic P1, since the atoms are all displaced. Only the average structure remains symmetrical.
  3. Experimentalists typically have their results in HKL coordinates of the unit cell, and would like to have a matching scale in the calculations. This is where the other issue connects to it: https://github.com/ISISNeutronMuon/MDANSE/issues/458 - this is the request from experimentalists who measured a specific cut in the reciprocal space, and would like to run the calculation for the same set of points.

So, maybe we need to add a convenience feature of MDANSE producing an output of the average positions of the atoms as a single structure file, so the users can use their favourite tool (spglib, phonopy?) to check the symmetry and find the primitive and Bravais unit cells. And this should also use the atom selection tool, in case there were guest atoms in the lattice undergoing diffusion which the user wanted to exclude from the averaging.