Somthing wrong with my test mesh-cut-off.

[fangyongxinxi]

Hello,

I run a script with a for-loop to test the mesh-cut-off with ATK12.8.r1. I can not find the right number, I wonder why. Is somthing wrong with my script , or my system ?

Thanks~

the input sript is below, and the total energy are:
cut-off  Etotal
30  -4647.2374
50  -4646.413
70 -4646.147
90 -4646.059
110 -4645.999
150 -4645.953
170 -4645.942
190 -4645.932
210 -4645.926
240 -4645.920
260 -4645.915
280 -4645.912

# input script
vector_a = [3.94206054094, -1.97163500648e-13, -1.97712127445e-14]*Angstrom
vector_b = [-1.97163500648e-13, 3.94206054094, -2.86317462121e-14]*Angstrom
vector_c = [-7.13930706432e-14, -1.0304864771e-13, 14.0108929811]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

elements = [Potassium, Potassium, Iron, Iron, Iron, Iron, Selenium, Selenium,
            Selenium, Selenium]

fractional_coordinates = [[  2.80143445e-14,   4.92916417e-14,  -8.54315314e-07],
                          [  5.00000000e-01,   5.00000000e-01,   5.00000506e-01],
                          [  9.48261082e-14,   5.00000000e-01,   2.50011259e-01],
                          [  5.00000000e-01,   2.12883543e-13,   2.50011259e-01],
                          [  1.70509499e-13,   5.00000000e-01,   7.49988543e-01],
                          [  5.00000000e-01,   2.05720897e-13,   7.49988543e-01],
                          [  2.22933031e-14,   3.41568038e-14,   3.69467062e-01],
                          [  5.00000000e-01,   5.00000000e-01,   8.69445676e-01],
                          [  5.26577799e-14,   1.87595793e-14,   6.30534560e-01],
                          [  5.00000000e-01,   5.00000000e-01,   1.30552332e-01]]


bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )
basis_set = [
    GGABasis.Potassium_SingleZetaPolarized,
    GGABasis.Iron_SingleZetaPolarized,
    GGABasis.Selenium_SingleZetaPolarized,
    ]
exchange_correlation = SGGA.PBE

# for loop for mcut
for mcut in (30,50,70,90,110,150,170,190,210,240,260,280):
    numerical_accuracy_parameters = NumericalAccuracyParameters(
        grid_mesh_cutoff=mcut*Hartree,
        k_point_sampling=(12, 12, 4),
        )

    calculator = LCAOCalculator(
        basis_set=basis_set,
        exchange_correlation=exchange_correlation,
        numerical_accuracy_parameters=numerical_accuracy_parameters,
        )

    bulk_configuration.setCalculator(calculator)
    bulk_configuration.update()
    nlsave('testpara.nc', bulk_configuration)

    # -------------------------------------------------------------
    # Total energy
    id = str(mcut)
    total_energy = TotalEnergy(bulk_configuration)
    nlsave('testpara.nc', total_energy,id)

[zh]

If you plot the data of total energy vs mesh cutoff, it may be more clear to see the convergence.

For the compound in your script file, it is better to use the DZP basis set.

[nori]

Your calculation seems right.
It is often difficult to get total energy numerically converged for mesh cutoff if GGA is used.
According to my experience, 500 to 1500 Ry is needed.

If you switch XC functional to LDA, you will get a good convergence with small  mesh cutoff

[fangyongxinxi]

Thanks for your answer, it helps me.