Regarding the convergence of SSCHA calculations

[hellolori]

Dear developers,

I have some lingering questions after going through the topic titled "How can I assess the convergence of my calculations?" (found at this link. The topic suggests that, in general, if the gradient approaches zero and the Kong Liu ratio exceeds 0.5, then my calculations have likely converged satisfactorily. I've conducted calculations for H3S three times using identical input files as shown below in Python:

# Import the sscha code
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities
# Import the cellconstructor library to manage phonons
import cellconstructor as CC, cellconstructor.Phonons
import cellconstructor.Structure, cellconstructor.calculators
# Import the DFT calculator
import cellconstructor.calculators
from ase.calculators.espresso import Espresso
# Import numerical and general pourpouse libraries
import numpy as np, matplotlib.pyplot as plt
import sys, os
#
# Initialize the DFT (Quantum Espresso) calculator for H3S
# The input data is a dictionary that encodes the pw.x input file namelist
input_data = {
    'control' : {
        # Avoid writing wavefunctions on the disk
        'disk_io' : 'None',
        # Where to find the pseudopotential
        'pseudo_dir' : '.',
        'outdir' : './outdir',
        'wfcdir' : './wfcdir',
        'tprnfor'  : True,
        'tstress' : True
    },
    'system' : {
        # Specify the basis set cutoffs
        'ecutwfc' : 35,   # Cutoff for wavefunction
        'ecutrho' : 350, # Cutoff for the density
        # Information about smearing (it is a metal)
        'input_dft' : 'blyp', # exchange-correlation functional
        'occupations' : 'smearing',
        'smearing' : 'mv',
        'degauss' : 0.02
    },
    'electrons' : {
        'conv_thr' : 1e-8
    }
}
#
# the pseudopotential for each chemical element
# In this case just Gold
pseudopotentials = {'H' : 'H.pbe-rrkjus_psl.1.0.0.UPF', 'S' : 's_pbe_v1.4.uspp.F.UPF'}
#
# the kpoints mesh and the offset
kpts = (8,8,8)
koffset = (1,1,1)
#
# Specify the command to call quantum espresso
command = '~/software/openmpi/4.0.5/bin/mpiexec -n 4 ~/software/QE/6.6_openmpi/bin/pw.x  -npool 2 -i PREFIX.pwi > PREFIX.pwo'
#
# Prepare the quantum espresso calculator
calculator = CC.calculators.Espresso(input_data,
                                     pseudopotentials,
                                     command = command,
                                     kpts = kpts,
                                     koffset = koffset)
TEMPERATURE = 300
N_CONFIGS = 128
MAX_ITERATIONS = 20
START_DYN = 'start_sscha'
NQIRR = 3
# Let us load the starting dynamical matrix
dyn = CC.Phonons.Phonons(START_DYN, NQIRR)
dyn.Symmetrize()
#
# Initialize the random ionic ensemble
ensemble = sscha.Ensemble.Ensemble(dyn, TEMPERATURE)
#
# Initialize the free energy minimizer
minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minim.set_minimization_step(0.01)
minim.meaningful_factor=1e-5
#
# Initialize the NVT simulation
relax = sscha.Relax.SSCHA(minim, calculator, N_configs = N_CONFIGS,
                          max_pop = MAX_ITERATIONS)
#
ioinfo = sscha.Utilities.IOInfo()
ioinfo.SetupSaving("minim_info")
relax.setup_custom_functions(custom_function_post = ioinfo.CFP_SaveAll)
# Run the NVT simulation (save the stress to compute the pressure)
relax.relax(get_stress = True)
#
# Now we can save the final dynamical matrix
# And print in stdout the info about the minimization
relax.minim.finalize()
relax.minim.dyn.save_qe("sscha_T{}_dyn".format(TEMPERATURE))

Here's the data obtained during these calculations:
(1) The 1st run

(2) The 2nd run

(3) The 3rd run

From the figures above, it's evident that the gradient approaches zero and the Kong Liu ratio exceeds 0.5 at the conclusion of all three calculations, indicating convergence in line with the aforementioned topic. However, the resulting phonon frequencies from these three calculations exhibit significant differences, particularly in the high-frequency range, as illustrated below:

I have the following inquiries:
(1) How can I determine which result among the three calculations is the most reliable?
(2) Is it more reasonable to average the outcomes of the three calculations, including phonon frequencies and Hessian matrix elements?
(3) It appears that the time spent on a calculation is positively correlated with the number of configurations. How should I select a suitable value for the number of configurations within a population? Are there any empirical methods for providing a reliable estimate before initiating calculations?
(4) How can we determine whether the parameter "minim.meaningful_factor" is suitable for this calculation? It seems challenging to gauge this based solely on the convergence of phonon frequencies. This is because phonon frequencies can vary significantly, even when using identical input, as demonstrated above.

Best regards!

[diegomartinez2]

1. The final effective sample size of the first calculation is bigger than the other two, this implies that is the "most reliable" calculation.
   On the other hand 128 configurations is a low parameter for a reliable calculation in general. In the 3th. tutorial you have this figure (this is valid for the example system, but in general it follows a similar shape):
   

[diegomartinez2]

In order to check the reliability of the output you can redo with a bigger number of configurations and check the evolution of the results.

[diegomartinez2]

2. Yes, you can average the outcomes, but if the number of configurations is OK the difference should be small.
3. With the standard random method, the configurations are that, random. Again, with a large number of random configurations you can "map" the BO energy better. There is also the alternative of using the Sobol sequence instead of the random generator method that should need a smaller number of configurations to get a good sampling.