# -------------------------------------------------------------
# Molecule Configuration
# -------------------------------------------------------------

# Define elements
elements = [Oxygen, Carbon]

# Define coordinates
cartesian_coordinates = [[-1.2987,  0.0585,  0.    ],
                         [ 0.    ,  0.    ,  0.    ]]*Angstrom

# Set up configuration
molecule_configuration = MoleculeConfiguration(
    elements=elements,
    cartesian_coordinates=cartesian_coordinates
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Carbon_DoubleZetaPolarized,
    GGABasis.Oxygen_DoubleZetaPolarized,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = GGA.PBE

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

molecule_configuration.setCalculator(calculator)
nlprint(molecule_configuration)
molecule_configuration.update()
nlsave('CO.nc', molecule_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
molecule_configuration = OptimizeGeometry(
        molecule_configuration,
        max_forces=0.05*eV/Ang,
        max_steps=200,
        max_step_length=0.2*Ang,
        trajectory_filename='CO.nc',
        disable_stress=True,
        optimizer_method=LBFGS(),
        )
nlsave('CO.nc', molecule_configuration)
nlprint(molecule_configuration)

# -------------------------------------------------------------
# Total Energy
# -------------------------------------------------------------
total_energy = TotalEnergy(molecule_configuration)
nlsave('CO.nc', total_energy)
nlprint(total_energy)