Puzzling results of phonon calculations of gold ( QE + SSCHA)

[hellolori]

Dear developers,

I've been diligently following the instructions to familiarize myself with SSCHA. I successfully replicated the outcomes presented in the "2.1 The free energy of gold: a simulation in the NVT ensemble" section of the "HANDS-ON-SESSION 1 - FIRST SSCHA SIMULATIONS: FREE ENERGY AND STRUCTURAL RELAXATIONS" chapter. I was able to reproduce both "Fig. 2.1.2: Minimization data of Gold" and "Fig. 2.1.3: All the SSCHA phonon frequencies as a function of the step in the NVT simulation," which are illustrated below:

Then, I switch from using the "ASE" calculator to "Quantum Espresso," as outlined in the "2.3 Ab initio calculation with the SSCHA code" section. I aimed to calculate the phonon dispersion of gold at 300 K. However, I've encountered some puzzling results. Although the calculations ran to completion, the outcomes appear to be quite unreasonable. For your reference, here is a snapshot of the unexpected results:

It's evident that the free energy increases gradually with each step of free energy minimization. Additionally, the phonon frequencies exhibit a significant softening trend (right panel). This is particularly evident in the phonon dispersion diagram displayed below. I've observed that the phonon frequencies derived from the QE+SSCHA approach are notably diminished when compared to those obtained using ASE+SSCHA and harmonic calculations.

The above results lead me to believe that there might be an error or issue in my calculation process. I would greatly appreciate any suggestions or guidance you can provide.

I've attached the relevant input files and the sequence of steps taken during the calculations for your review:

(1)nvt_para.py

# Import the sscha code
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities

# Import the cellconstructor library to manage phonons
import cellconstructor as CC, cellconstructor.Phonons
import cellconstructor.Structure, cellconstructor.calculators

# Import the DFT calculator
import cellconstructor.calculators
from ase.calculators.espresso import Espresso
# Import numerical and general pourpouse libraries
import numpy as np, matplotlib.pyplot as plt
import sys, os

# Initialize the DFT (Quantum Espresso) calculator for gold
# The input data is a dictionary that encodes the pw.x input file namelist
input_data = {
    'control' : {
        'disk_io' : 'None',
        'pseudo_dir' : '.',
        "tstress" : True
    },
    'system' : {
        'ecutwfc' : 45,   
        'ecutrho' : 45*4, 
        'occupations' : 'smearing',
        'smearing' : 'mv',
        'degauss' : 0.03
    },
    'electrons' : {
        'conv_thr' : 1e-7
    }
}

pseudopotentials = {'Au' : 'Au_ONCV_PBE-1.0.upf'}
# the kpoints mesh and the offset
kpts = (1,1,1)
koffset = (1,1,1)

# Specify the command to call quantum espresso
command = '~/software/openmpi/4.0.5/bin/mpiexec -n 24 ~/software/QE/6.6_openmpi/bin/pw.x  -npool 3  -i  PREFIX.pwi > PREFIX.pwo'
# Prepare the quantum espresso calculator
calculator = CC.calculators.Espresso(input_data,
                                     pseudopotentials,
                                     command = command,
                                     kpts = kpts,
                                     koffset = koffset)

TEMPERATURE = 300
N_CONFIGS = 50
MAX_ITERATIONS = 20
START_DYN = 'harmonic_dyn'
NQIRR = 13

# Let us load the starting dynamical matrix
gold_dyn = CC.Phonons.Phonons(START_DYN, NQIRR)

# Initialize the random ionic ensemble
ensemble = sscha.Ensemble.Ensemble(gold_dyn, TEMPERATURE)

# Initialize the free energy minimizer
minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minim.set_minimization_step(0.01)

# Initialize the NVT simulation
relax = sscha.Relax.SSCHA(minim, calculator, N_configs = N_CONFIGS,
                          max_pop = MAX_ITERATIONS)

# Define the I/O operations
ioinfo = sscha.Utilities.IOInfo()
ioinfo.SetupSaving("minim_info")
relax.setup_custom_functions(custom_function_post = ioinfo.CFP_SaveAll)

# Run the NVT simulation (save the stress to compute the pressure)
relax.relax(get_stress = True)

# If instead you want to run a NPT simulation, use
# The target pressure is given in GPa.
#relax.vc_relax(target_press = 0)

# You can also run a mixed simulation (NVT) but with variable lattice parameters
#relax.vc_relax(fix_volume = True)

# Now we can save the final dynamical matrix
# And print in stdout the info about the minimization
relax.minim.finalize()
relax.minim.dyn.save_qe("sscha_T{}_dyn".format(TEMPERATURE))

The above script is built based on "nvt_local.py" in the tutorial. I just add a targe "tstress" to " input_data", change the "command" to a parallel version, and use "CC.calculators.Espresso" instead of the "Espresso".

(2)Execution of the script

mpiexec -n 5   python  nvt_para.py >out.log

(3) The last section of the file "out.log" （The file "out.log" is also attached at the end of this message）

I eagerly await your response and suggestions on how to address this matter.

Thank you in advance for your assistance.

Best regards,

out.log.txt

[hellolori]

Thank you for your attention. I've resolved the issues.

[hellolori]

It seems that adding a tag "tprnfor = .true." to the input of QE helps.

[mesonepigreco]

Precisely.
The issue is that Quantum Espresso isn't displaying the forces in the output due to the unset tprnfor flag. As a result, when we use ASE to get structure.get_forces() for minimizing free energy, it gives us a value of 0. This means that during minimization (which relies on forces), the goal is to decrease phonon frequencies as much as possible to lower the zero-point energy. However, the actual free energy is calculated accurately (since Quantum Espresso provides energy values for structures), and it acknowledges the presence of potential energy, leading to a higher calculated value.

I hope this explanation proves helpful to others as well!

Best regards,
Lorenzo