Some questions about vibration entropy calculation

I used bond stiffness-length model with fitsvsl -er=5 -ns=3 -ms=0.02 command to calculate the vibration entropy of zincblende ZnOS. The structures included in the calculation were structure 0-6. I found that the fitting of bond stiffness-length was changed with the involved structures. For example, if I used only structure 0-5 to fit bond and calculate phase diagram, I got the results as shown in Fig.1(a)-(c). However, if I included structue 6, I got vey different results, as shown in Fig.2, from the previous results.

So, I have some questions about the vibration calculation.How many structures are usually included to get convergence? And how to select structures makes the calculation efficiently?

The other question is about the discontinuous transition in the phase diagram. Is it naturally or due to the calculation error? Based on the previous results, I guess it due to the fitting error.

The stiffness vs length method has its limitations: there is spread in the data points and we fit a linear (or polynomial) model through them. If the data does not fit that model well, you can get variations as you add more data.
You could investigate this with:
gnuplot fitsvsl.gnu
I should add that, in your case, I have to warn you that for a ionic solid, there are long-range electrostatic interactions, so the nearest-neighbor view of stiffness vs length method may be a problem.
Also: have you tried the -op -pc or -dd options?
Perhaps fitfc would be a better approach? Also there are linear response codes that would perform better in ionic solids (but typically much more computationally expensive!).

BTW, I don’t see the figures (did you press "upload"?)…

The figures were included in the tar.gz file.
I’ve tried to include more structures and -op=3. The bond stiffness and the phase diagram seems convergency.

The figures for length vs stiffness look very good to me. There will always be some scatter around the trends.
Your phase diagrams really indicate some phase transitions. You should make sure that you are calculating the correct 2-phase equilibria. Start at 0K from two adjacent ground states.

I think I start from two adjacent ground states. I define lat.in as zincblende structure and put O and S in the same atoms position by using "O,S" form. Then the maps generated structure 0 for zb-ZnO and 1 for zb-ZnS. After eci was generated, phb -gs1=0 -gs2=1 -dT=25 -dx=1e-2 -er=50 -k=8.617e-5 -ltep=5e-3 -o=ph_wo_vibration.out was employed to obtain the phase diagram. All command were followed by user guide and tutorial.

The phase diagram without vibration free energy looks good. However, when the vibration free energy is included, the phase diagram shows so call phase transitions. I guess two possible reasons resulting in it. One is due to fitsvsl -ms=0.02. In the high temperature, the strain may be larger than the maximum strain 2%. The so call phase transition is due to the short of large strain data. The other possible reason is that the number of structure included into stiffness-bond fitting still lack. I will try to figure out the reason. What’s your opinion?

Another possibility is that the quasi harmonic approx is not so good at very high T. Try re-fitting using only vol_0 (this will revert to a plain harmonic approx, which is sometimes more reasonable).
Alternatively, if you want to use the quasiharmonic approx, specify the bulk modulus explicitly with the -b of -bf option of svsl (calculate it separately in a different directory. Maybe ZnO and ZnS and interpolate linearly). This may work better because svsl, by default, uses the nearest neighbor FC to calculate the bulk mod, which is not the best approximation, especially for an ionic solid!
More generally, look at the fvib for individual structure to see if they make sense.