I have a little problem concerning hexagonal volume fitting using VASP. I have followed these steps:

  1. Relax the structure from a given volume using ISIF =4

  2. Copy CONTCAR to POSCAR and relax it again

  3. Run with ISMEAR = -5 and without relaxation

  4. Repeat for different volumes (+10% +5% 0% -5% -10%)

  5. Use your script for EOS fitting.

It worked successfully with cubic systems. But for a hexagonal system (e.g. Carbon), I got the figure below. I don't know why I always get wrong results for the first volume point. Even with 7 volume points the first point always wrong.Equation of state fitting

  • 1
    $\begingroup$ What point in this figure is the first point? When you say its "wrong", what is your desired margin of error? $\endgroup$
    – Tyberius
    May 23, 2020 at 19:24
  • 1
    $\begingroup$ +1. Seems like a good question. I agree with Tyberius that the problem could be a little bit better described. As someone who doesn't work as much in this field, I'm not sure what the "correct" behavior should look like. Should it be monotonic (meaning the 3 points with smallest volume, are "wrong" here), or should it be somewhat flat (meaning that the point with largest volume here, is wrong)? $\endgroup$ May 23, 2020 at 20:04
  • 2
    $\begingroup$ I think the problem is that the variations are to high. Normally the cell variations are around 5-7%. Higher than that can be considered as deformations. (I just edited the firs comment). $\endgroup$
    – Camps
    May 23, 2020 at 21:21
  • 2
    $\begingroup$ How are you varying the volume? And by hexagonal carbon you mean graphite? Does the first point also have a low energy when you increase the volume slightly (e.g. make the first point 36 Å$^3$)? Finally, does still occur with a different hexagonal crystal (e.g. metallic Ti, $\beta$-quartz)? $\endgroup$ May 23, 2020 at 23:58
  • 2
    $\begingroup$ The expected behavior should be parabolic for energy vs volume, which is what the OP is trying to say. Although I don't know the goal of this. Anyway, I would use ISIF = 2 to obtain energies. Another thing you could try is change a/b lattice constant sequentially instead of the volume if you can - i have had better success with it. If nothing works, maybe check if any of pymatgen's functions would help. $\endgroup$
    – gogo
    May 24, 2020 at 1:58

3 Answers 3


Normally this should follow a nice equation of state (Rose, Murnigham etc.). I'd say the curve looks okay: $35 \mathrm{A}^3$ seems to be the questionable outlier.

There is not enough information on the setup here, but I will give two thoughts anyway:

  1. Has the k-point convergence been tested?
  2. HCP lattices require gamma centered k-point grids: https://cms.mpi.univie.ac.at/vasp/vasp/hexagonal_lattices.html

There is little information in your question. It would be useful if you add your INCAR files. It seems to me that the problem could be in the fact that for such a cell there could be a “a free” parameter (c/a). So, you might be changing the volumen and the c/a ratio at the same time. Also, check that your are removing the density and wavefuncton files as a change in volumen changes the basis set


First of all, thanks to everyone who contributed to answer and discuss this question, your responses helped me to found the solution. Which is the vdW correction was missing in my calculations. Since that WS2 is a layered material, vdW interlayer correction must be added. So, I have added IVDW=10 in INCAR file to include the DFT-D2 method and this gaves me good results.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .