# %% (6,0) Carbon Nanotube

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

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

# Define coordinates
fractional_coordinates = [[ 0.659867810732,  0.5           ,  0.666666666667],
                          [ 0.638449585341,  0.579933905366,  0.833333333333],
                          [ 0.659867810732,  0.5           ,  0.333333333333],
                          [ 0.579933905366,  0.638449585341,  0.666666666667],
                          [ 0.5           ,  0.659867810732,  0.833333333333],
                          [ 0.638449585341,  0.579933905366,  0.166666666667],
                          [ 0.579933905366,  0.638449585341,  0.333333333333],
                          [ 0.420066094634,  0.638449585341,  0.666666666667],
                          [ 0.361550414659,  0.579933905366,  0.833333333333],
                          [ 0.5           ,  0.659867810732,  0.166666666667],
                          [ 0.420066094634,  0.638449585341,  0.333333333333],
                          [ 0.340132189268,  0.5           ,  0.666666666667],
                          [ 0.361550414659,  0.420066094634,  0.833333333333],
                          [ 0.361550414659,  0.579933905366,  0.166666666667],
                          [ 0.340132189268,  0.5           ,  0.333333333333],
                          [ 0.420066094634,  0.361550414659,  0.666666666667],
                          [ 0.5           ,  0.340132189268,  0.833333333333],
                          [ 0.361550414659,  0.420066094634,  0.166666666667],
                          [ 0.420066094634,  0.361550414659,  0.333333333333],
                          [ 0.579933905366,  0.361550414659,  0.666666666667],
                          [ 0.638449585341,  0.420066094634,  0.833333333333],
                          [ 0.5           ,  0.340132189268,  0.166666666667],
                          [ 0.579933905366,  0.361550414659,  0.333333333333],
                          [ 0.638449585341,  0.420066094634,  0.166666666667]]

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


# %% Set ForceFieldCalculator

# %% ForceFieldCalculator
potentialSet = Tersoff_C_2010()
calculator = TremoloXCalculator(parameters=potentialSet)

# %% Set Calculator
bulk_configuration.setCalculator(calculator)
bulk_configuration.update()

# %% OptimizeGeometry
optimized_configuration = OptimizeGeometry(
    configuration=bulk_configuration,
    max_forces=0.01*eV/Angstrom,
    constraints=[
        BravaisLatticeConstraint()
    ],
    trajectory_filename='CNT_results.hdf5'
)

nlsave('CNT_results.hdf5', optimized_configuration)


# %% DynamicalMatrix
dynamical_matrix = DynamicalMatrix(
    configuration=optimized_configuration,
    filename='CNT_results.hdf5',
    object_id='dm',
    calculator=calculator,
    repetitions=(1, 1, 5)
)
dynamical_matrix.update()


# %% ProjectedPhononDensityOfStates
qpoints = MonkhorstPackGrid(
    na=1,
    nb=1,
    nc=1000
)

projectedphonondensityofstates = ProjectedPhononDensityOfStates(
    dynamical_matrix=dynamical_matrix,
    qpoints=qpoints,
    energies=numpy.linspace(0.0, 210.0, 1001)*meV,
    spectrum_method=GaussianBroadening(
        broadening=0.5*meV
    )
)
nlsave('CNT_results.hdf5', projectedphonondensityofstates)