# ------------------------------------------------------------- 
# Bulk Configuration 
# ------------------------------------------------------------- 

bulk_configuration = nlread('Pd100_CO.nc', BulkConfiguration)[-1] 

# Add tags 
bulk_configuration.addTags('CO',          [4, 5]) 
bulk_configuration.addTags('Pd',          [0, 1, 2, 3]) 
bulk_configuration.addTags('Selection 0', [0, 1]) 

# ------------------------------------------------------------- 
# Calculator 
# ------------------------------------------------------------- 

# Get the calculator settings from non-BSSE calculation. 
calculator = bulk_configuration.calculator() 
basis_set = calculator.basisSet() 
exchange_correlation = calculator.exchangeCorrelation() 
numerical_accuracy_parameters = calculator.numericalAccuracyParameters() 
poisson_solver = calculator.poissonSolver() 

# Create BSSE calculator. 
bsse_calculator = counterpoiseCorrected(LCAOCalculator, ["CO", "Pd"]) 
calculator = bsse_calculator( 
    basis_set=basis_set, 
    exchange_correlation=exchange_correlation, 
    numerical_accuracy_parameters=numerical_accuracy_parameters, 
    poisson_solver=poisson_solver, 
    ) 

bulk_configuration.setCalculator(calculator) 
nlprint(bulk_configuration) 
bulk_configuration.update() 
nlsave('Pd100_CO_cp.nc', bulk_configuration) 

# ------------------------------------------------------------- 
# Optimize Geometry 
# ------------------------------------------------------------- 
constraints = [0, 1] 

bulk_configuration = OptimizeGeometry( 
        bulk_configuration, 
        max_forces=0.05*eV/Ang, 
        max_steps=200, 
        max_step_length=0.2*Ang, 
        constraints=constraints, 
        trajectory_filename='Pd100_CO_cp.nc', 
        disable_stress=True, 
        optimizer_method=LBFGS(), 
        ) 
nlsave('Pd100_CO_cp.nc', bulk_configuration) 
nlprint(bulk_configuration) 

# ------------------------------------------------------------- 
# Total Energy 
# ------------------------------------------------------------- 
total_energy = TotalEnergy(bulk_configuration) 
nlsave('Pd100_CO_cp.nc', total_energy) 
nlprint(total_energy)