We are using GULP to study a zone-boundary phase transition in a molecular crystal. As a sanity check, we performed calculations in a twofold supercell, so that the modes that, in the true cell, occur at either the zone centre or the zone boundary should occur in the supercell at the zone centre.
It is straightforward to separate these by eye into zone-centre and zone-boundary modes. Indeed, when we compare against the single-cell calculation, we find that the eigenvalues match exactly, and divide between the zone centre and zone boundary exactly as we expect.
Unfortunately, though, the eigenvectors do not appear to match. Indeed, the eigenvectors in the smaller cell appear to correspond to motion that, given our potentials, would have a far higher energetic cost than the eigenvalues suggest, so we believe that somehow these are incorrect, although the input files are exactly analogous.
We’d be very grateful for any suggestions you might have to help fix this! Full details and input/output files are given below. We have until recently been using v 5.1, but we have just re-run on 6.1 and confirmed that our problem is still as described below. Many thanks in advance.
We are calculating the phonons of the plastic crystal adamantane with GULP. Adamantane (C10H16) consists of rigid, almost spherical molecules and has two solid phases. At low temperatures, it has the simple tetragonal space group P4Ě…21c; at high temperatures, the face-centred cubic space group Fm3Ě…m.
In order to do a direct comparison between phonons of the two phases, we transformed the tetragonal phase to (pseudo)cubic coordinates and calculated the phonons in GULP. The tetragonal cell has 2 molecules in the unit cell, while the pseudocubic one has 4. As a sanity check, we wanted to confirm whether the these modes are the same as the ones produced by GULP of the same system in the tetragonal setting. Specifically, the phonon modes in the tetragonal setting at the point M=[0.5, 0.5, 0.0] should map back to modes at Gamma in the cubic setting.
Attached are the phonons at Gamma in the cubic setting, and those at Gamma and M in the tetragonal setting. They are produced with GULP’s output castep_phon option. The frequencies at Gamma in the cubic setting are indeed reproduced at points Gamma and M in the tetragonal setting. However, the eigenvectors of the tetragonal M point do not seem to line up with the equivalent ones at the cubic Gamma point.
For example, mode 1 of the tetragonal M point should be equivalent to mode 7 of the cubic Gamma point. The frequencies agree, but visualisation in Jmol shows that in the cubic setting it’s a molecular librational mode, while in the tetragonal setting it’s a stretching mode that involves distorting the molecules.
Our intuition is that the eigenvectors of the modes at M in the tetragonal cell are incorrect: we have deliberately set the harmonic stretch and bend force constants very high, so at these energies there should only be rigid body movements. This is corroborated by the fact that there are clearly 12 (in the tetragonal cell) and 24 (in the cubic cell) low-frequency modes, which account for the 3 translational and 3 rotational modes for each molecule. The internal modes only start around 900 cm–1. On the other hand, the eigenvectors in the cubic cell are believable.
We do not understand why the tetragonal cell would produce such low-energy stretch modes; but more importantly, it is strange that the modes in the two different settings have the same frequencies, but completely different eigenvectors. We hoped it was an issue with the visualisation in Jmol, but a double-check with different software produced the same disparity.
This doesn’t seem to be an inconsistency in the phase factor between different cells, since the molecule in the centre of the cell is distorted in the same way as the one on the vertex. Moreover, it doesn’t seem to be due to the notorious different possible definitions of r, since the workflow from GULP to CASTEP-style phonons to Jmol is one that we use frequently and that usually works.
One possible reason for the discrepancy that should be checked is in regard to how the eigenvectors from GULP are used, since there are 2 conventions on the form of eigenvectors. The eigenvectors in GULP literally come from the matrix that diagonalises the dynamical matrix (i.e. the phased Hessian divided by the masses). However, in some places the eigenvectors are reweighted by the inverse square root of the mass (and renormalised). If the wrong eigenvectors are assumed when animating the motion then this would turn a rigid body motion into something looks like the internal degrees are being coupled to external ones. For example, consider a dummy molecule of HI (chosen for extreme mass difference) with a stiff spring constant (similar to your adamantane case) with the GULP input as follows:
phon noden eigen
cart
H core 0.0 0.0 0.0
I core 2.0 0.0 0.0
harm
H core I core 100.0 2.0 0.0 3.0
The resulting eigenvectors are:
Frequencies (cm-1) and Eigenvectors :
Frequency -0.0000 0.0000 0.0000 0.0000 0.0000 5209.4189
1 x -0.088857 -0.000000 0.000000 -0.000000 0.000000 -0.996044
1 y 0.000000 1.000000 0.000000 -0.000000 0.000000 0.000000
1 z 0.000000 0.000000 1.000000 -0.000000 0.000000 -0.000000
2 x -0.996044 0.000000 0.000000 -0.000000 0.000000 0.088857
2 y -0.000000 0.000000 0.000000 1.000000 0.000000 0.000000
2 z -0.000000 0.000000 0.000000 0.000000 1.000000 0.000000
If we look at the first mode, this has a frequency of zero, and is a rigid translation in x. However, if animated, it would appear that the H-I bond is stretching as well due to the difference in the eigenvectors for H and I. If we divide both values by the square root of the mass then you get that both eigenvector components are -0.088416 (prior to renormalisation).
These values would then make sense in terms of the motion.
Bottom line is that I think Jmol expects the eigenvectors divided by the square root of the masses, rather than the pure eigenvectors.
Hope that helps.
Julian
PS I’ll try to add an option that switches the values output between the 2 conventions asap so that you can test.
[Reverted since the description of this test system might prove useful]
I do appreciate that this is a complicated system. It wasn’t at first obvious how to condense it into a simpler test case but I think I have now managed to do so. All files are in the same folder as before; the new ones have names beginning test.
You will see that in this diatomic model system the lowest-frequency phonon (–1.58 cm–1) at (0, 0, 0.5) has a substantial stretching component, which again shouldn’t be possible given the high bond-stretching frequency.
Once again, performing the calculation in a supercell gives the same frequencies but more sensible eigenvectors.
As one might expect, the dynamical matrix has large components in the direction of the harmonic bonds (i.e., in the 3,3; 3,6; 6,3; and 6,6 entries) and small ones elsewhere. I thought we might therefore be running into rounding error, making the eigenvectors arbitrary. However, increasing the Buckingham A by a factor of 100 and/or decreasing the harmonic k by a factor of 10 gives similar results: naturally the frequency changes a little, but there is still a low-frequency mode at (0, 0, 0.5) with a large stretching component.
For what it’s worth, I have also noted that, in the original adamantane cell, using the proper symmetry (as opposed to the P1 cell I showed in my previous post) makes no difference to the anomalous results.
Dear Anthony,
Apologies - I thought I’d replace my old reply (since it wasn’t actually a useful solution) with what I think the real answer is. I didn’t realise you were about to post a reply yourself and so my change has messed up the context for your reply - sorry.
Regards,
Julian
Aha! That makes a lot of sense. If you are able to add a quick flag to change this, that would be great, but otherwise it should be pretty easy for me to change at this end too.
I should be able to add a keyword/option to change this tomorrow (it’s evening here now) and update GULP-6.1. When this is done I’ll post a note here to let everyone know…
PS It wouldn’t make a difference when all the atoms are C, so I’ll also look at this example once I’ve made the change.
Great, many thanks again. My model does show the same effect with every atom being the same, so I think it would be worth checking when you have a moment.
In fact, the more I think about this the less I am convinced that this is the solution. The stretch modes we are seeing even in the original adamantane system involve C–C stretches so shouldn’t be affected by a factor of the mass.
I’ve now had chance to look at the test problem and so here are my conclusions:
Short answer: Jmol and our brains struggle to cope with complex eigenvectors because we’re tied to real space (which is not unreasonable since we live there and Jmol is designed for molecules).
Long answer: The expression for the real space motion in terms of the complex eigenvectors is given according to expressions such as equation (18) in the file at:
What this shows is the real space part that we would like to visualise as motion depends not only on the eigenvectors (and mass), but also on a complex phase factor which causes the real space projection to oscillate with time. In order to visualise the motion it is therefore necessary to allow for this phase on top of the eigenvectors themselves & this is what is missing when we try to interpret the output (and Jmol is unaware of since it doesn’t know about the Brillouin zone).
If we take the specific case of the test where there is a periodic array of near-rigid diatomics in 1-D then there are 2 basis functions (the z coordinates of the atoms) and 2 vibrations along the chain. Therefore there will always have to be an in and out of phase mode (translation and stretch at gamma). The stretch is the mode at 2128 cm^-1 which shows little dispersion. The mode at 1.58i cm^-1 at 0.5 in this k direction is still actually an acoustic mode even though the real space part of the eigenvector looks like it would cause the bond to stretch because the values are different for the 2 atoms. However, if you evaluate the magnitude of the coefficient squared for each atom it is 0.5 and the real space direction is the same. It’s just that one atom’s motion is split between real and reciprocal space, and the other isn’t. This is just a difference due to the phase factor term because the atoms are at different relative positions in real space and so have different phase factors. If you transform the eigenvectors of the first atom according to the relative phase factor then it should rotate the components in the complex plane to match the other atom (that’s the theory at least).
Bottom line: It’s easier just to view eigenvectors at gamma.
There is now a modified version of GULP-6.1 available where you should be able to choose the type of eigenvectors you’d like in the output and associated files. The new option to add is:
Many thanks once again for your help with this – that has pinpointed our problem exactly.
Of course you are right that, given the correct eigenvalues, this had to be a problem with how the eigenvectors were being interpreted. In fact, though, Jmol does know about Brillouin zones, so it isn’t simply that we’re ignoring the phase factor entirely. The problem turns out to be the ambiguity in defining the displacements associated with each mode: in the \exp(i\mathbf{Q}\cdot\mathbf{r}_j) term, is \mathbf{r}_j the lattice vector associated with the relevant cell (as in the document you linked to above), or the position vector of the atom itself? It seems that GULP is assuming the first version and Jmol the second. Indeed, in my model example with atoms at z=\pm0.2, dividing through by \exp(2\pi i \times\frac12\times \pm0.2) by hand gives phases that Jmol interprets correctly.
This is precisely the “notorious different possible definitions of r” that in my original post I was so sure I had considered already – sorry!! I believe I was mistakenly thinking of the CASTEP to Jmol route which we indeed use very frequently. That is, I think CASTEP uses the “Jmol” rather than “GULP” definition.
So: if you are willing to add more flags to GULP, one that would be especially helpful to me would be one that converts between these two conventions. Of course if this is not easy or you don’t have time, it is not too tricky to divide by \exp(i\mathbf{Q}\cdot\mathbf{r}_j) (where \mathbf{r}_j is now the position vector of an atom in the original cell coordinates) by hand.
If I am right about CASTEP, it might also make sense to have this flag activated by default when outputting eigenvectors in CASTEP format? I will check and confirm my suspicion about which convention it uses.
Glad things are starting to make sense. I obviously underestimated JMol as well, since I assumed (as a non-user) that the “mol” might exclude the “sol” case.
On the rj aspect, I was assuming that this represented the atomic positions since determines the relative phase for different sites. I’ll certainly investigate whether I can correct by the phase factor as well as an option. As you say, it may well be that this is expected for the CASTEP format in which case I’ll make it the default. As soon as I’ve made progress on the phase issue I’ll post this.
For what it’s worth I have now confirmed that CASTEP indeed uses the full atomic position vector; see N2-phonon.phonon which comes from CASTEP itself. This is also implied by Keith Refson’s writeup here.
Out of interest, do you have a preferred viewer for visualising GULP (or other) vibrational modes?
Thanks for confirming the expectation from CASTEP.
I’ve now uploaded a new version of GULP-6.1 that should by default output the eigenvectors that make sense to Jmol (I hope). In other words, dividing by the square-root of the mass and the phase factor for the atom are the default. The old eigenvectors can be obtained using the option:
eigvec original phased
Let me know if this doesn’t correct things as expected.
Great! I can confirm that this now behaves as expected. For future reference, Jmol expects eigv phase (i.e., the new default) and eigv original (not the new default – if my problem was the only reason you changed this default it might be worth reverting?).
For what it’s worth, something slightly odd is still going on with how Jmol chooses the amplitude of motion to display away from the gamma point (degenerate modes have radically different scales), but this seems to be very much on their end; the eigenvectors GULP outputs all look sensible to me. I will post back here for the benefit of future readers if I sort that out, but I don’t think it is anything you need worry about from a GULP perspective. For my immediate use case the main thing I need is eigenvector files that will work with our own analysis scripts in the same way as those generated by CASTEP, and I think I now have exactly that!
Once again, thank you very much for your time and effort, and for a wonderful program!
Thanks for all your feedback as well - I appreciate it.
One quick thing just to check regarding the defaults (which I don’t mind changing it is makes more sense):
There are 2 separate options for eigv - original vs mass-unweighted and phased vs phased.
The old default was:
eigv original phased
and the new one is:
eigv mass unphased
In other words, both options are flipped together (which I think makes sense). I just want to confirm that the new default with both options changed is the one that is correct for Jmol (aside from any amplitude issues)?
Many thanks,
Julian
Hang on, are you sure the old default is as you described it? Jmol expects eigv original phased, and it looks to me as if the previous default was eigv original unphased. I don’t think my problem needed the mass change at all, though of course it may be useful to other users.
Hi Anthony,
Thanks for clarifying that the mass change isn’t needed for Jmol. As for the default, it’s probably just confusion over phased vs unphased since both are in some senses are phased. In the latest version both settings are flipped relative to the previous ones. I’ll have a think about what should be the new default, but I’ll certainly add a note as to what should be specified for Jmol to get correct output.
Best regards,
Julian
I realised over the weekend that I’d flipped the logic in the code for the phased vs unphased and so I’ve now corrected this in GULP. I’ve also changed the default to be what Jmol would expected “eigvec original unphased” (after the correction). Sorry for the mix up.
Best regards,
