How to start with thermal conductivity calculation

[biglinn]

I find that calculating thermal conductivity needs ML potential to get starting dynamical matrices, as well as relaxing structure in the 8th tutorial.
Since I have got dynamical matrices by SSCHA before, so could I use these as strating matrices instead of calculating by ML potential?
My script is like this:

import cellconstructor as CC, cellconstructor.Phonons
import sscha, sscha.Ensemble
import numpy as np

# Fix the seed so that we all generate the same ensemble
np.random.seed(0)

# Load the dynamical matrix
dyn = CC.Phonons.Phonons("sscha_dyn_population10_", 4)
#dyn.ForcePositiveDefinite()
dyn.Symmetrize()

ensemble = sscha.Ensemble.Ensemble(dyn,300)
ensemble.generate(20)

dyn_hessian, d3_tensor = ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')

(Little question: What difference does hessian_dyn have in terms of my starting matrices?)

After getting hessian_dyn and third order FC(d3.npy), do subsequent calculation:

from __future__ import print_function
from __future__ import division

import numpy as np
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.ForceTensor
import cellconstructor.ThermalConductivity
import time

dyn_prefix = 'hessian_dyn'
nqirr = 4

SSCHA_TO_MS = cellconstructor.ThermalConductivity.SSCHA_TO_MS
RY_TO_THZ = cellconstructor.ThermalConductivity.SSCHA_TO_THZ
dyn = CC.Phonons.Phonons(dyn_prefix, nqirr)

supercell = dyn.GetSupercell()

fc3 = CC.ForceTensor.Tensor3(dyn.structure,
dyn.structure.generate_supercell(supercell), supercell)

d3 = np.load("d3.npy")
fc3.SetupFromTensor(d3)
fc3 = CC.ThermalConductivity.centering_fc3(fc3)

mesh = [2,2,2]
smear = 0.03/RY_TO_THZ

tc = CC.ThermalConductivity.ThermalConductivity(dyn, fc3,
kpoint_grid = mesh, scattering_grid = mesh, smearing_scale = None,
smearing_type = 'constant', cp_mode = 'quantum', off_diag = False)

temperatures = np.linspace(200,1200,10,dtype=float)
start_time = time.time()
tc.setup_harmonic_properties(smear)
tc.write_harmonic_properties_to_file()

tc.calculate_kappa(mode = 'SRTA', temperatures = temperatures,
write_lifetimes = True, gauss_smearing = True, offdiag_mode = 'wigner',
kappa_filename = 'Thermal_conductivity_SRTA', lf_method = 'fortran-LA')

print('Calculated SSCHA kappa in: ', time.time() - start_time)

tc.calculate_kappa(mode = 'GK', write_lineshapes=False,
ne = 1000, temperatures = temperatures,
kappa_filename = 'Thermal_conductivity_GK')

print('Calculated SSCHA kappa in: ', time.time() - start_time)
# Save ThermalConductivity object for postprocessing.
tc.save_pickle()

However, there is not corrcet result:

[mesonepigreco]

@DjordjeDangic maybe can help

Regarding the dynamical matrix, you should not use the Hessian for the Thermal Conductivity calculation but rather the final SSCHA matrix.
In this workflow, the Hessian is only needed to compute the third-order force constant tensor (the d3 in the script).

If the Hessian has imaginary frequencies, that could be the cause of the NaN in the result.

[biglinn]

Thanks! I got it.
But the script in the tutorial has some steps, first generating a set of ensembles, next relaxing structure, then generating new ensembles and saving dynamical matrices.
Is it a standrad workflow? Could I directly use my calculated matrices as the final SSCHA matrices?

[biglinn]

Or is the step of relaxing structure necessary?
Because I find that the script's main function is to get bulk modulus for doing relaxing by vc_relax.

[mesonepigreco]

Hi, what tutorial are you referring to?
For the thermal conductivity, you can directly use the final SSCHA dynamical matrices and the third-order force constant, no need to do anything else.

[biglinn]

Hi, it's the 8th tutorial(Thermal conductivity) on the website. The first script is following:

import numpy as np
from quippy.potential import Potential
from ase import Atoms
import ase.io
from ase.eos import calculate_eos
from ase.units import kJ
from ase.phonons import Phonons as AsePhonons
from ase.constraints import ExpCellFilter
from ase.optimize import BFGS, QuasiNewton
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.Structure
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax


# This function will use ASE to give us a starting
# guess for dynamical matrices
def get_starting_dynamical_matrices(structure_filename,
        potential, supercell):
    atoms = ase.io.read(structure_filename)
    atoms.set_calculator(potential)
    ecf = ExpCellFilter(atoms)
    qn = QuasiNewton(ecf)
    qn.run(fmax=0.0005)

    structure = CC.Structure.Structure()
    structure.generate_from_ase_atoms(atoms, get_masses = True)
    dyn = CC.Phonons.compute_phonons_finite_displacements(structure,
    potential, supercell = supercell)
    dyn.Symmetrize()
    dyn.ForcePositiveDefinite()

    eos = calculate_eos(atoms)
    v0, e0, B = eos.fit()
    bulk = B / kJ * 1.0e24

    return dyn, bulk


# Our input variables
temperature = 100.0
nconf = 1000
max_pop = 1000

# Load in Tersoff potential
pot = Potential('IP Tersoff', param_filename=
'../06_the_SSCHA_with_machine_learning_potentials/ip.parms.Tersoff.xml')
supercell = tuple((4*np.ones(3, dtype=int)).tolist())
dyn, bulk = get_starting_dynamical_matrices(
'../06_the_SSCHA_with_machine_learning_potentials/POSCAR', pot, supercell)


# Generate the ensemble and the minimizer objects
ensemble = sscha.Ensemble.Ensemble(dyn, T0=temperature,
supercell = dyn.GetSupercell())
ensemble.generate(N = nconf)
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minimizer.min_step_dyn = 0.1
minimizer.kong_liu_ratio = 0.5
minimizer.meaningful_factor = 0.001
minimizer.max_ka = 1000

# Relax structure
relax = sscha.Relax.SSCHA(minimizer, ase_calculator = pot,
N_configs = nconf, max_pop = max_pop, save_ensemble = True)
relax.vc_relax(static_bulk_modulus = bulk,
ensemble_loc = "directory_of_the_ensemble")

# Generate ensemble for third-order FC with the relaxed dynamical matrices
new_ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0=temperature,
supercell = relax.minim.dyn.GetSupercell())
new_ensemble.generate(N = nconf*5)
new_ensemble.compute_ensemble(pot, compute_stress = True,
stress_numerical = False, cluster = None, verbose = True)
# We minimize the free energy with this new ensemble
new_minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
new_minimizer.minim_struct = False
new_minimizer.set_minimization_step(0.1)
new_minimizer.meaningful_factor = 0.001
new_minimizer.max_ka = 10000
new_minimizer.init()
new_minimizer.run()
new_minimizer.dyn.save_qe('final_dyn')
# Update weights with a new dynamical matrice
new_ensemble.update_weights(new_minimizer.dyn, temperature)

# Calculate Hessian and the third order tensor (return_d3 = True)
dyn_hessian, d3_tensor = new_ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')
Here we used 4x4x4 supercell. You will need to converge results with respect to the size of the supercell. A good check for the convergence could be the decay of the second and third order force constants with the distance. Now that we have second and third order force constants, we can calculate lattice thermal conductivity. For this we provide following script:

from __future__ import print_function
from __future__ import division

import numpy as np
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.ForceTensor
import cellconstructor.ThermalConductivity
import time

dyn_prefix = 'final_dyn'
nqirr = 8

SSCHA_TO_MS = cellconstructor.ThermalConductivity.SSCHA_TO_MS
RY_TO_THZ = cellconstructor.ThermalConductivity.SSCHA_TO_THZ
dyn = CC.Phonons.Phonons(dyn_prefix, nqirr)

supercell = dyn.GetSupercell()

fc3 = CC.ForceTensor.Tensor3(dyn.structure,
dyn.structure.generate_supercell(supercell), supercell)

d3 = np.load("d3.npy")
fc3.SetupFromTensor(d3)
fc3 = CC.ThermalConductivity.centering_fc3(fc3)

mesh = [10,10,10]
smear = 0.03/RY_TO_THZ

tc = CC.ThermalConductivity.ThermalConductivity(dyn, fc3,
kpoint_grid = mesh, scattering_grid = mesh, smearing_scale = None,
smearing_type = 'constant', cp_mode = 'quantum', off_diag = False)

temperatures = np.linspace(200,1200,10,dtype=float)
start_time = time.time()
tc.setup_harmonic_properties(smear)
tc.write_harmonic_properties_to_file()

tc.calculate_kappa(mode = 'SRTA', temperatures = temperatures,
write_lifetimes = True, gauss_smearing = True, offdiag_mode = 'wigner',
kappa_filename = 'Thermal_conductivity_SRTA', lf_method = 'fortran-LA')

print('Calculated SSCHA kappa in: ', time.time() - start_time)

tc.calculate_kappa(mode = 'GK', write_lineshapes=False,
ne = 1000, temperatures = temperatures,
kappa_filename = 'Thermal_conductivity_GK')

print('Calculated SSCHA kappa in: ', time.time() - start_time)
# Save ThermalConductivity object for postprocessing.
tc.save_pickle()

Do we have to get final matrices like this?

[biglinn]

You say that we can directly use final SSCHA matrices. Does the matrices refer to those by free energy minimization?

[DjordjeDangic]

As it is written in the tutorial you need only the final auxiliary SSCHA force constants ('final_dyn') and the third-order force constants ('d3.npy'). Do not use Hessian force constants ('hessian_dyn').

In the tutorial, we are relaxing structure prior to the calculation of the thermal conductivity because we want to calculate thermal conductivity at 0 pressure and temperature at 100 K. Later, we assume that structure and force constants do not significantly change with temperature so we are using the same force constants to calculate thermal conductivity at elevated temperatures (higher than 100 K).

Let's say you want to calculate lattice thermal conductivity of PbTe at 300 K at ambient pressure. Then you will first need to do SSCHA relaxation of PbTe at 300 K and 0 GPa. Once your relaxation finishes you will get dynamical matrices ('final_dyn' in tutorial). Then you need to calculate third order force constants ('d3.npy' in tutorial). Using these two ingredients you can calculate lattice thermal conductivity. Whether you do this using machine learning potential or DFT does not matter. However, if you are using machine learning potential you need to make sure it is faithful representation of interatomic interaction.

[biglinn]

Thanks! I see what you mean.
But I did a test. I don't realx structure, like this:

import cellconstructor as CC, cellconstructor.Phonons
import sscha, sscha.Ensemble
import numpy as np

# Fix the seed so that we all generate the same ensemble
np.random.seed(0)

# Load the dynamical matrix
dyn = CC.Phonons.Phonons("sscha_dyn_population5_", 3)
#dyn.ForcePositiveDefinite()
dyn.Symmetrize()

ensemble = sscha.Ensemble.Ensemble(dyn,300)
ensemble.generate(20)

dyn_hessian, d3_tensor = ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')

I just read the matrices sscha_dyn_population5_ to produce d3.npy, then the result of thermal conductivity Thermal_conductivity_SRTA is following:

#  Temperature (K)   Kappa xx (W/mK)   Kappa yy (W/mK)   Kappa zz (W/mK)   Kappa xy (W/mK)   Kappa yz (W/mK)   Kappa zx (W/mK)
   2.000000000000e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   3.111111111111e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   4.222222222222e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   5.333333333333e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   6.444444444444e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   7.555555555556e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   8.666666666667e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   9.777777777778e+02   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   1.088888888889e+03   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00
   1.200000000000e+03   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00   0.000000000000e+00

According to my understanding, relaxing structure is to calculate thermal conductivity at particular pressue and temperature. So not doing relaxing should also correspond to a certain pressure and temperature.
So, this way should also be right? But my results of thermal conductivity are all zero.

[DjordjeDangic]

And the code did not print out any warning or error? Can you check your phonon frequencies, group velocities and phonon lifetimes? That would tell us where the issue is.

Also you are calculating third-order force constants with 20 configurations according to your script. That is not nearly enough.

[biglinn]

I don't receive prompts of warning or erro.
Lifetime is really strange, as follows:

All Phonon_transport_properties_XX files are like this. Lifetime is all zeros.

And I try to increase number of configurations to 100, 1000. The results don't change.

[biglinn]

So, I wonder if it is necessary to relax structure with a potential (such as ML potential).

If so, does calculating thermal conductivity have nothing to do with doing SSCHA simulations(free energy minimization)?
Because I find that the script starts with a POSCAR, instead of using calculated dynamical matrices.

[biglinn]

Hi, @DjordjeDangic
I take two calculations with relaxing structure using ML potential.
One is according to the script of tutorial to calculate bulk modulus, relax.vc_relax(static_bulk_modulus = bulk, ensemble_loc = "data"). The complete script is as follows:

import numpy as np
from deepmd.calculator import DP
from ase import Atoms
import ase.io
from ase.eos import calculate_eos
from ase.units import kJ
from ase.constraints import ExpCellFilter
from ase.optimize import BFGS, QuasiNewton
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.Structure, cellconstructor.calculators
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax

def get_starting_dynamical_matrices(structure_filename,
        potential, supercell):
    atoms = ase.io.read(structure_filename)
    atoms.set_calculator(potential)
    ecf = ExpCellFilter(atoms)
    qn = QuasiNewton(ecf)
    qn.run(fmax=0.0005)

    structure = CC.Structure.Structure()
    structure.generate_from_ase_atoms(atoms, get_masses = True)
    dyn = CC.Phonons.compute_phonons_finite_displacements(structure,
    potential, supercell = supercell)
    dyn.Symmetrize()
    dyn.ForcePositiveDefinite()

    eos = calculate_eos(atoms)
    v0, e0, B = eos.fit()
    bulk = B / kJ * 1.0e24

    return dyn, bulk

model_path ='./iter15_000.pb'
calculator =DP(model=model_path)
supercell = (2,2,2)
dyn, bulk = get_starting_dynamical_matrices('POSCAR', calculator, supercell)

# Our input variables
temperature = 300.0
nconf = 100
max_pop = 1000

# Generate the ensemble and the minimizer objects
ensemble = sscha.Ensemble.Ensemble(dyn, T0=temperature)
ensemble.generate(N = nconf)
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minimizer.min_step_dyn = 0.1
minimizer.kong_liu_ratio = 0.5
minimizer.meaningful_factor = 0.001
minimizer.max_ka = 1000

# Relax structure
relax = sscha.Relax.SSCHA(minimizer, calculator, N_configs = nconf, max_pop = max_pop, save_ensemble = True)
#relax.vc_relax(target_press = 0)
relax.vc_relax(static_bulk_modulus = bulk, ensemble_loc = "data")

# Generate ensemble for third-order FC with the relaxed dynamical matrices
new_ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0=temperature)
new_ensemble.generate(N = nconf*5)
new_ensemble.compute_ensemble(calculator, compute_stress = True,stress_numerical = False, cluster = None, verbose = True)
# We minimize the free energy with this new ensemble
new_minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
new_minimizer.minim_struct = False
new_minimizer.set_minimization_step(0.1)
new_minimizer.meaningful_factor = 0.001
new_minimizer.max_ka = 10000
new_minimizer.init()
new_minimizer.run()
new_minimizer.dyn.save_qe('final_dyn')
# Update weights with a new dynamical matrice
new_ensemble.update_weights(new_minimizer.dyn, temperature)

# Calculate Hessian and the third order tensor (return_d3 = True)
dyn_hessian, d3_tensor = new_ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')
import numpy as np
from deepmd.calculator import DP
from ase import Atoms
import ase.io
from ase.eos import calculate_eos
from ase.units import kJ
from ase.constraints import ExpCellFilter
from ase.optimize import BFGS, QuasiNewton
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.Structure, cellconstructor.calculators
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax

def get_starting_dynamical_matrices(structure_filename,
        potential, supercell):
    atoms = ase.io.read(structure_filename)
    atoms.set_calculator(potential)
    ecf = ExpCellFilter(atoms)
    qn = QuasiNewton(ecf)
    qn.run(fmax=0.0005)

    structure = CC.Structure.Structure()
    structure.generate_from_ase_atoms(atoms, get_masses = True)
    dyn = CC.Phonons.compute_phonons_finite_displacements(structure,
    potential, supercell = supercell)
    dyn.Symmetrize()
    dyn.ForcePositiveDefinite()

    eos = calculate_eos(atoms)
    v0, e0, B = eos.fit()
    bulk = B / kJ * 1.0e24

    return dyn, bulk

model_path ='./iter15_000.pb'
calculator =DP(model=model_path)
supercell = (2,2,2)
dyn, bulk = get_starting_dynamical_matrices('POSCAR', calculator, supercell)

# Our input variables
temperature = 300.0
nconf = 100
max_pop = 1000

# Generate the ensemble and the minimizer objects
ensemble = sscha.Ensemble.Ensemble(dyn, T0=temperature)
ensemble.generate(N = nconf)
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minimizer.min_step_dyn = 0.1
minimizer.kong_liu_ratio = 0.5
minimizer.meaningful_factor = 0.001
minimizer.max_ka = 1000

# Relax structure
relax = sscha.Relax.SSCHA(minimizer, calculator, N_configs = nconf, max_pop = max_pop, save_ensemble = True)
#relax.vc_relax(target_press = 0)
relax.vc_relax(static_bulk_modulus = bulk, ensemble_loc = "data")

# Generate ensemble for third-order FC with the relaxed dynamical matrices
new_ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0=temperature)
new_ensemble.generate(N = nconf*5)
new_ensemble.compute_ensemble(calculator, compute_stress = True,stress_numerical = False, cluster = None, verbose = True)
# We minimize the free energy with this new ensemble
new_minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
new_minimizer.minim_struct = False
new_minimizer.set_minimization_step(0.1)
new_minimizer.meaningful_factor = 0.001
new_minimizer.max_ka = 10000
new_minimizer.init()
new_minimizer.run()
new_minimizer.dyn.save_qe('final_dyn')
# Update weights with a new dynamical matrice
new_ensemble.update_weights(new_minimizer.dyn, temperature)

# Calculate Hessian and the third order tensor (return_d3 = True)
dyn_hessian, d3_tensor = new_ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')

The other is just doing general vc_relax, relax.vc_relax(target_press = 0, ensemble_loc = "data"). The complete script is as follows:
import numpy as np

from deepmd.calculator import DP
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.Structure, cellconstructor.calculators
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax, sscha.Utilities

model_path ='./iter15_000.pb'
calculator =DP(model=model_path)

the_structure = CC.Structure.Structure()
the_structure.read_generic_file("POSCAR")
supercell = (2,2,2)
dyn = CC.Phonons.compute_phonons_finite_displacements(the_structure, calculator, supercell = supercell)
dyn.Symmetrize()
dyn.ForcePositiveDefinite()

# Our input variables
temperature = 300.0
nconf = 100
max_pop = 1000

# Generate the ensemble and the minimizer objects
ensemble = sscha.Ensemble.Ensemble(dyn, T0=temperature)
ensemble.generate(N = nconf)
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minimizer.min_step_dyn = 0.1
minimizer.kong_liu_ratio = 0.5
minimizer.meaningful_factor = 0.001
minimizer.max_ka = 1000

# Relax structure
relax = sscha.Relax.SSCHA(minimizer, calculator, N_configs = nconf, max_pop = max_pop, save_ensemble = True)
#relax.vc_relax(target_press = 0)
relax.vc_relax(target_press = 0, ensemble_loc = "data")

# Generate ensemble for third-order FC with the relaxed dynamical matrices
new_ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0=temperature)
new_ensemble.generate(N = nconf*5)
new_ensemble.compute_ensemble(calculator, compute_stress = True,stress_numerical = False, cluster = None, verbose = True)
# We minimize the free energy with this new ensemble
new_minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
new_minimizer.minim_struct = False
new_minimizer.set_minimization_step(0.1)
new_minimizer.meaningful_factor = 0.001
new_minimizer.max_ka = 10000
new_minimizer.init()
new_minimizer.run()
new_minimizer.dyn.save_qe('final_dyn')
# Update weights with a new dynamical matrice
new_ensemble.update_weights(new_minimizer.dyn, temperature)

# Calculate Hessian and the third order tensor (return_d3 = True)
dyn_hessian, d3_tensor = new_ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')

The second calculation is successfully done.
But the first is not. The final_dynX have some small imaginary frequency, like this:

If calculating thermal conductivity next, there are erros:

/public/home/biglinn/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/ForceTensor.py:255: ComplexWarning: Casting complex values to real discards the imaginary part
  self.tensor[i_block, 3*na1:3*na1+3, 3*na2: 3*na2+3] = tensor[3*na1 : 3*na1+3, 3*nat2_sc : 3*nat2_sc + 3]
/public/home/biglinn/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/ThermalConductivity.py:2646: ComplexWarning: Casting complex values to real discards the imaginary part
  self.gvels[ikpt], self.ddms[ikpt] = self.get_group_velocity(kpt, self.freqs[ikpt], self.eigvecs[ikpt])
/public/home/biglinn/miniconda3/envs/sscha/lib/python3.10/site-packages/cellconstructor/ThermalConductivity.py:3175: RuntimeWarning: invalid value encountered in sqrt
  w_q = np.sqrt(w2_q)

So why does this situation appear?

[mesonepigreco]

If you have julia installed in your system, and the last version of cellconstructor and python-sscha, then it should work automatically.
If not, the code throws a warning at the beginning, saying that Julia is not available and that it falls back on the old implementation.
If you followed the installation instructions on the website, you should have already installed also julia in the conda sscha environment, and it should work out-of-the-box.

To run in parallel, you can submit it with

mpirun -np NPROC python myscript.py

The mpi version must coincide with the one used for mpi4py, so check if the code prints the same output multiple times (in that case, it is not running in parallel, and you should try a different mpi version).
To switch between MPI versions, replace the mpirun with the specific command, such as mpirun.mpich or mpirun.openmpi or something similar, depending on your system.

[biglinn]

About the last version, what I'm using is 1.4.1 cellconstructor and 1.4.1 python-sscha with 0.6.2 julia installed. Could this version work in parallel?
And I work in a SLURM cluster, so can I submit job with srun -n NPROC python myscript.py?

[biglinn]

As you said, the code prints the same output muitiple times, like this:

Also some prompts:

The application appears to have been direct launched using "srun",
but OMPI was not built with SLURM's PMI support and therefore cannot
execute. There are several options for building PMI support under
SLURM, depending upon the SLURM version you are using:

  version 16.05 or later: you can use SLURM's PMIx support. This
  requires that you configure and build SLURM --with-pmix.

  Versions earlier than 16.05: you must use either SLURM's PMI-1 or
  PMI-2 support. SLURM builds PMI-1 by default, or you can manually
  install PMI-2. You must then build Open MPI using --with-pmi pointing
  to the SLURM PMI library location.

Please configure as appropriate and try again.

Does it need to import some parallel modules in myscrpit.py?

[mesonepigreco]

As said, the mpi4py was built with a different version of MPI than the one used by SLURM.
This is normal if the code was installed with Anaconda.
In this case, you should use the mpi version of Anaconda. So, instead of running Python with srun, you should use mpirun.
However, this means that MPI will be much less optimized on the cluster. To fix this, you can ask the cluster managers how to have mpi4py work with your Anaconda environment where SSCHA is installed.

[biglinn]

Oh, that sounds like a complicated problem. Thanks! I'll take a try.

[biglinn]

Hi, @DjordjeDangic

If the thermal conductivity is two small(about an order of magnitude smaller than experiment), what's the possible reason?
In the tutorial, there are several parameters mentioned for convergence, smear, mesh etc. For calculating appropriate thermal conductivity, are there any requirements for these parameters?

[biglinn]

About the number of configurations in minimization, I have tested 10000 configurations from the first 500. But the thermal conductivity is still not converged, keeping increasing.
I feel strange. Does it need more configurations? Or is there somethhing else at work?

[DjordjeDangic]

It's impossible to say without looking at the data. However, 10000 configurations is not unheard of. Are results converging, is the difference between 10k and 9k smaller than between 9k and 8k?

Check if you are looking at the primitive cell of SnSe!

[biglinn]

Yeah, I'm looking at the primitive cell, like this:

I have test 2000, 4000, 6000, 10000 and 20000. They are respectively as follows:
The top row from left to right is 2000, 4000 and 6000, and The bottom row is 10000 and 20000.

The difference between 10k and 6k is even smaller that that between 6k and 4k.
But until 20k, the valuse still increase a little. Should I use more configurations?

[DjordjeDangic]

I think it should be converged, since I have performed several calculations and they are all like the figure above. By the way, is there some parameter to judge if results are converged? The phonon dispersion is as follows:

And following is one of Phonon_transport_properties_X, the first three lines are zeros.

How many phonon bands do you see in your graph? Because I definitely see more than 12! So you can not have only 4 atoms per primitive cell.

This is not the way to present a convergence study. How are we supposed to compare different calculations when they are on five separate graphs? But yeah, I would agree with you that results do not seem converged if I am assuming correctly labels for each graph. Also, are you using the same smearing and q-point grids for each of these calculations?

[biglinn]

I think it should be converged, since I have performed several calculations and they are all like the figure above. By the way, is there some parameter to judge if results are converged? The phonon dispersion is as follows:

And following is one of Phonon_transport_properties_X, the first three lines are zeros.

I used the conventional cell before, so there are more than 12 phonon bands. Now, I have changed to the primitive cell. Its phonon dispersion is like this:

Pardon me, I should have put them on one graph. But I can be sure that labels are correct.
And I keep smearing and q-grid the same, since I only change the number of configurations in different calculations. That's why I think more configurations may be needed.