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

# Define coordinates
cartesian_coordinates = [[ 1.034527,  0.034736, -0.847537],
                         [ 0.060835,  0.006745, -0.350482],
                         [-0.520065,  0.892764, -0.635075],
                         [-0.482833, -0.89338 , -0.662936],
                         [ 0.364387, -0.008234,  1.010835],
                         [-0.456852, -0.032631,  1.485195]]*Angstrom

# Set up configuration
molecule_configuration = MoleculeConfiguration(
    elements=elements,
    cartesian_coordinates=cartesian_coordinates
)
# Add a tag to distinguish the hydroxyl-H from the methyl-H.
molecule_configuration.addTags('HO', [5])

# Set up bonds between the atoms.
molecule_configuration.findBonds()

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

# Set up a new TremoloXPotentialSet
potential_set = TremoloXPotentialSet(name='Methanol_bonded')

# Add particle types.
potential_set.addParticleType(
    ParticleType.fromElement(Carbon, charge=0.1166*elementary_charge)
)

potential_set.addParticleType(
    ParticleType.fromElement(Oxygen, charge=-0.6497*elementary_charge)
)

# Add particle types for methyl hydrogen (without tag).
potential_set.addParticleType(
    ParticleType.fromElement(Hydrogen, charge=0.0372*elementary_charge)
)

# Add particle types for hydroxyl hydrogen (with HO-tag).
potential_set.addParticleType(
    ParticleType.fromElement(Hydrogen, charge=0.4215*elementary_charge, tags='HO')
)

# Set up a new bond potential for C-C bonds and add it to the potential set.
bond_potential = HarmonicBondPotential(
    particleType1=ParticleIdentifier('C'),
    particleType2=ParticleIdentifier('H'),
    k=14.743795*eV/Ang**2,
    r_0=1.090*Ang,
)
potential_set.addPotential(bond_potential)

# Set up a new bond potential for C-H bonds and add it to the potential set.
bond_potential = HarmonicBondPotential(
    particleType1=ParticleIdentifier('C'),
    particleType2=ParticleIdentifier('O'),
    k=13.87651*eV/Ang**2,
    r_0=1.410*Ang,
)
potential_set.addPotential(bond_potential)

# Set up a new bond potential for C-H bonds and add it to the potential set.
bond_potential = HarmonicBondPotential(
    particleType1=ParticleIdentifier('H', tag='HO'),
    particleType2=ParticleIdentifier('O'),
    k=23.9803*eV/Ang**2,
    r_0=0.960*Ang,
)
potential_set.addPotential(bond_potential)

# Set up a new angle potential for H-C-C angles and add it to the potential set.
angle_potential = HarmonicAnglePotential(
    particleType1=ParticleIdentifier('H'),
    particleType2=ParticleIdentifier('C'),
    particleType3=ParticleIdentifier('O'),
    k=2.1682*eV,
    theta0=1.9111355,
)
potential_set.addPotential(angle_potential)

# Set up a new angle potential for H-C-H angles and add it to the potential set.
angle_potential = HarmonicAnglePotential(
    particleType1=ParticleIdentifier('H'),
    particleType2=ParticleIdentifier('C'),
    particleType3=ParticleIdentifier('H'),
    k=1.51774*eV,
    theta0=1.9111355,
)
potential_set.addPotential(angle_potential)

# Set up a new angle potential for H-C-H angles and add it to the potential set.
angle_potential = HarmonicAnglePotential(
    particleType1=ParticleIdentifier('C'),
    particleType2=ParticleIdentifier('O'),
    particleType3=ParticleIdentifier('H', tag='HO'),
    k=2.385*eV,
    theta0=1.893682,
)
potential_set.addPotential(angle_potential)
# Set up a new torsion potential for H-C-C-H groups and add it to the potential set.
# This potential has a multiplicity of 1.
torsion_potential = CosineTorsionPotential(
    particleType1=ParticleIdentifier('H'),
    particleType2=ParticleIdentifier('C'),
    particleType3=ParticleIdentifier('O'),
    particleType4=ParticleIdentifier('H', tag='HO'),
    k=[0.007227*eV],
    n=[3],
    delta=[0.0]
)
potential_set.addPotential(torsion_potential)

# Add a coulomb solver to the potential.
potential_set.setCoulombSolver(CoulombDSF(r_cut=15.0*Ang))

# Create a new TremoloXCalculator with this potential.
calculator = TremoloXCalculator(parameters=potential_set)

molecule_configuration.setCalculator(calculator)