filename = 'ElectroOpticalTensorWorkflow.hdf5'

# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------

# Set up lattice
lattice = SimpleTetragonal(4.070385412109747*Angstrom, 4.115226913895633*Angstrom)

# Define elements
elements = [Barium, Titanium, Oxygen, Oxygen, Oxygen]

# Define coordinates
fractional_coordinates = [[ 0.028650982641, -0.028665746616,  0.087957144726],
                          [ 0.544026489212,  0.455959167129,  0.605611994272],
                          [ 0.515281929165,  0.484708002041,  1.05622848206 ],
                          [ 0.514348270569, -0.004594374366,  0.570135591311],
                          [ 0.00458017807 ,  0.48564221155 ,  0.570135494531]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
k_point_sampling = MonkhorstPackGrid(
    na=10,
    nb=10,
    nc=10,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=105.0*Hartree,
    k_point_sampling=k_point_sampling,
    occupation_method=FermiDirac(100.0*Kelvin*boltzmann_constant),
    )

iteration_control_parameters = IterationControlParameters(
    tolerance=1e-05,
    max_steps=1000,
    )

calculator = LCAOCalculator(
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    iteration_control_parameters=iteration_control_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave(filename, bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
constraints = [SpaceGroupConstraint()]

bulk_configuration = OptimizeGeometry(
    bulk_configuration,
    max_forces=0.01*eV/Ang,
    max_stress=0.05*GPa,
    max_steps=2000,
    max_step_length=0.2*Ang,
    constraints=constraints,
    trajectory_filename=None,
    optimizer_method=LBFGS(),
)
nlsave(filename, bulk_configuration)
nlprint(bulk_configuration)

# -------------------------------------------------------------
# Optical Spectrum
# -------------------------------------------------------------
kpoint_grid = MonkhorstPackGrid(
    na=25,
    nb=25,
    nc=25,
    )

optical_spectrum = OpticalSpectrum(
    configuration=bulk_configuration,
    kpoints=kpoint_grid,
    energies=numpy.linspace(0, 0, 1)*eV,
    broadening=0.1*eV,
    bands_below_fermi_level=1000,
    bands_above_fermi_level=1000,
    )
nlsave(filename, optical_spectrum)

# -------------------------------------------------------------
# Born Effective Charge
# -------------------------------------------------------------
born_effective_charge = BornEffectiveCharge(
    configuration=bulk_configuration,
    kpoints_a=MonkhorstPackGrid(30, 10, 10),
    kpoints_b=MonkhorstPackGrid(10, 30, 10),
    kpoints_c=MonkhorstPackGrid(10, 10, 30),
    atomic_displacement=0.01*Angstrom,
    )
nlsave(filename, born_effective_charge)

# -------------------------------------------------------------
# Dynamical Matrix
# -------------------------------------------------------------
polar_phonon_splitting_parameters = PolarPhononSplittingParameters(
    optical_spectrum=optical_spectrum, born_effective_charge=born_effective_charge)
dynamical_matrix = DynamicalMatrix(
    bulk_configuration,
    filename=filename,
    object_id='dynamical_matrix',
    repetitions=(5, 5, 5),
    atomic_displacement=0.01*Angstrom,
    acoustic_sum_rule=True,
    finite_difference_method=Central,
    force_tolerance=1e-08*Hartree/Bohr**2,
    processes_per_displacement=1,
    log_filename_prefix='forces_displacement_',
    use_wigner_seitz_scheme=False,
    polar_phonon_splitting_parameters=polar_phonon_splitting_parameters,
    )
dynamical_matrix.update()

# -------------------------------------------------------------
# Dielectric Tensor
# -------------------------------------------------------------
dielectric_tensor = DielectricTensor(
    configuration=bulk_configuration,
    dynamical_matrix=dynamical_matrix,
    optical_spectrum=optical_spectrum,
    born_effective_charge=born_effective_charge,
)
nlsave(filename, dielectric_tensor)

# -------------------------------------------------------------
# Susceptibility Derivatives
# -------------------------------------------------------------
susceptibility_derivatives = SusceptibilityDerivatives(
    bulk_configuration,
    filename=filename,
    object_id='susceptibility_derivatives',
    atomic_displacement=0.01*Angstrom,
    finite_difference_method=Central,
    broadening=0.1*eV,
    bands_below_fermi_level=1000,
    bands_above_fermi_level=1000,
    processes_per_displacement=1,
    log_filename_prefix='susceptibility_displacement_',
)
susceptibility_derivatives.update()

# -------------------------------------------------------------
# Raman Spectrum
# -------------------------------------------------------------
raman_spectrum = RamanSpectrum(
    bulk_configuration,
    dynamical_matrix,
    susceptibility_derivatives,
    qpoint=[0.0, 0.0, 0.0],
    phonon_modes=None,
    method=PolarizationDependent,
    polarization_in=[1, 0, 0],
    polarization_out=[1, 0, 0],
    broadening=0.25 * meV,
    )
nlsave(filename, raman_spectrum)

# -------------------------------------------------------------
# Electro Optical Tensor
# -------------------------------------------------------------
electro_optical_tensor = ElectroOpticalTensor(
    configuration=bulk_configuration,
    dielectric_tensor=dielectric_tensor,
    raman_spectrum=raman_spectrum,
    second_harmonics_generation_susceptibility=None,
)
nlsave(filename, electro_optical_tensor)
nlprint(electro_optical_tensor)