# Set minimal log verbosity
setVerbosity(MinimalLog)

# Define elements
elements = [Carbon, Nitrogen, Hydrogen, Hydrogen, Carbon, Hydrogen, Hydrogen,
            Oxygen, Oxygen, Hydrogen]

# Define coordinates
cartesian_coordinates = [[ 0.004212789037,  0.000971957487,  0.047351665502],
                         [ 1.466200629874, -0.05214671952 ,  0.084234583105],
                         [-0.342946287509, -0.002631605963, -1.00793020604 ],
                         [-0.413982078001, -0.900683251432,  0.54677392013 ],
                         [-0.503977231839,  1.23896512228 ,  0.720474876604],
                         [ 1.863040371746,  0.721978981544, -0.497567022566],
                         [ 1.806267691758,  0.033734459151,  1.070750195227],
                         [-0.210544217266,  2.422528604664,  0.232223616555],
                         [-1.293023389163,  1.146950982333,  1.766460567445],
                         [ 1.786364781973, -0.967997544378, -0.30807170096 ]]*Angstrom

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

# Add bonds
bonds = [[0, 1],
         [0, 2],
         [0, 3],
         [0, 4],
         [1, 5],
         [1, 6],
         [1, 9],
         [4, 7],
         [4, 8]]
configuration.setBonds(bonds)


# Create LCAOCalculator
solvation_parameters = CosmoSolvationParameters(
    solvent_dielectric_constant=78.3553,
    solvent_surface_tension=72.8*dyne/cm
)
calculator = LCAOCalculator(
    solvation_parameters=solvation_parameters
)

# Set Calculator
configuration.setCalculator(calculator)
configuration.update()
nlsave('glycine_zwitter_ion.hdf5', configuration)

# Optimize the geometry
optimized_configuration = OptimizeGeometry(
    configuration=configuration,
    trajectory_filename='glycine_zwitter_ion.hdf5'
)
nlsave('glycine_zwitter_ion.hdf5', optimized_configuration)

# Get the solvent surface
solvent_surface = SolventSurface(
    configuration=optimized_configuration
)
nlsave('glycine_zwitter_ion.hdf5', solvent_surface)