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

# Set up lattice
vector_a = [3.84001408591, 0.0, 0.0]*Angstrom
vector_b = [0.0, 3.84001408591, 0.0]*Angstrom
vector_c = [0.0, 0.0, 5.4306]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

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

# Define coordinates
fractional_coordinates = [[ 0.5  ,  0.   ,  0.125],
                          [ 0.5  ,  0.5  ,  0.375],
                          [ 0.   ,  0.5  ,  0.625],
                          [ 0.   ,  0.   ,  0.875]]

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

# Add tags
bulk_configuration.addTags('doping_0')

# Add external potential
external_potential = AtomicCompensationCharge([
    ('doping_0', 0.000400390214212),
    ('doping_1', -0.000400390214212)
    ])

bulk_configuration.setExternalPotential(external_potential)

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

potentialSet = StillingerWeber_Si_1985()
calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setVerletListsDelta(0.25*Angstrom)

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('dynmat_bulk.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Dynamical Matrix
# -------------------------------------------------------------
dynamical_matrix = DynamicalMatrix(
    configuration=bulk_configuration,
    repetitions=(3, 3, 1),
    atomic_displacement=0.01*Angstrom,
    acoustic_sum_rule=True,
    symmetrize=True,
    finite_difference_method=Central,
    force_tolerance=1e-08*Hartree/Bohr**2,
    processes_per_displacement=1,
    log_filename_prefix='displacement_',
    use_wigner_seitz_scheme=False,
    )
nlsave('dynmat_bulk.hdf5', dynamical_matrix)

# -------------------------------------------------------------
# Vibrational Mode
# -------------------------------------------------------------
vibrational_mode = VibrationalMode(
    configuration=bulk_configuration,
    dynamical_matrix=dynamical_matrix,
    kpoint_fractional=[0, 0, 0],
    mode_indices=None,
    )
nlsave('dynmat_bulk.hdf5', vibrational_mode)