Author Topic: Phonon Band Structure Calculations-Negative Frequency Band Structure  (Read 3724 times)

0 Members and 1 Guest are viewing this topic.

Offline zhenlan

  • Regular QuantumATK user
  • **
  • Posts: 8
  • Country: dk
  • Reputation: 0
    • View Profile
Hi,
I am currently working on calculating the phonon band structure for a particular system. However, I encountered a discrepancy between the results obtained using ATK (QuantumATK) and VASP:

ATK Calculation: The phonon band structure shows negative frequencies, suggesting dynamic instability.
VASP Calculation: The phonon band structure does not exhibit these negative frequencies.
For your reference, I have attached the phonon band structures generated by both ATK and VASP. Additionally, I have included the input script used for the ATK calculations below:

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

# Set up lattice
lattice = Hexagonal(6.53017*Angstrom, 21.0*Angstrom)

# 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,
            Carbon, Carbon, Carbon, Carbon]

# Define coordinates
fractional_coordinates = [[ 0.      ,  0.      ,  0.420238],
                          [ 0.095238,  0.47619 ,  0.420238],
                          [ 0.285714,  0.428571,  0.420238],
                          [ 0.523809,  0.619047,  0.420238],
                          [ 0.142857,  0.714286,  0.420238],
                          [ 0.380952,  0.904761,  0.420238],
                          [ 0.571429,  0.857143,  0.420238],
                          [ 0.809524,  1.047619,  0.420238],
                          [ 0.238095,  1.190476,  0.420238],
                          [ 0.428571,  1.142857,  0.420238],
                          [ 0.666667,  1.333334,  0.420238],
                          [ 0.857143,  1.285714,  0.420238],
                          [ 0.714286,  1.571429,  0.420238],
                          [ 0.952381,  1.761905,  0.420238],
                          [ 0.      ,  0.      ,  0.579762],
                          [ 0.190476,  0.238095,  0.579762],
                          [ 0.142857,  0.428571,  0.579762],
                          [ 0.333333,  0.666666,  0.579762],
                          [ 0.571429,  0.714286,  0.579762],
                          [ 0.761905,  0.952381,  0.579762],
                          [ 0.047619,  0.809524,  0.579762],
                          [ 0.285714,  0.857143,  0.579762],
                          [ 0.47619 ,  1.095238,  0.579762],
                          [ 0.714286,  1.142857,  0.579762],
                          [ 0.904762,  1.380953,  0.579762],
                          [ 0.428571,  1.285714,  0.579762],
                          [ 0.619048,  1.52381 ,  0.579762],
                          [ 0.857143,  1.571429,  0.579762]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
k_point_sampling = MonkhorstPackGrid(
    na=7,
    nb=7,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=50.0*Hartree,
    k_point_sampling=k_point_sampling,
    occupation_method=MethfesselPaxton(1000.0*Kelvin*boltzmann_constant, 1),
    )

algorithm = PulayMixer(
    restart_strategy=AdaptiveHistoryRestart(
        effective_rank_fraction=0.0,
        ),
    )

iteration_control_parameters = IterationControlParameters(
    tolerance=1e-06,
    algorithm=algorithm,
    )

#----------------------------------------
# Grimme DFTD3
#----------------------------------------
correction_extension = GrimmeDFTD3(
    exchange_correlation=GGA.PBE,
    maximum_neighbour_distance=30.0*Ang,
    include_three_body_term=False,
    )

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

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('TBG-21-relax.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------

constraints = [FixStrain(x=False, y=False, z=True)]

bulk_configuration = OptimizeGeometry(
    bulk_configuration,
    max_forces=0.001*eV/Ang,
    max_stress=0.1*GPa,
    max_steps=200,
    max_step_length=0.2*Ang,
    constraints=constraints,
    trajectory_filename='TBG-21-relax_trajectory.hdf5',
    trajectory_interval=5.0*Minute,
    restart_strategy=RestartFromTrajectory(),
    optimizer_method=LBFGS(),
    enable_optimization_stop_file=True,
)
nlsave('TBG-21-relax.hdf5', bulk_configuration)
nlprint(bulk_configuration)


# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Analysis from File
# -------------------------------------------------------------
path = '/home/energy/zhenlan/ATK/TBG/21/TBG-21-relax.hdf5'
configuration = nlread(path, object_id='BulkConfiguration_1')[0]

# -------------------------------------------------------------
# Dynamical Matrix
# -------------------------------------------------------------
dynamical_matrix = DynamicalMatrix(
    configuration,
    filename='TBG-21-dynmat.hdf5',
    object_id='dynamical_matrix',
    repetitions=(3, 3, 1),
    use_symmetry=False,
    atomic_displacement=0.01*Angstrom,
    acoustic_sum_rule=True,
    finite_difference_method=Central,
    force_tolerance=1e-08*Hartree/Bohr**2,
    processes_per_displacement=4,
    log_filename_prefix='forces_displacement_',
    use_wigner_seitz_scheme=False,
    )
dynamical_matrix.update()

# -------------------------------------------------------------
# Phonon Bandstructure
# -------------------------------------------------------------
phonon_bandstructure = PhononBandstructure(
    configuration=configuration,
    dynamical_matrix=dynamical_matrix,
    route=['G', 'M', 'K', 'G'],
    points_per_segment=100,
    number_of_bands=All
    )
nlsave('TBG-21-dynmat.hdf5', phonon_bandstructure)

If anyone has encountered similar issues or has insights on how to align the settings between ATK and VASP more closely, I would greatly appreciate your input.

Best regards,
ZY

Offline Anders Blom

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 5575
  • Country: dk
  • Reputation: 96
    • View Profile
    • QuantumATK at Synopsys
Replied by email, hopefully some suggestions worked.

1)   Run the same calculation with a simple Forcefield, which will take just a few minutes or even less.
2)   Run it with PlaneWave for a more direct comparison to VASP. More time-consuming for sure though…
3)   We don’t have much experience with Methfessel-Paxton; I doubt it would be a major problem, but perhaps at least use a lower temperature like 300 K, and perhaps just stick to FermiDirac.
4)   When experimenting with various settings, use a cheap DFT LCAO basis set like SingleZetaPolarized to make it run faster; any problems with negative frequencies are likely to show up there as well as when using default PseudoDojo Medium.