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

# Define coordinates
cartesian_coordinates = [[ 0.987524, -0.003838, -1.18508 ],
                         [-0.020839, -0.022441, -0.751621],
                         [-0.569312,  0.833015, -1.166929],
                         [-0.514525, -0.932892, -1.115519],
                         [ 0.020891,  0.022496,  0.751587],
                         [ 0.568699, -0.833372,  1.166882],
                         [ 0.515082,  0.932608,  1.115691],
                         [-0.987522,  0.004424,  1.184988]]*Angstrom

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

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

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

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

# Add particle types for Carbon and Hydrogen.
potential_set.addParticleType(
    ParticleType.fromElement(Carbon, charge=-0.3*elementary_charge)
)
potential_set.addParticleType(
    ParticleType.fromElement(Hydrogen, charge=0.1*elementary_charge)
)

# Set up a new bond potential for C-C bonds and add it to the potential set.
bond_potential = HarmonicBondPotential(
    particleType1='C',
    particleType2='C',
    k=13.44287*eV/Ang**2,
    r_0=1.526*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='C',
    particleType2='H',
    k=14.74379*eV/Ang**2,
    r_0=1.090*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='H',
    particleType2='C',
    particleType3='C',
    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='H',
    particleType2='C',
    particleType3='H',
    k=1.51774*eV,
    theta0=1.9111355,
)
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='H',
    particleType2='C',
    particleType3='C',
    particleType4='H',
    k=[0.0067455*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)
