In the third example of sscha, the symmetry of the q-point does not match in the third step.

[WanshengBai]

Dear Developers,

I hope this message finds you well. I am writing to report an issue I encountered while performing SSCHA calculations, specifically during the temperature cycling step where I'm facing a q-point symmetry mismatch problem.

I have modified the force calculation interface in my program to use the QE interface. Below is my complete program:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
#  SSCHA_exercise_PhaseTransition.py
#

import os
import sys
import warnings

# Ignore Julia warnings
warnings.filterwarnings("ignore", message="Julia extension not available")

# ==================== Direct monkey patch solution ====================

print("Applying direct monkey patch...")

import cellconstructor.calculators

# Save original method
original_calculate = cellconstructor.calculators.Espresso.calculate

# Create global counter
_qe_counter = 0
_qe_main_dir = "qe_three"

# Create main directory
if not os.path.exists(_qe_main_dir):
    os.makedirs(_qe_main_dir)
    print(f"Created main directory: {_qe_main_dir}")

def patched_calculate(self, structure):
    """Patched calculate method"""
    global _qe_counter, _qe_main_dir
    _qe_counter += 1
    
    # Save original directory
    original_dir = self.directory
    
    # Create organized directory path
    qe_dir = os.path.join(_qe_main_dir, f"qe{_qe_counter}")
    if not os.path.exists(qe_dir):
        os.makedirs(qe_dir)
    
    print(f"[Direct patch] Created directory for calculation #{_qe_counter}: {qe_dir}")
    
    # Set new directory
    self.set_directory(qe_dir)
    
    try:
        # Call original calculate method
        return original_calculate(self, structure)
    except Exception as e:
        print(f"QE calculation failed: {e}")
        # Create error marker file
        with open(os.path.join(qe_dir, "ERROR.txt"), "w") as f:
            f.write(f"Calculation failed: {e}\n")
        raise
    finally:
        # Restore original directory
        self.set_directory(original_dir)

# Apply patch
cellconstructor.calculators.Espresso.calculate = patched_calculate

print("✓ Direct monkey patch applied successfully")

# ==================== Original script content ====================

# Import the cellconstructor stuff
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.ForceTensor
import cellconstructor.Structure
import cellconstructor.Spectral

# Import the modules to run the sscha
import sscha, sscha.Ensemble, sscha.SchaMinimizer
import sscha.Relax, sscha.Utilities

import spglib
from ase.visualize import view
from ase.data import atomic_numbers  # Added this import

# Import Matplotlib to plot
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import timeit

# ==================== QE Calculator Configuration ====================

# QE calculator configuration
pseudo = {
    "Sn": "Sn.pbe-dn-rrkjus_psl.1.0.0.UPF",
    "Te": "Te.pbe-n-rrkjus_psl.1.0.0.UPF"
}

input_data = {
    "control": {
        "calculation": 'scf',
        "pseudo_dir": '/netdisk/home/wshbai/pbe',
        "verbosity": 'high',
        "disk_io": "none",
        "tstress": True,
        "tprnfor": True
    },
    "system": {
        "occupations": "smearing",
        "smearing": "gauss",
        "ecutwfc": 82,
        "ecutrho": 820,
        "degauss": 0.001
    },
    "electrons": {
        "conv_thr": 1e-6
    }
}

# Use original calculator directly (patch already applied)
calc = CC.calculators.Espresso(
    input_data=input_data,
    pseudopotentials=pseudo,
    command="/public/software/intel/2020.4/intelpython3/bin/mpirun -np 24 /netdisk/public/software/espresso/7.0/pw.x -in PREFIX.pwi > PREFIX.pwo",
    kpts=(2, 2, 2),
    koffset=(0, 0, 0)
)

# Set calculator label and directory
calc.set_label("snte_sscha")
calc.set_directory(".")

# ==================== Main logic of original script ====================

# Setting the variables:
# Setting the temperature in Kelvin:
Temperature = 0
# Setting the number of configurations:
configurations = 10
# Setting the names and location of the files:
Files_dyn_SnTe = "Rh3Tl2S2.dyn"
# Set the number of irreducible q (related to the supercell size):
nqirr = 3
# Setting the frequencies output file:
File_frequencies = "frequencies.dat"
# Setting the dynamical matrix output filename:
File_final_dyn = "final_sscha_T{}_".format(int(Temperature))
sobol = False
sobol_scatter = False

# Load the dynamical matrix for the force field
ff_dyn = CC.Phonons.Phonons("Rh3Tl2S2.dyn", 3)

# Note: We are using real QE calculator instead of ToyModel
# Therefore, we no longer need to set ff_calculator, but use the calc defined above

# Define the temperatures, from 50 to 300 K, 6 temperatures
temperatures = np.linspace(50, 300, 6)

lowest_hessian_mode = []
lowest_sscha_mode = []

# Perform a simulation at each temperature
t_old = Temperature

print(f"Note: All QE calculation directories will be organized under '{_qe_main_dir}' directory")

for Temperature in temperatures:
    print(f"\n=== Starting calculation for temperature {Temperature} K ===")
    
    # Load the starting dynamical matrix
    dyn = CC.Phonons.Phonons(File_final_dyn.format(int(t_old)), nqirr)

    # Prepare the ensemble
    ensemble = sscha.Ensemble.Ensemble(dyn, T0 = Temperature,
                                  supercell = dyn.GetSupercell())

    # Prepare the minimizer
    minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
    minim.min_step_struc = 0.05
    minim.min_step_dyn = 0.002
    minim.kong_liu_ratio = 0.5
    minim.meaningful_factor = 0.000001
    #minim.root_representation = "root4"
    #minim.precond_dyn = False
    #minim.minim_struct = True
    #minim.neglect_symmetries = True
    minim.enforce_sum_rule = True  # Lorenzo's solution to the error

    # Prepare the relaxer (through many population)
    # Use real QE calculator instead of ToyModel
    relax = sscha.Relax.SSCHA(minim, ase_calculator=calc,
                  N_configs=configurations, max_pop=5)

    # Relax
    relax.relax(sobol=sobol, sobol_scramble=sobol_scatter)

    # Save the dynamical matrix
    relax.minim.dyn.save_qe(File_final_dyn.format(int(Temperature)))

    # Fixed part: correctly get atomic numbers
    structure = relax.minim.dyn.structure
    lattice = structure.unit_cell
    positions = structure.coords
    
    # Method 1: if structure.atoms are atomic symbol strings
    if isinstance(structure.atoms[0], str):
        numbers = [atomic_numbers[atom] for atom in structure.atoms]
    # Method 2: if structure.atoms are already atom objects
    else:
        numbers = [atom.Z for atom in structure.atoms]
    
    # Call spglib with correct format
    symm = spglib.get_spacegroup((lattice, positions, numbers), 0.005)
    print('Current SG = ', symm,' at T=',int(Temperature))

    # Recompute the ensemble for the hessian calculation
    ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0 = Temperature,
                            supercell = dyn.GetSupercell())
    ensemble.generate(configurations, sobol=sobol, sobol_scramble=sobol_scatter)
    
    # Use real QE calculator to compute energies and forces
    ensemble.get_energy_forces(calc, compute_stress=False)
    # gets the energies and forces from calc

    # update weights!!!
    ensemble.update_weights(relax.minim.dyn, Temperature)
    # Get the free energy hessian
    dyn_hessian = ensemble.get_free_energy_hessian(include_v4=False)
    # free energy hessian as in Bianco paper 2017
    dyn_hessian.save_qe("hessian_T{}_".format(int(Temperature)))

    # Get the lowest frequencies for the sscha and the free energy hessian
    w_sscha, pols_sscha = relax.minim.dyn.DiagonalizeSupercell() # dynamical matrix
    # Get the structure in the supercell
    superstructure = relax.minim.dyn.structure.generate_supercell(relax.minim.dyn.GetSupercell())

    # Discard the acoustic modes
    acoustic_modes = CC.Methods.get_translations(pols_sscha, superstructure.get_masses_array())
    w_sscha = w_sscha[~acoustic_modes]

    lowest_sscha_mode.append(np.min(w_sscha) * CC.Units.RY_TO_CM) # Convert from Ry to cm-1

    w_hessian, pols_hessian = dyn_hessian.DiagonalizeSupercell() # recomputed dyn for hessian
    # Discard the acoustic modes
    acoustic_modes = CC.Methods.get_translations(pols_hessian, superstructure.get_masses_array())
    w_hessian = w_hessian[~acoustic_modes]
    lowest_hessian_mode.append(np.min(w_hessian) * CC.Units.RY_TO_CM) # Convert from Ry to cm-1

    t_old = Temperature

# Prepare file to save results
freq_data = np.zeros( (len(temperatures), 3))
freq_data[:, 0] = temperatures
freq_data[:, 1] = lowest_sscha_mode
freq_data[:, 2] = lowest_hessian_mode

# Save results to file
np.savetxt("{}_hessian_vs_temperature.dat".format(configurations),
            freq_data, header = "T [K]; SSCHA mode [cm-1]; Free energy hessian [cm-1]")

hessian_data = np.loadtxt("{}_hessian_vs_temperature.dat".format(configurations))

plt.figure(dpi = 120)
plt.plot(hessian_data[:,0], hessian_data[:,1], label = "Min SCHA freq", marker = ">")
plt.plot(hessian_data[:,0], hessian_data[:,2], label = "Free energy curvature", marker = "o")
plt.axhline(0, 0, 1, color = "k", ls = "dotted") # Draw the zero
plt.xlabel("Temperature [K]")
plt.ylabel("Frequency [cm-1]")
plt.legend()
plt.tight_layout()
plt.savefig('{}_Temp_Freq.png'.format(configurations))
#plt.show()

plt.figure(dpi = 120)
plt.plot(hessian_data[:,0], np.sign(hessian_data[:,2]) * hessian_data[:,2]**2,
                    label = "Free energy curvature", marker = "o")
plt.axhline(0, 0, 1, color = "k", ls = "dotted") # Draw the zero
plt.xlabel("Temperature [K]")
plt.ylabel("$\omega^2$ [cm-2]")
plt.legend()
plt.tight_layout()
plt.savefig('{}_Temp_Omeg.png'.format(configurations))
plt.show()

print("SSCHA phase transition calculation completed!")

When I run the program, I encounter the following error:

# ---------------- NEW MINIMIZATION STEP --------------------
Step ka =  416
Minimization step, force computed: 10
Good step found with 1.1763612623183063e-05, try increment


Number of symmetries before the step:  48
Harmonic contribution to free energy =     -94.94131462 meV
Anharmonic contribution to free energy = -2451338.80946862 +-   69506.43374791 meV
Free energy = -2451433.75078324 +-   69506.43374791 meV
FC gradient modulus = 2012236.60426818 +-   11782.88016138 bohr^2
Struct gradient modulus =     165.99226433 +- 982835024505.53881836 meV/A
Kong-Liu effective sample size =  9.763199576132429


The gw gradient satisfy the convergence condition.
Check the stopping criteria: Running =  True

 # ---------------- NEW MINIMIZATION STEP --------------------
Step ka =  417
Minimization step, force computed: 10
Total frequencies (excluding translations):
0) -0.0000  | 1.9597 cm-1
1) -0.0000  | 1.9597 cm-1
2) -0.0000  | 1.9597 cm-1
3) 35.0149  | 37.0251 cm-1
4) 35.0149  | 37.0251 cm-1
5) 35.0149  | 37.0251 cm-1
6) 35.0149  | 37.0251 cm-1
7) 35.0149  | 37.0251 cm-1
8) 35.0149  | 37.0251 cm-1
9) 44.7433  | 40.5105 cm-1
10) 44.7433  | 40.5105 cm-1
11) 44.7433  | 40.5105 cm-1
12) 58.0392  | 57.5593 cm-1
13) 58.0392  | 57.5593 cm-1
14) 58.0392  | 57.5593 cm-1
15) 79.3637  | 75.9498 cm-1
16) 79.3637  | 75.9498 cm-1
17) 79.3637  | 75.9498 cm-1
18) 79.3637  | 75.9498 cm-1
19) 79.3637  | 75.9498 cm-1
20) 79.3637  | 75.9498 cm-1
21) 78.3740  | 78.7802 cm-1
22) 78.3740  | 78.7802 cm-1
23) 78.3740  | 78.7802 cm-1
24) 78.3740  | 78.7802 cm-1
25) 78.3740  | 78.7802 cm-1
26) 78.3740  | 78.7802 cm-1
27) 78.3740  | 78.7802 cm-1
28) 78.3740  | 78.7802 cm-1
29) 103.4367  | 103.4885 cm-1
30) 103.4367  | 103.4885 cm-1
31) 103.4367  | 103.4885 cm-1
32) 103.4367  | 103.4885 cm-1
33) 103.4367  | 103.4885 cm-1
34) 103.4367  | 103.4885 cm-1
35) 103.4367  | 103.4885 cm-1
36) 103.4367  | 103.4885 cm-1
37) 106.8183  | 106.5371 cm-1
38) 106.8183  | 106.5371 cm-1
39) 106.8183  | 106.5371 cm-1
40) 106.8183  | 106.5371 cm-1
41) 119.4381  | 118.8126 cm-1
42) 119.4381  | 118.8126 cm-1
43) 119.4381  | 118.8126 cm-1
44) 119.4381  | 118.8126 cm-1

Immaginary frequencies found! Redoing the step.
Step too large (scalar = 21873.499921534683 | kl_ratio = 0.0), reducing to 1.541469327374884e-05
Total frequencies (excluding translations):
0) -0.0000  | 1.9597 cm-1
1) -0.0000  | 1.9597 cm-1
2) -0.0000  | 1.9597 cm-1
3) 35.0449  | 37.0251 cm-1
4) 35.0449  | 37.0251 cm-1
5) 35.0449  | 37.0251 cm-1
6) 35.0449  | 37.0251 cm-1
7) 35.0449  | 37.0251 cm-1
8) 35.0449  | 37.0251 cm-1
9) 44.6710  | 40.5105 cm-1
10) 44.6710  | 40.5105 cm-1
11) 44.6710  | 40.5105 cm-1
12) 58.0213  | 57.5593 cm-1
13) 58.0213  | 57.5593 cm-1
14) 58.0213  | 57.5593 cm-1
15) 79.2957  | 75.9498 cm-1
16) 79.2957  | 75.9498 cm-1
17) 79.2957  | 75.9498 cm-1
18) 79.2957  | 75.9498 cm-1
19) 79.2957  | 75.9498 cm-1
20) 79.2957  | 75.9498 cm-1
21) 78.3803  | 78.7802 cm-1
22) 78.3803  | 78.7802 cm-1
23) 78.3803  | 78.7802 cm-1
24) 78.3803  | 78.7802 cm-1
25) 78.3803  | 78.7802 cm-1
26) 78.3803  | 78.7802 cm-1
27) 78.3803  | 78.7802 cm-1
28) 78.3803  | 78.7802 cm-1
29) 103.4373  | 103.4885 cm-1
30) 103.4373  | 103.4885 cm-1
31) 103.4373  | 103.4885 cm-1
32) 103.4373  | 103.4885 cm-1
33) 103.4373  | 103.4885 cm-1
34) 103.4373  | 103.4885 cm-1
35) 103.4373  | 103.4885 cm-1
36) 103.4373  | 103.4885 cm-1
37) 106.8136  | 106.5371 cm-1
38) 106.8136  | 106.5371 cm-1
39) 106.8136  | 106.5371 cm-1
40) 106.8136  | 106.5371 cm-1
41) 119.4257  | 118.8126 cm-1
42) 119.4257  | 118.8126 cm-1
43) 119.4257  | 118.8126 cm-1
44) 119.4257  | 118.8126 cm-1

Immaginary frequencies found! Redoing the step.
Step too large (scalar = 21873.499921534683 | kl_ratio = 0.0), reducing to 1.346597491406021e-05


Number of symmetries before the step:  48
Harmonic contribution to free energy =     -94.84309440 meV
Anharmonic contribution to free energy = -2452753.80682993 +-   68290.21210504 meV
Free energy = -2452848.64992433 +-   68290.21210504 meV
FC gradient modulus = 2012236.60426818 +-   11782.88016138 bohr^2
Struct gradient modulus =     165.99226433 +- 982835024505.53881836 meV/A
Kong-Liu effective sample size =  9.726418653553724


The gw gradient satisfy the convergence condition.
Check the stopping criteria: Running =  True

 # ---------------- NEW MINIMIZATION STEP --------------------
Step ka =  418
Minimization step, force computed: 10
Good step found with 1.346597491406021e-05, try increment


Number of symmetries before the step:  48
Harmonic contribution to free energy =     -94.84309440 meV
Anharmonic contribution to free energy = -2452753.80682993 +-   68290.21210504 meV
Free energy = -2452848.64992433 +-   68290.21210504 meV
FC gradient modulus = 1962767.89559535 +-   11782.88016138 bohr^2
Struct gradient modulus =    1049.82725851 +- 5401008940212.83789062 meV/A
Kong-Liu effective sample size =  9.726418653553724


The gw gradient satisfy the convergence condition.
Check the stopping criteria: Running =  True

 # ---------------- NEW MINIMIZATION STEP --------------------
Step ka =  419
Minimization step, force computed: 10
Reciprocal lattice vectors:
[[-0.15651077  0.15651077  0.15651077]
 [ 0.15651077 -0.15651077  0.15651077]
 [ 0.15651077  0.15651077 -0.15651077]]
Passed q star:
[[-0.07825539 -0.07825539  0.07825539]
 [ 0.07825539 -0.07825539 -0.07825539]
 [-0.07825539  0.07825539 -0.07825539]
 [ 0.07825539  0.07825539  0.07825539]]
QE q star:
[[-0.07825539 -0.07825539  0.07825539]
 [-0.07825539  0.07825539  0.07825539]]
Traceback (most recent call last):
  File "/netdisk/home/wshbai/pro/test-sscha/SnTe/sscha/10.py", line 196, in <module>
    relax.relax(sobol=sobol, sobol_scramble=sobol_scatter)
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/sscha/Relax.py", line 389, in relax
    self.minim.run(custom_function_pre = self.__cfpre__,
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/sscha/SchaMinimizer.py", line 1300, in run
    timer.execute_timed_function(self.minimization_step, custom_function_gradient)
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/cellconstructor/Timer.py", line 112, in execute_timed_function
    ret = function(*args, timer=new_timer,**kwargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/sscha/SchaMinimizer.py", line 642, in minimization_step
    timer.execute_timed_function(self.dyn.Symmetrize, use_spglib = self.use_spglib)
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/cellconstructor/Timer.py", line 112, in execute_timed_function
    ret = function(*args, timer=new_timer,**kwargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/cellconstructor/Phonons.py", line 3281, in Symmetrize
    timer.execute_timed_function(qe_sym.SymmetrizeFCQ, fcq, self.q_stars, asr = asr, verbose = verbose)
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/cellconstructor/Timer.py", line 114, in execute_timed_function
    ret = function(*args, **kwargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/cellconstructor/symmetries.py", line 1048, in SymmetrizeFCQ
    self.ApplyQStar(fcq[q0_index : q0_index + q_len, :,:], np.array(q_stars[i]))
  File "/netdisk/home/wshbai/.conda/envs/sscha/lib/python3.11/site-packages/cellconstructor/symmetries.py", line 838, in ApplyQStar
    raise ValueError("Error, the passed q star does not match the one computed by QE")
ValueError: Error, the passed q star does not match the one computed by QE

The error occurs during the symmetry application step in the SSCHA minimization process. From the error message, I can see that the program expects a q-star with 4 q-points, but QE computes only 2 q-points, leading to a mismatch.

I would be extremely grateful if you could help me resolve this issue and provide any suggestions for improving the calculation.

Thank you for your time and assistance.

Sincerely,
Wansheng Bai

[mesonepigreco]

The fact that this happens mid-minimization, rather than at the first step, suggests that something bizarre occurred during the minimization. Notably, this occurs after finding imaginary frequencies. The very strange part is that some of the frequencies are very close to 0 cm-1.

This likely occurs because you have a very small number of configurations (10 in this case), which, especially if the system has an instability, is too little at low temperatures.

Also, you are computing the Hessian with the same 10 configurations, which will likely give you random noise (you probably need at least several hundred).

The easiest solution is to increase the number of configurations to 50/100, which are more reasonable numbers for the minimization, then increase to 500-1000 for the free energy hessian. Another possible solution is to minimize higher temperatures before, as their minimization is usually more stable.

[WanshengBai]

Dear SSCHA Developers,

Hello!

First and foremost, I sincerely thank you for your valuable suggestions. However, I still have some uncertainties regarding the implementation. My understanding is that to address the issue of insufficient q-point sampling, adjustments can be made in two aspects:

1. During the initial free energy minimization step, instead of starting at 0 K, it might be better to use a slightly higher initial temperature.
2. In the final temperature cycling step, increasing the number of sampled configurations may help improve convergence.

Regarding the second point, I would like to confirm: Is increasing the number of configurations achieved precisely in the temperature cycling step by adjusting the configurations parameter, as shown in the following code snippet?

# Recompute the ensemble for the hessian calculation  
ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0 = Temperature,  
                            supercell = dyn.GetSupercell())  
ensemble.generate(configurations, sobol=sobol, sobol_scramble=sobol_scatter)  

If possible, I would greatly appreciate it if you could provide more detailed explanations on the following two points:

• What considerations should be taken into account when selecting a non-zero initial temperature for the initial free energy minimization step?
• What is the specific method for correctly increasing the number of configurations in the temperature cycling step, and how does it impact convergence?

Your professional advice is crucial for optimizing our computational workflow and deepening our understanding of the methodology.

Once again, thank you for your efforts and time. I look forward to your reply.

Best regards,
Wansheng Bai