RAMAN Spectrum on monolayer MoS2

[AsifShah]

Dear Admin, I am trying to calculate Raman Spectrum of monolayer MoS2. But I only see one peak around 396 cm^-1 whereas the experimental peaks are around 385 and 403 cm^-1. Kindly clarify the issue. I also tried changing polarization direction to 100. That gives three peaks but incorrect also. The script is as follows:

Code

# %% Molybdenite (1)

# Set up lattice
lattice = Hexagonal(3.1604*Angstrom, 40.0*Angstrom)

# Define elements
elements = [Molybdenum, Sulfur, Sulfur]

# Define coordinates
fractional_coordinates = [[ 0.333333333333,  0.666666666667,  0.5           ],
                          [ 0.666666666667,  0.333333333333,  0.460348625   ],
                          [ 0.666666666667,  0.333333333333,  0.539651375   ]]

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

molybdenite_1_name = "Molybdenite (1)"


# %% Set LCAOCalculator

# %% LCAOCalculator

# ----------------------------------------
# Basis Set
# ----------------------------------------
basis_set = [
    BasisGGASG15.Sulfur_High,
    BasisGGASG15.Molybdenum_High,
]

k_point_sampling = MonkhorstPackGrid(na=13, nb=13)

numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=120.0 * Hartree, k_point_sampling=k_point_sampling
)

calculator = LCAOCalculator(
    basis_set=basis_set,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    checkpoint_handler=NoCheckpointHandler,
)


# %% Set Calculator

molybdenite_1.setCalculator(calculator)

molybdenite_1.update()

nlsave('Raman_SG15High.hdf5', molybdenite_1)


# %% OptimizeGeometry

restart_strategy = RestartFromTrajectory(
    trajectory_filename='Raman_SG15High.hdf5',
    object_id='optimize_trajectory',
)

optimized_configuration = OptimizeGeometry(
    configuration=molybdenite_1,
    max_forces=0.01 * eV / Angstrom,
    constraints=[FixStrain(True, True, True)],
    trajectory_filename='Raman_SG15High.hdf5',
    trajectory_object_id='optimize_trajectory',
    optimize_cell=False,
    restart_strategy=restart_strategy,
)

nlsave('Raman_SG15High.hdf5', optimized_configuration, object_id='optgeom')


# %% DynamicalMatrix

dynamical_matrix = DynamicalMatrix(
    configuration=optimized_configuration,
    filename='Raman_SG15High.hdf5',
    object_id='dm',
    calculator=calculator,
)
dynamical_matrix.update()


# %% SusceptibilityDerivatives

kpoints = MonkhorstPackGrid(na=23, nb=23, nc=2)

susceptibility_derivatives = SusceptibilityDerivatives(
    configuration=optimized_configuration,
    filename='Raman_SG15High.hdf5',
    object_id='sd',
    kpoints=kpoints,
)

susceptibility_derivatives.update()


# %% RamanSpectrum

raman_spectrum = RamanSpectrum(
    configuration=optimized_configuration,
    dynamical_matrix=dynamical_matrix,
    susceptibility_derivatives=susceptibility_derivatives,
    polarization_in=(0.0, 0.0, 1.0),
    polarization_out=(0.0, 0.0, 1.0),
    number_of_angles=None,
    polarization_orthogonal=None,
)
nlsave('Raman_SG15High.hdf5', raman_spectrum)

[Anders Blom]

Sorry, but it's an impossible question to answer short of myself running all calculations to check that every parameter is using correct sensible values (enough k-points, maybe needs more?), if the structure is correctly optimized (why not optimize the cell size?), etc, etc.

[AsifShah]

Ah, just follow up. It turned out I had to choose the polarization average rather than polarization in/out. That worked, qualitatively atleast.