constrain the SSCHA simulation to the R3m symmetry

[ZhichaoLiu-0306]

Dear developer，

I try to reproduce Figure 3 in your new paper https://iopscience.iop.org/article/10.1088/1361-648X/ac066b/meta. I have a problem with constraining symmetry. The initial dynamical matrix "ffield_dynq*" with the symmetry of Fm-3m (225) is in Tutorials-sscha\SnTe_ToyModel. The cell parameters of SnTe (R3m: 160) are from the vc-relax by using QE. However, after running the commands "AdjustToNewCell" and "GetStrainMatrix", the symmetry of R-3m (166) is obtained, which is different from the symmetry of R3m I aimed. So, I am confused in the application method of "AdjustToNewCell" and "GetStrainMatrix". How can I constrain the system to the R3m symmetry if I want to use the SnTe toy model force field?

Here the script is:

############################################################################################
from __future__ import print_function
from __future__ import division
import sys,os
import ase
from numpy import *
import numpy as np
import cellconstructor as CC
import cellconstructor.Structure
import cellconstructor.Phonons
import cellconstructor.symmetries as SYM
import cellconstructor.Methods as Methods
from ase.calculators.espresso import Espresso
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities
import fforces as ff
import fforces.Calculator
import spglib
ff_dyn = CC.Phonons.Phonons("ffield_dynq", 3)
ff_dyn.ForcePositiveDefinite()
ff_dyn.Symmetrize()
ff_calculator = ff.Calculator.ToyModelCalculator(ff_dyn)
ff_calculator.type_cal = "pbtex"
ff_calculator.p3 = 0.036475
ff_calculator.p4 = -0.022
ff_calculator.p4x = -0.014
print("ff_dyn sym:   ", spglib.get_spacegroup(ff_dyn.structure.get_ase_atoms()))

alat = 8.62054978 #(Bohr unit  ecutwfc 50 ecutrho 500 "gaussian" 1D-2) (Sn.pbe-dn-kjpaw_psl.1.0.0.UPF  Te.pbe-n-kjpaw_psl.1.0.0.UPF)
R3m_unit_cell = alat * np.array([[0.00598958, 0.70709410, 0.70709410],
                                 [0.70709410, 0.00598958, 0.70709410],
                                 [0.70709410, 0.70709410, 0.00598958]])
R3m = CC.Structure.Structure(2)
R3m.atoms = ["Sn", "Te"]
R3m.coords[0,:] = np.array([-0.0038462659, -0.0038462659, -0.0038462659]) # in cartesian axis
R3m.coords[1,:] = np.array([0.6841114166, 0.6841114166, 0.6841114166])
R3m.unit_cell = R3m_unit_cell
R3m.has_unit_cell = True
R3m.masses = {"Sn": 108197.545983522, 
              "Te": 116300.285296078}   # Ry units
R3msuper=R3m.generate_supercell([2, 2, 2])
print("aimed sym:   ", spglib.get_spacegroup(R3msuper.get_ase_atoms()))

ff_dyn.AdjustToNewCell(R3m.unit_cell, symmetrize = True)
print("operate 1:   ", spglib.get_spacegroup(ff_dyn.structure.get_ase_atoms()))
dyn_sscha = ff_dyn.GetStrainMatrix(R3m.unit_cell, 0)
print("operate 2:   ", spglib.get_spacegroup(dyn_sscha.structure.get_ase_atoms()))
dyn_sscha.save_qe("begin_sscha_T0_")

Best regards
Z. Liu

[mesonepigreco]

Hi Z. Liu,
AdjustToUnitCell changes only the unit cell (adjusting the q points and the cartesian vectors). However, since the original ff_dyn was a cubic structure, the atomic positions are not those of the R3m structure. To get the R3m symmetry, you should also change the atomic coordinates to those of your cell.

Inserting the line

ff_dyn.structure.coords[:,:] = R3m.coords[:,:]

After

ff_dyn.AdjustToNewCell(R3m.unit_cell, symmetrize = True)

You will get the correct symmetry group.

Bests,
Lorenzo

[chenyubi14]

Dear Lorenzo,

I am reaching out to inquire further about the symmetry constraint in your paper (https://iopscience.iop.org/article/10.1088/1361-648X/ac066b/meta). You demonstrated the application of finite temperature structure search, leading to the discovery of a new phase with Cc symmetry.

However, we can do even more: finite temperature structure search.
To find the real ground state, we release all the symmetry constrains in our simulation and perform a full relaxation with the SSCHA at T = 100 K. We discovered a new phase of Cc symmetry defined in a 1 × 2 × 1 supercell of the original cubic cell.

I am curious about how to release the symmetry constraint. In my understanding, the symmetry has been preserved during minimization. The above code ff_dyn.AdjustToNewCell(R3m.unit_cell, symmetrize = True) assigns symmetry R3m, which is maintained after minimization. My questions are:

1. Could you provide me with an example code of releasing symmetry constraints?
2. If no symmetry is constrained, how did you eventually get Cc symmetry, rather than C1 no symmetry? Do you get it by spglib.get_spacegroup()?

Thanks,
Yubi

[mesonepigreco]

Hi,
To check if more symmetries are gained during a SSCHA optimization you can use spglib.get_spacegroup as you mentioned. Indeed, you need to change the symmetry threshold otherwise you will always obtain the lowest constrained symmetries.

A quantitative test that can be done to assert if the symmetry group is just close or exactly the one you see with a higher symmetry threshold is that the distance from the symmetry group should go to zero as 1/sqrt(N) where N is the number of configurations employed in the SSCHA minimization.

For this purpose, you can get the symmetrized structure:

# Get the new structure
THRESHOLD = 0.05
new_structure = minim.dyn.structure.copy()
spglib_symmetries = spglib.get_symmetry(new_structure.get_ase_atoms(), THRESHOLD)
symmetries = CC.symmetries.GetSymmetriesFromSPGLIB(spglib_symmetries)
new_structure.impose_symmetries(symmetries)

# Check the distance with the new structure
distance = np.linalg.norm( (new_structure.coords - minim.dyn.structure.coords).flatten()) 

You can iterate the last part by repeating the last minimization with less configurations

TOTAL_CONFIGS = 1024
N_CFGS = [32, 64, 128, 256, 512, 1024] # Assuming an ensemble of 1024 configurations
dist_data = np.zeros((len(N_CFGS), 2))
for i, ncfgs in enumerate(N_CFGS):
    mask = np.zeros(TOTAL_CONFIGS, dtype=bool)
    mask[:ncfgs] = True
    new_ensemble = ensemble.split(mask)

    # Redefine the minimization
    minim = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
    minim.meaningful_factor = 0.001
    minim.init()
    minim.run()

    distance = np.linalg.norm( (new_structure.coords - minim.dyn.structure.coords).flatten())
    dist_data[i, 0] = ncfgs
    dist_data[i, 1] = distance

np.save_txt("distance_from_high_symmetry.txt", dist_data, header = "N configs; distance [A]")

Then you can analyze if the data in distance_from_high_symmetry decrease as 1/sqrt(N).
Please, note that I just wrote this script on the GitHub, so it may need some tweaks to properly work.

[chenyubi14]

Dear Lorenzo,

Thank you for your prompt response. It is very helpful to check the symmetry. However, I'm still unclear about the method for releasing the symmetry constraints during the simulation.

(1) Do I need to manually set the symmetry of ff_dyn.structure to be C1? Could you provide me with an example code?
(2) After turning off the symmetry, how can I determine whether the supercell structure is more stable?

I apologize if these questions seem basic; I'm eager to fully understand this process.

Thanks,
Yubi

[mesonepigreco]

Hi,

• (1) : that is good, but then you must adjust also the dynamical matrix as the number of irreducible q point changes:

dyn.structure.coords += np.random.normal(0.0, 0.001, size=dyn.structure.coords.shape)
dyn.AdjustQStar()

• (2) Just compare the free energies of the phases you want to compare (scaled on the same mole/number of atoms).

Bests,
Lorenzo