Regarding obtaining detailed information about the new structure

[Prutthi]

Dear SSCHA Expert,

I am reaching out to convey my enthusiastic interest in SSCHA. I am relatively new to using the SSCHA code, and I have recently completed an example involving the strucure relaxation with temperature of an Au. I have some inquiries regarding obtaining detailed information about the new structure resulting from the calculation:

1. How can I retrieve comprehensive details about the new structure?
2. When SSCHA performs calculations on a supercell structure, is the reported new structure detail based on the supercell structure or the unit cell structure?

Your guidance and insights on these matters would be greatly appreciated.

Sincerely,

Prutthipong

[mesonepigreco]

Hi Prutthipong,
We are very happy to see your enthusiasm! Thanks for reaching out.

The informations are stored in the dynamical matrices files.
In the tutorial it is saved on a file with the following line:

relax.minim.dyn.save_qe("sscha_T300_dyn")

These files (13 in case of gold on a 4x4x4 supercell) contain the structure in the primitive cell dynamical matrix in Fourier space on a 4x4x4 q-mesh. The 13 files are different points in the 4x4x4 q-mesh that are inequivalent by symmetry (there is no symmetry that convert one of this 13 q points into the other). This 13 number changes depending on symmetries, supercell and structure.

You can load the final dynamical matrix and extract the structure (in the primitive cell) with

NQIRR = 13
dyn = CC.Phonons.Phonons("sscha_T300_dyn", NQIRR)

structure = dyn.structure

You can then print the volume, the coordinates of the atoms with:

print("Volume:", structure.get_volume())
print("Atomic positions:", structure.coords)

You can even convert the structure into an ASE atoms and visualize in 3D:

import ase, ase.visualize
ase_structure = structure.get_ase_atoms()
ase.visualize.view(ase_structure)

You can save as a cif using ASE or any other format you like and parse and analyze also with different softwares.

Regarding the supercell: the calculation is performed in the supercell, however the average positions of the atoms is still defined in the primitive cell, therefore, the dynamical matrix stores the average positions of atoms (in the primitive cells). To generate the structure in the supercell you can use:

supercell_structure = structure.generate_supercell((4,4,4)) # get a 4x4x4 supercell
ase.visualize.view(supercell_structure.get_ase_atoms())

To extract a random ensemble at the correct temperature (with atoms in the dispaced configurations that sample thermal and quantum fluctuations) you must use the dynamical matrix dyn we loaded from the output of the SSCHA caculation:

random_structures = dyn.ExtractRandomStructures(size=10, T = 1000) # Extract 10 structures in the supercell at 1000 K
random_struct = random_structures[0]
ase.visualize.view(random_struct.get_ase_atoms())

Then you can use these structure to do extra calculation (like the bandgap that will be renormalized with electron-phonon coupling effects).

Hope this helps,
Bests,
Lorenzo

[Prutthi]

Hi Lorenzo,

I appreciate your assistance. Currently, I have implemented your suggestion for the following step: "NQIRR = 13 dyn = CC.Phonons.Phonons("sscha_T300_dyn", NQIRR)."

Nevertheless, I would appreciate clarification on how to select the suitable value for NQIRR, which in this case is set to 13, in the context of the SSCHA code. To be more specific, I'm uncertain about the prerequisite steps, such as whether I should carry out electron-phonon or phonon-phonon calculations using the Quantum espresso code, in order to determine the appropriate value of NQIRR before advancing with the SSCHA code.

Best,

Prutthipong

[mesonepigreco]

Hi, nqirr is just the number of dynamical matrices printed in output and saved. It is cellconstructor when computing the phonons that automatically calculate the nqirr value in this line of the tutorial (at the beginning):

gold_harmonic_dyn = CC.Phonons.compute_phonons_finite_displacements(gold_structure_relaxed, calculator, supercell = (4,4,4))

There the nqirr is fixed by the structure and the supercell. To see in the code what it is you can ask the length of the q_stars:

print("nqirr = ", len(gold_harmonic_dyn.q_stars))

If you use quantum-espresso to compute the harmonic phonons ab initio (with the ph.x program, for example), then the nqirr is printed in the dynamical matrix ending with 0, together with the list of the irreducible q-points.

The only prerequisite step of a SSCHA calculation is having the starting dynamical matrices, those can either be obtained directly from cellconstructor (as in the line above that computes the dynamical matrix of gold) or by quantum espresso employing the ph.x program.
In both cases, nqirr is given in output.

[Prutthi]

Hi, I want to express my sincere gratitude for your valuable suggestions and support. Your guidance has significantly enhanced my comprehension of how to utilize the SSCHA code.

[mesonepigreco]

Yes of course,

You need to use your phonon calculations as a starting guess for the dynamical matrix and then you can relax the structure using temperature, and analyze if imaginary phonons remains even after you consider thermal fluctuations. However, the SSCHA relaxation will keep the symmetries of the original calculation, so if your instability is locked by a symmetry, the SSCHA will not relax in that direction (but you can study if the instability disappears with temperature by computing the phonons as a function of temperature).
To relax the instability, you should distort a bit your original structure to break the symmetries.

Here there is a tutorial where we show how to investigate a phase transition as a function of temperature on a structure that has imaginary phonons:

• http://sscha.eu/Tutorials/tutorial_03_secondorder_phase_transitions/

This is a short version of the full SnTe example we reported in the paper of the code (Section 6 of https://iopscience.iop.org/article/10.1088/1361-648X/ac066b/meta ). Here we fully characterize the phase transition of this toy model material both from the high-symmetry phase and the low symmetry phase.

If you want a deeper theoretical background, there is also the video lecture on Second Order phase transitions by Raffaello Bianco here, http://sscha.eu/lectures/

All the best,
Lorenzo

[Prutthi]

Hi Lorenzo,

I deeply appreciate your continued support and for sharing your publications and video lectures. Your explanations have significantly improved my understanding.

However, I still have a lingering question regarding the Au example. Specifically, I'm interested in the atomic positions and lattice parameters. To explore this, I used the following code to print the volume and atomic coordinates:

print("Volume:", structure.get_volume())
print("Atomic positions:", structure.coords)
Upon examining the results of this example, I have observed a significant shift in the atomic positions of Au, and it appears that there is an overlap of Au atoms in this context:

Volume: 72.57895213056294
Atomic positions: [[0. 0. 0. ]
[0. 2.08564426 2.08564426]
[2.08564426 0. 2.08564426]
[2.08564426 2.08564426 0. ]]
Interestingly, the lattice parameter appears to have remained almost unchanged. In my understanding, when a structure is computed under temperature conditions, the lattice parameter should experience thermal expansion. I'm uncertain about what might be missing in my calculations.

If you wouldn't mind, could you kindly provide some guidance and insights regarding this matter?

Best regards,

Prutthipong

[mesonepigreco]

Hi, indeed you should see the thermal expansion.

In particular, the line:

# Run the NVT simulation
relax.relax(get_stress = True)

Keeps the lattice fixed. It is always good, especially if starting from the harmonic case, to first relax the system at fixed volume, and then relax also the volume.

The part of the code that relax the volume at target pressure is:

relax.vc_relax(target_press = 0)

vc_relax means "variable cell" relaxation, coherently as it is in quantum espresso.
This case we are constraining the pressure to 0 GPa (ambient pressure).
The SSCHA also allows to perform a variable cell relaxation bug fixing the total volume (changing only the angles and relative size of lattice parameters), with the following:

relax.vc_relax(fix_volume=True)

In the case of gold, however, since the cell is cubic, the only free parameter of the lattice is the total volume, therefore this command is equivalent to a relax.

Are you using vc_relax to simulate the system? In principle, if you relax the lattice at several temperatures you should get something like the image shown in the tutorial:

Regarding your structure, you can visualize it with:

import ase, ase.visualize
ase.visualize.view(structure.get_ase_atoms())

Remember to extract the structure from the minimization:

structure = relax.minim.dyn.structure

Hope this helps,
Best,
Lorenzo

[Prutthi]

Hi Lorenzo,

I want to express my sincere gratitude for your prompt response and the invaluable support you've provided. I'm enthusiastic about implementing your suggestions.

Regarding the previous question about the noticeable shift in atomic positions compared to the initial structure, I have observed that employing the conventional cell of Au may not yield satisfactory results. However, when using the primitive cell of Au, I found that the atomic positions of Au matched those of the initial structure (primitive cell). This leads me to believe that employing the primitive cell is highly advantageous for the SSCHA code, if I've understood correctly.

With this newfound insight, and in light of your mention of an instability being locked by symmetry, it appears that SSCHA will not relax in that specific direction. This seems to be the root cause of the overlap of Au atoms in certain contexts. Hence, it is indeed crucial that a crystal structure generated using the SSCHA code does not remain bound by symmetry constraints. Have I interpreted this correctly?

Best regards,
Prutthipong

[mesonepigreco]

Hi, in general symmetry constrains helps the SSCHA to optimize the structure much more effectively, so it is a good thing.
You can understand a posteriori if some of the systems wants to escape from the high-symmetry solution by evaluating the free energy Hessian and looking for imaginary phonons.

If the free energy hessian phonons are positive definite in the whole spectrum, then the symmetries choosen for the SSCHA are good.

To simulate the free energy hessian, there is the tutorial on the website:
Tutorial Free Energy Hessian

Indeed, you can always displace the atoms to lose all symmetries and let the SSCHA optimize the configurations.
However, take into account that you may need to increase the number of configurations with lower symmetries, and the final results is affected by stochastic noise (SSCHA is employing a Monte Carlo technique to optimize the structure).
This means that the atomic positions have an error associated to them that depends on the number of configurations and symmetries.

In your case, it is very strange that Au atoms overlap.
For example, in the output you attached in a previous response you have:

Atomic positions: [[0. 0. 0. ]
[0. 2.08564426 2.08564426]
[2.08564426 0. 2.08564426]
[2.08564426 2.08564426 0. ]]

Here, atoms are not overlapping, but are in the standard FCC configurations.
They do not move during the optimization since everything is locked by symmetries (regardless if you use the primitive or the conventional cell).
For example if I ask cellconstructor to print the atomic coordinates in crystal units, I get:

>>>  CC.Methods.covariant_coordinates(structure.unit_cell, structure.coords) 
array([[0. , 0. , 0. ],
       [0. , 0.5, 0.5],
       [0.5, 0. , 0.5],
       [0.5, 0.5, 0. ]])

As you can see, the relative position with respect to the unit cell vector is still the one fixed by symmetry.
If you see a shift in the absolute value of the coordinates, it is because the volume and lattice parameter is changing, not because the atoms moved.

If you want to brake the symmetries, you can manually shift an atomic position (then regenerate the stars of q points in the dynamical matrix with the new symmetry). Here an example:

# Shift an atom to violate the FCC symmetries
dyn.structure.coords[0, 0] += 0.05
# Reorganize the star of q points with the new symmetries
dyn.AdjustQStar()

However, I would break the symmetries only if the free energy Hessian displays an imaginary mode, as this relaxation will be computationally more expensive.

Bests,
Lorenzo