Outline:
1) Get vasprun.xml from a VASP phonon calculation.
2) Use Phonopy to find the irreducible representations of each phonon mode.
3) Separately find the Raman/IR selection rules for active modes using tabulated space group information (e.g. using the Bilbao Crystallographic Server).
4) Compare 2) and 3) above to find which of your modes are Raman/IR-active!
Quantum ESPRESSO is quite masterful/missterful/themterful/zemterful/pronounterful at phonon-related mumbo jumbo. The ph.x routine is really useful in finding phonon modes and their symmetries, and calculating their dispersions, among many other things. However, sometimes, a problem might necessitate the usage of VASP. However, when it comes to phonons, VASP is the equivalent of Greg. Bruh, you can’t even calculate Raman intensities without using some long-winded approach like raman-sc. Although I did not try hard enough to figure out how to use raman-sc on many nodes – which, unfortunately, were needed for my problem – I am VERY grateful it exists. Otherwise we would have had the final form of Greg.
Anyway, I had to use VASP for a specific problem involving spin-orbit coupling (SOC). Adding SOC to my setup did not make the Quantum ESPRESSO version converge despite trying several different things. AND I had to provide some experimental collaborators with something to supplement their Raman measurements by figuring out which of the DFT-calculated phonon modes are Raman-active. So, here’s how I went about doing it (the same steps apply for finding IR-active modes):
First, run a phonon calculation. I like the DFTP method outlined in Danny Rehn’s tutorial (which is also the source of the following code). Just make sure you obtain the calculation’s vasprun.xml file.
Then, use Phonopy to get the modes’ irreducible representations. I know, I know, Phonopy has a reputation for doing other things, but it works here too. To do this, first generate a FORCE_CONSTANTS file using your vasprun.xml:
phonopy --fc vasprun.xmlThen create a file symm.conf with the following code, to compute the irreducible representations (irreps) at the Gamma point (0, 0, 0). The tolerance is set to 0.001, and it is used to identify degenerate frequencies. See the Phonopy manual here for more information.
IRREPS = 0 0 0 1e-3
SHOW_IRREPS = .TRUE.Then, run the following to get the irreps:
phonopy --readfc --dim="1 1 1" symm.confOpen the file irreps.yaml. Each frequency will have an “ir_label” (e.g., A1g, Eg, or even “None”). I believe the “ir” in “ir_label” stands for “irreducible representation” (and not “infrared”). We use these labels to compare against selection rules.
To find the correct selection rules, head over to the Bilbao Crystallographic Server. Click “Raman and Hyper-Raman scattering“, and then “SAM” (Spectral Active Modes). Here’s the link.
Enter your crystal’s space group or upload the crystal structure. For space group information, you can use existing knowledge of the material (maybe from the Materials Project (see “International Number“), or just any of the other examples I discuss in this post, including Phonopy).
If asked for the Wyckoff Positions, you can use what you know or just find them on Materials Project page for the material. Or you could get this information from the Phonopy output when “–symmetry” is run (as in this post): Some atoms will have their Wyckoff positions listed in the output (e.g., ‘i’, ‘d’). Simply make a list of all of these for the input structure, and check the corresponding boxes in the Bilbao SAM page. For example, for a monolayer Nb₃Cl₈ structure with space group 164, Phonopy gave me Wyckoff positions ‘i‘ for atoms 1, 7 and 17, and ‘d‘ for atoms 13 and 15. When I entered 164 in the Bilbao SAM website, my Wyckoff position options were: 12j, 6i, 6h, 6g, 3f, 3e, 2d, 2c, 1b, 1a. So, I checked the boxes for 6i and 2d. Recall that each letter corresponds to a unique arrangement of atoms or ions with specific symmetry properties, and the numbers represent the multiplicity of a particular Wyckoff position. Only about the letters are needed in this case.
You will then see selection rules for Raman and IR-active modes. Something like the screenshot below:
The table cells with a number (i.e., not a ‘dot’) are active modes. So, just find which of the modes in irreps.yaml match these using ir_label. Aaaaaaand… you are done!
Very helpful. This has greatly accelerated my progress. Thank you for writing this out and sharing!!
I am still struggling to get the right character tables for all of my modes. For most phonopy prints:
“Not found. Try adjusting tolerance value in IRREPS.”
I tried this with a few different values in symm.conf but it did not work so far.
You are most welcome! Would you mind looking at my response to Sanger He’s comment below? Sorry there’s no easy answer.
Hi, just a quick update. What helped for my system was turning the EDIFF dial way down. At first I tried with EDIFF=1e-7 which seemed tight, then I was advised by a colleague to go lower. EDIFF=1e-9 was better (more modes symmetries identified), and EDIFF=1e-11 did the job. I thought that I would get overconvergence or numerical errors but much to my surprise, it worked nicely.
Now I am trying to figure out a supercell calculation (1x1x2). VASP reports twice as many modes, but from phonopy I am getting the same amount of modes and irreps as for the primitive cell. I tried adjusting the second command to:
phonopy –readfc –dim=”1 1 2″ symm.conf
but it does not do the job.
That about EDIFF is very interesting. I guess the more-stringent convergence criterion makes the forces between atoms converge more accurately? But I can totally see it happening the other way around as well. An EDIFF convergence test might be a good safety measure.
As for the supercell calculation, I could be wrong, but I think phonopy calculates modes based on the primitive cell regardless. Do consult the documentation! Also, I’m not sure what specifically you were referring to when you wrote “VASP reports twice as many modes…”.
Good luck!
Here is the clarification. I am using phonopy to post-process a VASP finite differences calculation of a supercell (1x1x2). If n is the number of modes for the primitive cell, VASP now reports 2*n frequencies. When I post process with phonopy, it only works after I add the following line to symm.conf:
PRIMITIVE_AXES = AUTO
Now it reports again n frequencies. Each ir label is the same as for the primitive cell. Perhaps all of the additional modes “fold” into the primitive cell but I am not sure that this is generally true. I can imagine some modes that cannot be reduced to the primitive cell. I am thinking eg. about flexural modes in a free-standing graphene monolayer.
I am writing this mainly for other people who stumble on this article and have similar problems. Don’t worry if you do not have all the solutions. Thank you once again for this useful post and replies!!
In my irreps.yaml file there are too many “none” at the ir_reps symbol, in the output it says “Not found. Try adjusting tolerance value in IRREPS.”, but this won’t work, even when decreasing IRREPS the space group and the point group will degenerate from P3M1(156) to Cm(8), how to solve this problem?
Hey there! It is a bit hard for me to say much without more specifics… this can be very problem-dependent and require a lot of custom troubleshooting. As far as I know, there is nothing wrong with “none” and “not found”. As this could be very dependent on the symmetry you start with. So I would recommend you run your DFTP calculation for an input structure that has the expected space group (compare with the literature), after adjusting symmetry tolerances as needed (see this post: http://thauwa.me/2022/09/15/how-to-ensure-vasp-will-detect-expected-symmetry-in-poscar-input/). And make sure the DFTP calculation detects the expected symmetry by investigating OUTCAR (there might be mismatches between the initial and final compounded space group). These should match what phonopy detects. When I encountered this, I had to play around with symmetry threshold parameters before and after DFTP to get what I expected from the literature. Good luck!