Application of Rotational Sum Rule for SnS monolayer

hiPhive
1

Dear Sir,

         I have calculated phonon dispersion of SnS monolayer using 5x5x1 supercell  with K mesh 5x5x1. I am getting negative frequencies around Gamma point. To remove the error I tried to apply the rotational sum rules for this structure. I have SPOSCAR, FORCE_CONSTANTS obtained from Phonopy, POSCAR files. I am using using the following script to obtain ForceConstantPotential:

…
from ase.io import read
from hiphive import ClusterSpace, ForceConstantPotential
from hiphive.cutoffs import estimate_maximum_cutoff
from hiphive.utilities import extract_parameters
from hiphive import ForceConstants, ClusterSpace, ForceConstantPotential
from hiphive import enforce_rotational_sum_rules
from hiphive.utilities import extract_parameters

prim = read(‘POSCAR’)
supercell = read(‘SPOSCAR’)

# Define a cluster space using the largest cutoff you can
max_cutoff = estimate_maximum_cutoff(supercell) - 0.01
cutoffs = [8.0] # only second order needed
cs = ClusterSpace(prim, cutoffs)

# import the phonopy force constants using the correct supercell also provided by phonopy
fcs =ForceConstants.read_phonopy(supercell,‘FORCE_CONSTANTS’)

# Find the parameters that best fits the force constants given you cluster space
parameters = extract_parameters(fcs, cs)

fcp = ForceConstantPotential(fcs, parameters)
…
The error I am getting is as follows:
…
Primitive cell:
** Formula: S2Sn2**
** Cell:**
** [ 4.16278 0.00000 0.00000]**
** [ 0.00000 4.36945 0.00000]**
** [ 0.00000 0.00000 20.45051]**
** Basis:**
** Sn [ 0.25000 0.38112 0.36012]**
** Sn [ 0.75000 0.88112 0.49733]**
** S [ 0.75000 0.97299 0.37231]**
** S [ 0.25000 0.47299 0.48514]**

Crystal symmetry:
** Spacegroup: Pmn2_1 (31)**
** Unique site: 2**
** Symmetry operations: 4**
** symprec: 1.00e-05**

Cutoffs:
** Maximum cutoff: 8.0**
** Found 4 center atoms with 48 images totaling 52 atoms**

Clusters:
** Clusters: {2: 162}**
** Total number of clusters: 162**

Orbits:
** Orbits: {2: 26}**
** Total number of orbits: 26**

Eigentensors:
** Eigentensors: {2: 198}**
** Total number of parameters: 198**

Constraints:
** Acoustic: True**
** Number of degrees of freedom: {2: 188}**
** Total number of degrees of freedom: 188**
Force constant reconstruction error order 2: 6.8425%
Traceback (most recent call last):
** File “/home/safikul_akd/Anaconda_physics/lib/python3.9/site-packages/hiphive/safikul.py”, line 23, in **
** fcp = ForceConstantPotential(fcs, parameters)**
** File “/home/safikul_akd/Anaconda_physics/envs/phonopy/lib/python3.10/site-packages/hiphive/force_constant_potential.py”, line 44, in init**
** self._prim = cs.primitive_structure.copy()**
AttributeError: ‘SortedForceConstants’ object has no attribute 'primitive_structure’
…
Sir, could you please help me by providing the full script to generate ForceConstantPotential and to apply Born-Huang rotational sum rules?

Kind Regards,
Safikul

2

This is commonly asked question and a good starting point is the answer given in the hiphive-FAQ

You are running fcp = ForceConstantPotential(fcs, parameters), but if you look at the basic tutorial or documentation for ForceConstantPotential you will see that it should be initialized with a ClusterSpace and parameters, e.g. fcp = ForceConstantPotential(cs, parameters).

3

Dear Sir,
I have generated fc2_rotinv.hdf5 file. Sir, how can I generate phonon dispersion curve using this file in phonopy?

Kind Regards,
Safikul

4

Dear Sir,

           Thank you for your help. I have solved this issue and plotted phonon dispersion curve and it looks fine. There is no negative frequencies around Gamma point.

Kind Regards,
Safikul

5

Dear Sir,

         I am using the following script to add rotational sum rules for SnS monolayer

…
from ase.io import read
from hiphive import ClusterSpace, ForceConstantPotential
from hiphive.cutoffs import estimate_maximum_cutoff
from hiphive.utilities import extract_parameters
from hiphive import ForceConstants, ClusterSpace, ForceConstantPotential
from hiphive import enforce_rotational_sum_rules
from hiphive.utilities import extract_parameters

prim = read(‘POSCAR’)
supercell = read(‘SPOSCAR’)

# Define a cluster space using the largest cutoff you can
max_cutoff = estimate_maximum_cutoff(supercell) - 0.01
cutoffs = [6.0] # only second order needed
cs = ClusterSpace(prim, cutoffs)

# import the phonopy force constants using the correct supercell also provided by phonopy
fcs = ForceConstants.read_phonopy(supercell, ‘FORCE_CONSTANTS’)

# Find the parameters that best fits the force constants given you cluster space
parameters = extract_parameters(fcs, cs)
# Enforce the rotational sum rules
parameters_rot = enforce_rotational_sum_rules(cs, parameters, [‘Huang’, ‘Born-Huang’])

# use the new parameters to make a fcp and then create the force constants and write to a phonopy file
fcp = ForceConstantPotential(cs, parameters_rot)
fcs = fcp.get_force_constants(supercell)
fcs.write_to_phonopy(‘FORCE_CONSTANTS_NEW’, format=‘text’)
…

Sir, I am generating FORCE_CONSTANTS_NEW file by this script. Sir, my question is when I am plotting the phonon dispersion curve with this new file then whole phonon spectra are getting modified. Why is this happening? I need to modify only the ZA accoustic branch near Gamma point which has negative frequencies and quadratic nature? So, how will I do that? Also I need to play with fitting parameters to get correct graph. Sir, could you please suggest how I should modify this script to add fitting parameters alpha?

Kind Regards,
Safikul

6

The reason why all modes change is because this is a post-processing method to adjust the force-constants so that they obey the rotational sum rules.

In your script you set cutoff = [6.0], this means that all force-constants corresponding to distances larger than 6.0Å will be set to zero. This will of course change the dispersion so ideally this cutoff should be set as large as possible. Try to increase this to the maximum allowed value and see if it helps.

If you have a look at the documentation for enforce rotational sum rules function you can see that it takes an alpha parameter as input. I suggest trying a few different values as done here and see if it helps.

Additionally there are some more examples regarding rotational sum rules in the hiphive-examples repo that might be of help.

7

Dear Sir,

         Increasing the cutoff=[ ] has no effect on correcting the band dispersion. Yet all the modes change. To use the condition fitting with constraints Sir I have followed the following scripts. But it gives me the following errors:

Error 1…
Traceback (most recent call last):
File “/home/safikul_akd/Anaconda_physics/lib/python3.9/site-packages/hiphive/iit.py”, line 5, in
from hiphive.fitting import Optimizer
ModuleNotFoundError: No module named ‘hiphive.fitting’
Error2…
Traceback (most recent call last):
File “/home/safikul_akd/Anaconda_physics/lib/python3.9/site-packages/hiphive/iit.py”, line 6, in
from tools import compute_dispersion
ImportError: cannot import name ‘compute_dispersion’ from ‘tools’ (/home/safikul_akd/Anaconda_physics/lib/python3.9/site-packages/tools/init.py)
…

Script
…
import numpy as np
import matplotlib.pyplot as plt
from hiphive import ForceConstantPotential, StructureContainer
from hiphive import enforce_rotational_sum_rules
from hiphive.fitting import Optimizer
from tools import compute_dispersion
import seaborn as sns
sns.set_context(‘talk’)

# parameters
lw = 2
alpha_transparent = 0.8
alpha_vals = np.logspace(1, 5, 25)

# fit model
sc = StructureContainer.read(‘POSCAR’)
cs = sc.cluster_space
opt = Optimizer(sc.get_fit_data(), train_size=1.0)
opt.train()
print(opt)

# get normal FCP
parameters = opt.parameters
fcp_normal = ForceConstantPotential(cs, parameters)
prim = fcp_normal.primitive_structure
q_normal, f_normal = compute_dispersion(prim, fcp_normal)

# apply sum rules
data = dict()
for alpha in alpha_vals:
** parameters_rot = enforce_rotational_sum_rules(cs, parameters, [‘Huang’, ‘Born-Huang’], alpha=alpha)**
** fcp = ForceConstantPotential(cs, parameters_rot)**
** q, f = compute_dispersion(prim, fcp)**
** data[alpha] = q, f**

# plotting
fig = plt.figure(figsize=(12, 5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

cmap = plt.get_cmap(‘plasma’)
colors = [cmap(i) for i in np.linspace(0, 1, len(data))]

for ax in [ax1, ax2]:
** for c, (alpha, (q, f)) in zip(colors, data.items()):**
** ax.plot(q, f, c=c, lw=lw, alpha=alpha_transparent)**
** ax.plot(q_normal, f_normal, ‘–k’)**
** ax.plot(np.nan, np.nan, ‘–k’, label=‘No rotational’)**

ax1.legend()
ax1.set_xlim([0.0, q_normal[-1]])
ax2.set_xlim([0.0, 0.07])
ax2.set_ylim([-0.6, 3])
ax1.set_ylabel(‘Frequency (THz)’)

# colorbar
norm = plt.Normalize(vmin=np.log10(alpha_vals.min()), vmax=np.log10(alpha_vals.max()))
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm._A = []
cbar = plt.colorbar(sm)
cbar.set_label(‘log10 (alpha)’)

fig.tight_layout()
fig.savefig(‘pdf/rotational_sum_rules.pdf’)
plt.show()
…
Sir, could you please help me to generate the script to import fitting with constraints ? I have SPOSCAR,POSCAR,FORCE_CONSTANTS files obtained from phonopy. Sir,I have not understood the working principle of this software fully. I am reading the theories regarding this software. Sir, kindly help me regarding this.

Best Regards,
Safikul

8

I suggest taking a look at the hiphive website and the basic tutorial and in order to understand the main objects/classes in hiphive.

In the script you’re running you get an error due to from hiphive.fitting import Optimizer it should be from trainstation import Optimizer.

1 Like
9

Dear Sir,

         The above **enforce_rotational_sum_rules_py script** contains the term **sc= StructureContainer.read ('data/structure_container.sc')**. This structure data has been generated by hiphive package. But I have structure files POSCAR, SPOSCAR and FORCE_CONSTANTS file from Phonopy. Sir, I have tried to import this files into the above script but getting many errors. Sir, my question is how should I import these files from Phonopy to the above **enforce_rotational_sum_rules_py script** to apply rotational sum rules. Also I getting the error regarding the package tools :

…
Traceback (most recent call last):
** File “/home/safikul_akd/Anaconda_physics/lib/python3.9/site-packages/hiphive/iit.py”, line 6, in **
** from tools import compute_dispersion**
ImportError: cannot import name ‘compute_dispersion’ from ‘tools’ (/home/safikul_akd/Anaconda_physics/envs/phonopy/lib/python3.10/site-packages/tools/init.py)
…
Sir, please help me regarding this.

Kind Regards,
Safikul

10

I think the formatting of your posts are bit incorrect, makes it hard to read.

If you only want to test with different values of alpha, then like I said this is just an additional parameter to enforce_rotational_sum_rules function, e.g. take your working script and add something like

alpha=0.5
parameters_rot = enforce_rotational_sum_rules(cs, parameters, [‘Huang’, ‘Born-Huang’], alpha=alpha)
11

Sorry Sir, I will try to rectify the formatting of my posts. As per your suggestion Sir I have tested the script to apply rotational sum rules with different values of alpha and plotted phonon dispersion curve for each value of alpha. But again the whole phonon spectra gets modified compared to the phonon spectra obtained from phonopy and also all the spectra are same for different alpha values. I have increased the cutoff value. Sir, Could you please explain the possible reasons behind it and how I should proceed now.

Kind Regards,
Safikul

1 Like
12

Dear Sir,

          Sorry for I am disturbing you again and again. I have been confused. Applying rotational sum rules change all the modes. Sir, please suggest what I should do now to remove only the negative frequencies of ZA branch near Gamma point and will get the correct phonon dispersion curve.

Kind Regards,
Safikul

13

Then you probably need to try large values of alpha, you should be able to produce a figure that looks like this i.e. where ZA modes is gradually shifted as you change alpha https://hiphive.materialsmodeling.org/_images/rotational_sum_rules.svg

There is no guarantee that this is possible to achieve.

The other modes changing is likely either due

  1. The cutoff imposed
    or
  2. Rotational sum rules being enforced

To test if the issue is related to the former (1.) I suggest you follow the recommendation on the hiphive-faq:

We recommend that you compare the dispersion between the original force constants, the ones produced by hiphive, and the ones produced by hiphive after enforcing the rotational sum rules.

14

Respected Sir,

                  I am trying to get phonon dispersion curve of SnS monolayer using hiphive. For the 1st step structure generation I am getting the following error:

Error:
raise NotImplementedError(‘No EMT-potential for {0}’
NotImplementedError: No EMT-potential for Sn

The script I am using is as follows:
Script:
from ase…io import read
from ase…io import write
from ase.calculators.emt import EMT
from hiphive.structure_generation import generate_mc_rattled_structures

# parameters
n_structures = 20
cell_size = 5
rattle_std = 0.02
minimum_distance = 2.3

# setup
prim = read(‘POSCAR’)
atoms_ideal = prim.repeat(cell_size)

# generate structures
structures = generate_mc_rattled_structures(atoms_ideal, n_structures, rattle_std,minimum_distance)

for atoms in structures:
** atoms.calc = EMT()**
** atoms.get_forces()**

# save structures
write(‘prim.extxyz’, prim)
write(‘supercell_ideal.extxyz’, atoms_ideal)
write(‘supercells_rattled.extxyz’, structures)
…
I am earnestly requesting you to please help me regarding this.

Kind Regards,
Safikul

15

Dear Sir,

         I have two more questions:
  1. The supercell I want to incorporate is 5x5x1. How will I do that? Using the above script this gives the position files of the supercell 5x5x5.
  2. After DFT calculations with this supercell structures generated I have to calculate force with the help of ASE. Suppose I have kept all the DFT output files to the folder DFT_OUTPUT. Then how should I progress?

Sir, please help me regarding this.

Kind Regards,
Safikul

16
  1. In ase you can use prim.repeat((5, 5, 1))
  2. Once you have done DFT calculations for the rattled supercells, then you need to read the DFT forces in your python script from DFT_OUTPUT. Depending on which DFT code you use, this can probably be done with the ase.io module File input and output — ASE documentation
17

Dear Sir,

           Thank you so much for helping me so far. I am grateful. Sir, I tried to generate DFT forces from the output files. I am using VASP and the output files are stored in the folder DFT_OUTPUT. Sir, in the DFT_OUTPUT folder there are 9 subfolders named (disp-00001..disp-00009) containing the vasp output files. Sir, could you please provide me the sample of the python script to generate DFT forces?

Kind Regards,
Safikul

18

Dear Sir,

        I am using the following script to add DFT forces but getting errors. I have to change the part **' forces=read('DFT_OUTPUT')** but unable to get it. Sir, please suggest what should I do to add DFT forces?

Sincerely,
Safikul
Script:
"""
Construct a ForceConstantPotential from training data generated previously.

Runs in approximately 100 seconds on an Intel Core i5-4670K CPU.
"""
from ase import ase
from ase.io import read,write
from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
from hiphive.utilities import prepare_structures
from trainstation import Optimizer

# read structures containing displacements and forces
forces=read(‘DFT_OUTPUT’)
prim = read(‘prim.extxyz’)
atoms_ideal = read(‘supercell_ideal.extxyz’)
rattled_structures = read(‘supercells_rattled.extxyz’, index=’:’)

# set up cluster space
cutoffs = [5.0, 4.0, 4.0]
cs = ClusterSpace(prim, cutoffs)
print(cs)
cs.print_orbits()

# … and structure container
structures = prepare_structures(rattled_structures, atoms_ideal)
sc = StructureContainer(cs)
for structure in structures:
** sc.add_structure(structure)**
print(sc)

# train model
opt = Optimizer(sc.get_fit_data())
opt.train()
print(opt)
# construct force constant potential
fcp = ForceConstantPotential(cs, opt.parameters)
fcp.write(‘fcc-nickel.fcp’)
print(fcp)

19

Respected Sir,

                  Using the following script I can generate the force in text format from each vasprun.xml files. Sir, I am facing one difficulty. How can I insert these forces in supercells_rattled.extxyz,prim.extxyz,supercell_ideal.extxyz files for further calculation? Sir, please help me regarding this.

Sincerely,
Safikul

Script:
import numpy as np
from ase.io import read,write
import ase.io
atoms=ase.io.read(filename=‘vasprun.xml0’,index=-1,format=“vasp-xml”)
atoms.get_forces()
force=atoms.get_forces()
print(force)
np.savetxt(‘force0’,force)

20

I can not write the python scripts for you. You need

  1. primitive structure
  2. ideal supercell (no displacements from ideal lattice)
  3. supercell with displacement and forces (vasprun.xml)

You can use the hiphive.utilities.prepare_structure function to setup an ASE atoms-object that you can directly pass to the StructureContainer.

supercell_ideal = read('POSCAR_supercell')
supercell_rattled = read('vasprun.xml')
supercell_training = prepare_structure(supercell_rattled, supercell_ideal)
...
sc = StructureContainer(cs)
sc.add_structure(supercell_training)