Strange convergence of higher order force constants

[utf]

Summary

I am trying to converge my higher order force constants for use in SCPH calculations. However, when I increase the number of training structures the SPCH free energy does not converge smoothly.

Any suggestions would be greatly appreciated.

Generating training set data

To generate the training set, I did:

1. A long MD calculation with fairly coarse settings, gamma point only, larger energy convergence. The MD calculation ran for 10,000 steps.
2. Selected samples from MD every 50 steps to give 200 structures total, displaced the atoms by 0.01 Å, and calculated those structures using very tight settings (same as used for harmonic calculations which gave accurate results). The command I used was: displace.py --VASP=POSCAR.super -md ../md/vasprun.xml --random --mag=0.01 .

Note my system is FCC, with 10 atoms in the unit cell. I did force constant calculations in a 2x2x2 supercell which I found converged the harmonic phonon frequencies. I only did small random displacements on the MD as the smallest interatomic distance is already quite small (1.16 Å) but maybe this should be increased?

Converging SPCH free energy

To test convergence, I fit the force constants using a certain number of structures and calculate thermodynamic properties using SCPH (off diagonal elements and bubble correction are not included).

I expected due to the high symmetry and large number of samples, that thermodynamic properties would converge quickly but this is not the case. Below see the plot showing the convergence of QHA term in the free energy with supercell size (note I fitted force constants up to 6th order with relatively large cutoffs, I also tried smaller cutoffs but exactly the same trend was observed. I also tried only fitting up to 4th order but this made things worse). Finally, note that the validation error for fitting 6th order force constants with 200 samples is quite small (~1.5 %).

Convergence at 0 K

Convergence at 1000 K

Do you have any suggestions for what is going on? Is my MD calculation too coarse, is there an instability in the SCPH calculation, or is there something else that is important for convergence?

Input files

MD

Incase they are useful, here is my INCAR for the MD calculation:

 ENCUT = 500.000000
 SIGMA = 0.050000
 EDIFF = 1.00e-05
 GGA = PS
 ALGO = Normal
 ISYM = 0
 ISMEAR = 0
 NCORE = 20
 NELMIN = 5
 NWRITE = 2
 NELM = 200
 LCHARG = .FALSE.
 LWAVE = .FALSE.
 LASPH = .TRUE.

# molecular dynamics
 MDALGO = 2 #NVT
 NSW    = 10000
 IBRION = 0
 ISIF   = 2 #NVT
 SMASS  = 1
 POTIM  = 2.0
 TEBEG  = 500

Displacements

 ENCUT = 700.000000
 SIGMA = 0.050000
 EDIFF = 1.00e-08
 SYMPREC = 1.00e-03
 EDIFFG = -1.00e-04
 PREC = Accurate
 GGA = PS
 ALGO = Normal
 ISYM = 0
 ISMEAR = 0
 NCORE = 20
 NELMIN = 5
 NWRITE = 2
 NELM = 200
 NSW = 0
 IBRION = 1
 ISIF = 3
 LCHARG = .FALSE.
 LWAVE = .FALSE.
 LASPH = .TRUE.
 ADDGRID = .TRUE.
 LREAL = .FALSE.

Fitting force constants

...

&interaction
  NORDER = 5  # 1: harmonic, 2: cubic, ..
/

&cutoff
  *-* None 10 7 7 7
/

&optimize
  DFSET = DFSET_cubic
  FC2XML = ../zn222.xml # Fix harmonic IFCs
  LMODEL = enet
  L1_RATIO = 1
  L1_ALPHA = 1e-06
  CV = 0
  CV_MINALPHA = 1e-6
  CV_MAXALPHA = 0.05
  CV_NALPHA = 10
  STANDARDIZE = 1
/

Calculating SCPH

...

&kpoint
  2  # KPMODE = 2: uniform mesh mode
  30 30 30
/

&scph
  SELF_OFFDIAG = 0
  KMESH_INTERPOLATE = 2 2 2
  KMESH_SCPH = 4 4 4
/

&analysis
  FE_BUBBLE = 0
/

[ttadano]

When fitting the force constants to the data by using elastic net (LMODEL = enet), an optimal value of L1_ALPHA changes with the size of the training dataset. Usually, it decreases with increasing NDATA. Did you rerun cross-validation (CV > 0) when you increased the training data? If no, please do so and find an optimal value of L1_ALPHA for each NDATA.

If the above treatment does not help, please try using LMODEL = adaptive-lasso. Adaptive LASSO tends to give a more sparse solution than LASSO (enet). I observed a faster convergence of the free energy with adaptive LASSO.