For the same setup of rigid body, I was trying to calculate phonon eigenvectors (input attached) for a unit cell. I calculated these at 50 k-points. But I was confused of some outputs in .eig file (output attached),
GULP seem to be recognizing only the 4 Na atoms and dumps the coordinates of only these atoms at the beginning of the file. Shouldn’t this include COM coordinates of remaining 4 superoxide molecules too?
From first k point results in eig file, I understand that, for each mode there are total of 8 eigenvectors with four corresponding to Na atoms and 4 eigenvectors corresponding to four superoxide molecules. But from second k-point, for same modes, there seem to be more than 8 eigenvectors (i.e. more number of lines are dumped). Can you please explain why are there are more eigenvectors and what do they correspond to for each mode from second k-point? or Is this some kind of output formatting issue in GULP? Would be great if you can help me out in understanding this output format.
PS: In one unit cell there are 4 Na atoms and 4 superoxide molecule (i.e. 8 oxygens and 4 He atoms)
There is an omission in that the coordinates of the rigid molecule aren’t printed in the .eig file, though since it’s a GULP format I guess I can define it to include (or not include) anything I like. The coordinates should be available elsewhere in the output.
All the modes for the first K point are at the gamma point and so everything is real, whereas from the second K point onwards then they’re complex and so the amount of information doubles. I guess you could pad the first K point with zeros to make them all look the same. For rigid molecules the eigenvectors are given in terms of the COM motion and then the 3 rotations which is why there are 6 numbers (at gamma) for the molecules, versus 3 for the atoms.
Hope that helps,
Julian
Thanks for the confirmation! Yeah that was helpful. Now I’m able to get some reasonable results (i.e. DOS and eigenvectors) for NaO2 supercell with perfect crystal structure.
But the problem still persists when I try to do the same calculations (i.e. DOS and eigenvector decomposition) for higher temperatures in disordered regime. I seem to be getting really high mode frequencies/imaginary. I followed these steps - I used equilibrated NVE snapshot of a higher temperature (say 300 K where the system is completely disordered) as a new input structure to the GULP. I accordingly dumped the image flags from LAMMPS and updated “connect” command too. I kept all other settings to be the same.
My intuition - GULP performs LD calculations for a given structure. So I thought, to perform temperature dependent DOS and eigenvector calculations,
(i) First, I run a NVT (fix to a target temperature T) and follow it with NVE ensemble in LAMMPS until equilibration for each ensemble.
(ii) Use snapshots from equilibrated NVE ensemble as input structure/coordinates to GULP and perform the same calculations as we did for perfect crystal.
(iii) Repeat the above process with multiple equilibrated snapshots from NVE and calculate averaged values of DOS and eigenvectors.
But there seems to be something clearly wrong with my idea as I’m getting wrong mode frequencies/vectors.
So it would be really helpful if you can share you thoughts on this.
PS: I’ve attached my GULP input file which uses NVE snapshot at 300 K as input structure.
Hi Erick,
It sounds like you’re trying to do something close to free energy minimisation in lattice dynamics & so it’s worth reading the papers on this topic since the same issues crop up. The key point is that you should only compute phonon frequencies in the harmonic approximation for a structure that it is at an internal energy minimum with respect to the coordinates, otherwise the frequencies will be meaningless and could easily be imaginary. So you can use NVT to find the cell parameters appropriate to the temperature, but then you need to minimise the coordinates within the fixed cell before computing phonons, rather than using snapshots. Even within this modified framework, what will happen is as the cell expands with temperature then you may well hit modes that go imaginary, depending on the symmetry of your minimised structure (of course when using quenched MD then you should have broken this and could be fine), as the structure will lower symmetry as the temperature increases to maximise the entropy. This is the usual break down of quasiharmonic free energy minimisation with increasing temperature due to anharmonic effects, and happens at low temperatures for molecular crystals and soft materials.
Regards,
Julian
Thanks for the detailed discussion! Yes I’m trying to do something similar to free energy minimization in lattice dynamics.
But, from manual, I realized that GULP might not be having this feature for rigid bodies. I think, following error message for NaO2 system confirms the same? Can you please let me know if my understanding is correct; if so, can you please suggest any possible circumvent to be implemented in GULP (or any other)
PS: I’m currently using perfect crystal structure as a test run for free minimization. But would ideally try with snapshot/configuration at finite temperature obtained from LAMMPS. I’ve also attached the input file for following message.
You are correct - there are no third derivatives for rigid bodies at present (and not likely to be for a while as it’s not a trivial exercise!). For high symmetry systems you can always do a manual bit of free energy minimisation by varying the cell parameter, optimising at constant volume and computing the free energy at the temperature of interest. In your case, since you have a cubic cell (I think), this is trivial since it’s just a matter of a 1-D scan and plotting the free energy vs cell parameter to find the minimum. GULP will allow you to compute the free energy for multiple temperatures per run, so you can find the temperature-dependent cell parameter from 1 set of runs vs cell parameter.
I think the problem was with the image flags. When I set all the image flags to 0 (like in the attached file) and use the unwrapped snapshots from LAMMPS output, I seem to be getting reasonable results. I think this could be due to the “unwrapped” keyword - since I’ve bonded superoxide ions, I meant to use this keyword to be specific about which images oxygen atoms bond to. But specifying non zero image flags might double correct the coordinate values and result in really high frequencies. Please let me know what you think
PS: Like you suggested previously, I even tried with temperature dependent lattice constants (obtained by relaxing the structure in NPT ensemble for a specified temperature). But that still was off (i.e. high frequencies).
The coordinates should be available elsewhere in the output.