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

# Define coordinates
cartesian_coordinates = [[ 1.023885,  0.066068, -1.16103 ],
                         [ 0.063979, -0.004459, -0.657896],
                         [-0.781426, -0.083013, -1.33519 ],
                         [-0.063847,  0.004458,  0.657835],
                         [ 0.781272,  0.083007,  1.335485],
                         [-1.023863, -0.066061,  1.160796]]*Ang

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

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

# Set up improper dihedral indices for two SP2-carbon atoms.
improper_dihedral_indices = [[1, 0, 2, 3], [3, 4, 5, 1]]
molecule_configuration.setImproperDihedralIndices(improper_dihedral_indices)

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

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

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

# Add improper torsion potential around SP2-carbon atoms.
improper_torsion_potential = ImproperTorsionPotential(
    particleType1='C',
    particleType2='H',
    particleType3='H',
    particleType4='C',
    k=0.1*eV/Radians**2,
    psi0=0.0*Radians,
)
potential_set.addPotential(improper_torsion_potential)

# 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 cosine 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=2.0943951,
)
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.15719*eV],
    n=[2],
    delta=[3.14159]
)
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)