# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Molecule Configuration
# -------------------------------------------------------------

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

# Define coordinates
cartesian_coordinates = [[ 0.309666700862,  0.758525407442, -0.437934848174],
                         [ 1.739666598276,  0.758525407442, -0.437934848174],
                         [-0.16699993161 ,  0.758525407442, -1.786151680918],
                         [-0.476666632472, -1.167590026966,  0.674108416372],
                         [-0.16699993161 ,  1.926115434409,  0.236173568198],
                         [-1.906666529886, -1.167590026966,  0.674108416372],
                         [ 0.            , -2.335180053932, -0.            ],
                         [ 0.            , -1.167590026966,  2.022325249116]]*Angstrom

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

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

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


_potential = TersoffBrennerPairPotential(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('Si', []),
	A = 1830.8*eV,
	B = 471.18*eV,
	l = 2.4799*1/Angstrom,
	mu = 1.7322*1/Angstrom,
	Re = 2.35*Angstrom,
	R1 = 2.7*Angstrom,
	R2 = 3.0*Angstrom,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerPairPotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('H', []),
	A = 80.07*eV,
	B = 31.38*eV,
	l = 4.2075*1/Angstrom,
	mu = 1.7956*1/Angstrom,
	Re = 0.74*Angstrom,
	R1 = 1.1*Angstrom,
	R2 = 1.7*Angstrom,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerPairPotential3(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('H', []),
	A = 323.54*eV,
	B = 84.18*eV,
	l = 2.9595*1/Angstrom,
	mu = 1.6158*1/Angstrom,
	Re = 1.475*Angstrom,
	R1 = 1.7*Angstrom,
	R2 = 2.0*Angstrom,
	xRep1 = numpy.array([1, 2, 3, 4, 5, 6, 7]),
	FRep1 = numpy.array([ 2.01 ,  2.218,  1.908,  2.   ,  2.   ,  2.   ,  2.   ]),
	dxFRep1 = numpy.array([ 0.   , -0.038, -0.154,  0.   ,  0.   ,  0.   ,  0.   ]),
	dxxFRep1 = numpy.array([ 1.324, -1.4  ,  1.168,  0.   ,  0.   ,  0.   ,  0.   ]),
	xRep2 = numpy.array([ 0.,  1.]),
	FRep2 = numpy.array([ 0.,  0.]),
	dxFRep2 = numpy.array([ 0.,  0.]),
	dxxFRep2 = numpy.array([ 0.,  0.]),
	xAttr1 = numpy.array([1, 2, 3, 4, 5, 6, 7]),
	FAttr1 = numpy.array([ 1.86 ,  2.07 ,  1.868,  2.   ,  2.   ,  2.   ,  2.   ]),
	dxFAttr1 = numpy.array([-0.    ,  0.0204, -0.0576,  0.    ,  0.    ,  0.    ,  0.    ]),
	dxxFAttr1 = numpy.array([ 1.2192, -1.1784,  1.0224,  0.    ,  0.    ,  0.    ,  0.    ]),
	xAttr2 = numpy.array([ 0.,  1.]),
	FAttr2 = numpy.array([ 0.,  0.]),
	dxFAttr2 = numpy.array([ 0.,  0.]),
	dxxFAttr2 = numpy.array([ 0.,  0.]),
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerBOPairPotential(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('Si', []),
	delta = 0.635,
	eta = 0.78734,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerBOPairPotential(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('H', []),
	delta = 0.80469,
	eta = 1.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerBOPairPotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('Si', []),
	delta = 0.80469,
	eta = 1.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerBOPairPotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('H', []),
	delta = 0.80469,
	eta = 1.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('Si', []),
	particleType3 = ParticleIdentifier('Si', []),
	beta = 3,
	alpha = 5.1975*1/Angstrom**3,
	g_c = 0.0,
	g_d = 0.16,
	g_h = -0.59826,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('H', []),
	particleType3 = ParticleIdentifier('H', []),
	beta = 1,
	alpha = 3.0*1/Angstrom,
	g_c = 4.0,
	g_d = 0.0,
	g_h = 0.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('H', []),
	particleType3 = ParticleIdentifier('Si', []),
	beta = 1,
	alpha = 3.0*1/Angstrom,
	g_c = 4.0,
	g_d = 0.0,
	g_h = 0.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('Si', []),
	particleType3 = ParticleIdentifier('H', []),
	beta = 1,
	alpha = 3.0*1/Angstrom,
	g_c = 4.0,
	g_d = 0.0,
	g_h = 0.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential(
	particleType1 = ParticleIdentifier('H', []),
	particleType2 = ParticleIdentifier('Si', []),
	particleType3 = ParticleIdentifier('Si', []),
	beta = 1,
	alpha = 0.0*1/Angstrom,
	g_c = 0.7,
	g_d = 1.0,
	g_h = -1.0,
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential4(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('Si', []),
	particleType3 = ParticleIdentifier('H', []),
	beta = 3,
	alpha = 4.0*1/Angstrom**3,
	g_c = 0.0216,
	g_d = 0.27,
	x = numpy.array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.]),
	H = numpy.array([-0.04 , -0.04 , -0.276, -0.47 , -0.47 , -0.47 , -0.47 ]),
	dxH = numpy.array([-0.    , -0.1028, -0.2968,  0.    ,  0.    ,  0.    ,  0.    ]),
	dxxH = numpy.array([ 0.2056, -0.4112,  0.0232, -0.    ,  0.    ,  0.    ,  0.    ]),
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential4(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('H', []),
	particleType3 = ParticleIdentifier('Si', []),
	beta = 3,
	alpha = 4.0*1/Angstrom**3,
	g_c = 0.0216,
	g_d = 0.27,
	x = numpy.array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.]),
	H = numpy.array([-0.04 , -0.04 , -0.276, -0.47 , -0.47 , -0.47 , -0.47 ]),
	dxH = numpy.array([-0.    , -0.1028, -0.2968,  0.    ,  0.    ,  0.    ,  0.    ]),
	dxxH = numpy.array([ 0.2056, -0.4112,  0.0232, -0.    ,  0.    ,  0.    ,  0.    ]),
)
potentialSet.addPotential(_potential)
_potential = TersoffBrennerTriplePotential4(
	particleType1 = ParticleIdentifier('Si', []),
	particleType2 = ParticleIdentifier('H', []),
	particleType3 = ParticleIdentifier('H', []),
	beta = 3,
	alpha = 4.0*1/Angstrom**3,
	g_c = 0.0216,
	g_d = 0.27,
	x = numpy.array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.]),
	H = numpy.array([-0.04 , -0.04 , -0.276, -0.47 , -0.47 , -0.47 , -0.47 ]),
	dxH = numpy.array([-0.    , -0.1028, -0.2968,  0.    ,  0.    ,  0.    ,  0.    ]),
	dxxH = numpy.array([ 0.2056, -0.4112,  0.0232, -0.    ,  0.    ,  0.    ,  0.    ]),
)
potentialSet.addPotential(_potential)
calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setInternalOrdering("default")
calculator.setVerletListsDelta(0.25*Angstrom)

molecule_configuration.setCalculator(calculator)
molecule_configuration.update()