from ATK.KohnSham import *

# Restore previously created silicon structure
si_file = VNLFile("si.vnl")
si_fcc  = si_file.readAtomicConfigurations()["si_fcc"]

# Define basis set and mesh cut-off
basis = basisSetParameters(type=SingleZeta)
mesh  = electronDensityParameters(mesh_cutoff=75.0*Rydberg)

print '\nConverging bulk fcc Si in number of k-points'
print 'kpts\t\tTotal energy (eV)\t(rel. diff)'
print '---------------------------------------------'

convergence_criterion = 1.0e-5

kpts = 0
E = Eprev = 1.0*eV

# Loop determine number of k-points for convergence
while ( (kpts<1) or (abs( (Eprev-E)/E ) > convergence_criterion) ):
    # Increase number of k-points for each iteration
    kpts = kpts + 1
    Eprev = E
    kgrid = brillouinZoneIntegrationParameters(
        monkhorst_pack_parameters=3*[kpts]
        )
    dft_method = KohnShamMethod(
        basis_set_parameters=basis,
        electron_density_parameters=mesh,
        brillouin_zone_integration_parameters=kgrid
        )
    # Perform actual calculation
    if kpts > 1:
        # Use the previous scf as initial calculation
        scf = executeSelfConsistentCalculation(
            atomic_configuration=si_fcc,
            method=dft_method,
            initial_calculation=scf,
            )
    else:
        # For (1,1,1) there is no initial calculation to use
        scf = executeSelfConsistentCalculation(
            atomic_configuration=si_fcc,
            method=dft_method
            )
    E = calculateTotalEnergy(scf)
    if kpts > 1:
        print 3*(kpts,), '\t', E.inUnitsOf(eV), '\t\t', abs((Eprev-E)/E)
    else:
        print 3*(kpts,), '\t', E.inUnitsOf(eV)
        
print '\nConverged with',kpts,'k-points in each direction\n'