# Set up configuration
molecule_configuration = MoleculeConfiguration(
    elements=[Nitrogen, Hydrogen, Hydrogen, Hydrogen],
    cartesian_coordinates=[[ 0.,       0.,        0.124001],
                           [ 0.,       0.941173, -0.289336],
                           [ 0.81508, -0.470587, -0.289336],
                           [-0.81508, -0.470587, -0.289336]] * Angstrom)

# Define the calculator
calculator = HuckelCalculator()
molecule_configuration.setCalculator(calculator)

# Calculate and save the mulliken population
mulliken_population = MullikenPopulation(
    configuration=molecule_configuration)

nlsave('mulliken.hdf5', mulliken_population)

# print all occupations
nlprint(mulliken_population)

# print partial occupations of N
print('N ', mulliken_population.atoms(spin=Spin.Sum)[0])
print(' | ', mulliken_population.orbitals(spin=Spin.Sum)[0])
print('')

# print partial occupations of first H
print('H ', mulliken_population.atoms(spin=Spin.Sum)[1])
print(' | ', mulliken_population.orbitals(spin=Spin.Sum)[1])
print('')

# print Mulliken population of N-H bond
print('N-H bond ', mulliken_population.bond(0, 1))