How to slove " you are out of the statistical sampling"

[Guodonglin-cqust]

Based on the QE, the structure optimization and harmonic dyn are performanced.

the structure has 2 atoms, Se and Sn.

the 1000 configurations are produced using sscha at 300K, but the trend of free energy change is actually rising, when it should be decreasing.

# Setup the minimization parameter for the covariance matrix
minimizer.min_step_struc = 0.01
minimizer.min_step_dyn = 0.01
minimizer.gradi_op = "all"
minimizer.precond_dyn = True
minimizer.root_representation = "normal"
minimizer.kong_liu_ratio = 0.5 # Usually 0.5 is a good value
minimizer.meaningful_factor = 0.00000001

 # ---------------- NEW MINIMIZATION STEP --------------------
Step ka =  45
Minimization step, force computed: 1000
Step too large (scalar = 4.9118081490355765 | kl_ratio = 0.8962609884821452), reducing to 0.00444444444653562


Number of symmetries before the step:  6
Harmonic contribution to free energy =    -142.96643902 meV
Anharmonic contribution to free energy =   -7996.42539978 +-       0.25355388 meV
Free energy =   -8139.39183881 +-       0.25355388 meV
FC gradient modulus =   30078.57411008 +-     178.49024039 bohr^2
Struct gradient modulus =       0.00000000 +-       0.00000000 meV/A
Kong-Liu effective sample size =  488.9700931584258

Check the stopping criteria: Running =  True

 # ---------------- NEW MINIMIZATION STEP --------------------
Step ka =  46
Minimization step, force computed: 1000
Good step found with 0.00444444444653562, try increment


Number of symmetries before the step:  6
Harmonic contribution to free energy =    -142.96643902 meV
Anharmonic contribution to free energy =   -7996.42539978 +-       0.25355388 meV
Free energy =   -8139.39183881 +-       0.25355388 meV
FC gradient modulus =   30097.73439323 +-     178.56357103 bohr^2
Struct gradient modulus =       0.00000000 +-       0.00000000 meV/A
Kong-Liu effective sample size =  488.9700931584258

KL: 488.9700931584258 KL/N: 0.4889700931584258 KL RAT: 0.5
  According to your input criteria
  you are out of the statistical sampling.
Check the stopping criteria: Running =  False
ROOT  NAME: minim_1
SAVE  NAME: minim_1
FNAME NAME: None

Restoring the last good dynamical matrix.
Updating the importance sampling...


 * * * * * * * * 
 *             * 
 *   RESULTS   * 
 *             * 
 * * * * * * * * 


Minimization ended after 46 steps

Free energy =   -8139.39183881 +-       0.25355388 meV
FC gradient modulus =   30097.73439323 +-     178.56357103 bohr^2
Struct gradient modulus =       0.00000000 +-       0.00000000 meV/A
Kong-Liu effective sample size =  488.9700931584258

Total force on the centroids [eV/A]:
   0)       0.000000      0.000000      0.000711  +-       0.000000      0.000000      0.001366
   1)       0.000000     -0.000000     -0.000711  +-      -0.000000     -0.000000      0.001366

[DjordjeDangic]

Hi,

Now you take the output dynamical matrix and generate new 1000 configurations, calculate forces, energies, and stresses for those, and rerun the minimization. You do this until the force constant gradient modulus becomes 0. I would advise to use less configurations in the first few iterations (I usually start with 100/200).

[Guodonglin-cqust]

I refer to the tutorial https://perturbo-code.github.io/mydoc_tdep.html, to obtain the temperature-renormalized phonon dynamics matrix to further calculate the electrical transport properties. QE was used to optimize the structure, ph, x and sscha were used to generate the second-order kinetic matrix, and the temperature-dependent dyn was generated according to the tutorial. However, it was found that the free energy increased during the minni process and did not decrease at all. This problem has been troubling for 3 months, still can 't solve this problem, can you give some advice, thank you.

[mesonepigreco]

Hi, as @DjordjeDangic mentioned, you are still very far from the free energy minimum.
For how the SSCHA is defined, the gradient is the exact gradient of the free energy only in the limit for an infinite number of configurations, i.e., for any finite ensemble, the free energy and the gradient could not match.

The gradient, however, converge much faster than the free energy, therefore, it is possible that your free energy seems to increase during the minimization as an artifact due to the number of configurations that is too low.

You can ignore it and try to proceed with the following population until convergence.
If, however, the free energy continues to increase across different ensembles, this is then a sign of something wrong.
The most likely source of error is a mismatch between the units of energy and forces given to the SSCHA.
Are you using the automatic cluster submission? Are you using ASE to get energy forces automatically? Are you using a custom interface with the code?

You can refer to the official SSCHA tutorial for more details: www.sscha.eu. In particular, the first one explains all the procedures of the sscha minimization. The out-of-statistical sample is not an error; it is just a message telling you that a new ensemble is required (see the tutorial).

[Guodonglin-cqust]

Regarding the SSCAH code, I replaced QE with VASP for calculating forces and energies. I even modified Structure.py to generate QE structure files (e.g., scf_population1_1.dat) sorted by atom types. However, I ultimately observed that the free energy kept increasing with simulation steps, and realized I overlooked a critical issue: the atomic displacements were not aligned with the corresponding atom ordering. To resolve this, I reverted to the original Structure.py file (where QE outputs unsorted structures) and wrote a script to reorder VASP forces and energies according to QE's structure format. This successfully reduced free energy with steps. However, this approach requires manual script modifications for different material compositions and ratios, lacking universality. I want to ask: Can SSCAH also output displacements sorted by atom types (e.g., in u_population1_1.dat files) to ensure consistency with the sorted scf_population1_1.dat files? This would perfectly solve the SSCAH+VASP compatibility issue.
For the displacement files in SSCAH, which are handled in Ensemble.py, I’m unsure how to modify the code to ensure the displacements are ordered by atom type. Can you provide specific guidance on adjusting Ensemble.py to output sorted displacement files (e.g., u_population1_1.dat) that align with the atom-type-sorted structure files like scf_population1_1.dat?

[mesonepigreco]

Thanks. Maybe it would be better to implement a different saving function for VASP, like save_vasp and load_vasp, which you could then use to reorder the atoms by type and order back correctly.
This will make the two saving mechanisms incompatible, which is good, as people will not mistake one ensemble with the other (missing the atom ordering).

[Guodonglin-cqust]

This is a good idea, but I'm not familiar with ASE, so I probably can't write a saving function with universal applicability.