# -------------------------------------------------------------
# Molecule configuration, CH4
# -------------------------------------------------------------

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

# Define coordinates
cartesian_coordinates = [[ 11.61454603,  11.09244804,   9.99999993],
                         [ 10.67862603,  10.93144595,  10.52879495],
                         [  9.99999997,  11.75542099,  10.3236729 ],
                         [ 10.22904002,  10.00000001,  10.1942581 ],
                         [ 10.87074303,  10.87768505,  11.59731907]]*Angstrom

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
calculator = BrennerCalculator()

molecule_configuration.setCalculator(calculator)
molecule_configuration.update()

molecule_configuration = OptimizeGeometry(
        molecule_configuration,
        max_forces=0.001*eV/Ang,
        max_steps=200,
        max_step_length=0.5*Ang,
        trajectory_filename=None,
        optimize_cell=False,
        optimizer_method=LBFGS(),
        )

# -------------------------------------------------------------
# Vibrational modes
# -------------------------------------------------------------
# Compute Hessian by finite difference of forces.
dynamical_matrix = DynamicalMatrix(molecule_configuration, 'ch4.hdf5', object_id='dynamical_matrix')

# Vibrational analysis.
vib_mode = VibrationalMode(dynamical_matrix, mode_indices=range(6,15))
vib_mode.nlprint()