# Set up configuration
molecule_configuration = MoleculeConfiguration(
    elements=[Nitrogen, Hydrogen, Hydrogen, Hydrogen],
    cartesian_coordinates=[[ 6.13508 ,  5.790587,  2.75    ],
                           [ 6.13508 ,  6.73176 ,  2.336663],
                           [ 6.95016 ,  5.32    ,  2.336663],
                           [ 5.32    ,  5.32    ,  2.336663]]*Angstrom
    )

# Add metallic regions (without them the external potential would be constant)
metallic_region_0 = BoxRegion(
    0*Volt,
    xmin = 0*Angstrom, xmax = 12*Angstrom,
    ymin = 0*Angstrom, ymax = 12*Angstrom,
    zmin = 0*Angstrom, zmax = 1*Angstrom
)

metallic_region_1 = BoxRegion(
    1*Volt,
    xmin = 0*Angstrom, xmax = 12*Angstrom,
    ymin = 0*Angstrom, ymax = 12*Angstrom,
    zmin = 5*Angstrom, zmax = 6*Angstrom
)

metallic_regions = [metallic_region_0, metallic_region_1]
molecule_configuration.setMetallicRegions(metallic_regions)

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

# Calculate and save the effective potential
external_potential = ExternalPotential(molecule_configuration)
nlsave('results.nc', molecule_configuration)
nlsave('results.nc', external_potential)