# Set up lattice
lattice = FaceCenteredCubic(5.4306*Angstrom)

# Define elements
elements = [Silicon, Silicon]

# Define coordinates
fractional_coordinates = [[0.0,  0.0,  0.0 ],
                          [0.25, 0.25, 0.25]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
)

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

# Set up bonds for bonded force fields
bulk_configuration.findBonds(fuzz_factor=1.2)

potentialSet = TremoloXPotentialSet(name='VFF_silicon')
potentialSet.addParticleType(ParticleType(
    symbol='Si',
    mass=28.0855*atomic_mass_unit,
    charge=None,
    sigma=None,
    sigma14=None,
    epsilon=None,
    epsilon14=None,
    atomicNumber=14,
    tags=[]
))

# Add a VFFBondStrechingPotential for Si-Si bonds.
potential = VFFBondStretchingPotential(
    particleType1=ParticleIdentifier('Si'),
    particleType2=ParticleIdentifier('Si'),
    alpha=0.205289014474*eV/Angstrom**4,
    delta=5.52964056752*Angstrom**2,
    A=0.0*1/Angstrom**2,
    epsilon=0.0*Angstrom**2,
)
potentialSet.addPotential(potential)

# Add a VFFBondBendingPotential for Si-Si-Si angles.
potential = VFFBondBendingPotential(
    particleType1=ParticleIdentifier('Si'),
    particleType2=ParticleIdentifier('Si',),
    particleType3=ParticleIdentifier('Si',),
    alpha=0.0584121324987*eV/Angstrom**4,
    delta=1.84321352251*Angstrom**2,
)
potentialSet.addPotential(potential)
calculator = TremoloXCalculator(parameters=potentialSet)

bulk_configuration.setCalculator(calculator)