How to convert a second-order force constant matrix obtained by almode to FORCE_CONSTANTS in phonopy format

[YuwenCompPhys]

Now that I have the si.fcs and si.xml files, how do I convert them to FORCE_CONSTANTS

[ttadano]

You can use the ALM option in phonopy.
https://phonopy.github.io/phonopy/command-options.html

[YuwenCompPhys]

Thank you very much for your answer. I mean, now that I have obtained the renormalized force constant file xxx.xml at a finite temperature, I now want to convert the renormalized force constant file directly to FORCE_CONSTANS, can you give me a case of using the operation? Thank you!

[LiuHaoXuan0403]

Hi Zhangyuwen, I have the the same question. Could you tell how to get renormalized force constant file XXX.fcs from XXX.xml.

[ZLS1110]

I have the same question "How to covert the renormalized force constant file xxx.xml at a finite temperature to the force constant file with the format (FORCE_CONSTANS)?" I have seen "scph_to_qefc.py" in the tools which can add anharmonic correction to a QE force constant format. I want to ask if there is a similar code about coverting the file into the format (FORCE_CONSTANS) ?

[ttadano]

Here is an example script for converting an alamode *.xml file into FORCE_CONSTANTS format.

import sys
from lxml import etree
from ase.units import Rydberg, Bohr
import numpy as np

def parse_xml(fname_xml):
    try:
        # Attempt to parse the XML file directly with lxml
        tree = etree.parse(fname_xml)
        return tree
    except etree.XMLSyntaxError:
        # If a parsing error occurs with lxml,
        # attempt to repair by re-reading the file with the 'recover' parser
        repair_parser = etree.XMLParser(recover=True)
        tree = etree.parse(fname_xml, parser=repair_parser)
        return tree

def get_forceconstants_xml(fname_xml):

    xml = parse_xml(fname_xml)
    root = xml.getroot()

    natom_super = int(root.find('Structure/NumberOfAtoms').text)
    ntrans = int(root.find('Symmetry/NumberOfTranslations').text)
    natom_prim = natom_super // ntrans

    # parse mapping table
    map_p2s = np.zeros((ntrans, natom_prim), dtype=int)
    for elems in root.findall('Symmetry/Translations/map'):
        itran = int(elems.get('tran')) - 1
        iatom = int(elems.get('atom')) - 1
        map_p2s[itran, iatom] = int(elems.text) - 1


    # parse harmonic force constants
    fc2_compact = np.zeros((natom_prim, natom_super, 3, 3), dtype=float)

    for elems in root.findall('ForceConstants/HARMONIC/FC2'):
        atom1, xyz1 = [int(t)-1 for t in elems.get('pair1').split()]
        atom2, xyz2, icell2 = [int(t)-1 for t in elems.get('pair2').split()]
        fcsval = float(elems.text)

        fc2_compact[atom1, atom2, xyz1, xyz2] += fcsval

    fc2_compact *= Rydberg / (Bohr**2)

    return map_p2s, fc2_compact

def print_fc2_phonopy(map_p2s, fc2_compact):

    natom_prim, natom_super = fc2_compact.shape[:2]

    print("{:5d} {:5d}".format(natom_prim, natom_super))
    for i in range(natom_prim):
        for j in range(natom_super):
            print("{:5d} {:5d}".format(map_p2s[0,i]+1, j+1))
            for k in range(3):
                for l in range(3):
                    print("{:20.15f}".format(fc2_compact[i, j, k, l]), end='')
                print('')

def main():
    fname_xml = sys.argv[1]
    map_p2s, fc2_compact = get_forceconstants_xml(fname_xml)
    print_fc2_phonopy(map_p2s, fc2_compact)


if __name__ == '__main__':
    main()

Please save the above as a python file, for example convert.py, and run

python convert.py alamode.xml > FORCE_CONSTANTS

Please test if this works properly for your specific material. If not, please fix the issue.

[liudadie12138]

@ttadano Thank you for taking the time to answer my question, I mean the third-order IFC that needs to be re-normalized, because I want to calculate the thermal conductivity after reforming, but I checked the information and found that you replied to other authors that alamode currently only supports the second-order force constant after reforming + the old third-order force constant for thermal conductivity calculation, I didn't find a specific tutorial for this calculation method, I asked someone and he told me that it should be modified in the RTA process
&general
PREFIX=pc rta
MODE=RTA
NKD=4
FCSXML=out3.xml
FC2XML=out2 800.xm
TMIN=800
TMAX=800
DT=50
Here FCSXML uses the fourth-order force constant of opt fitting after cv, and FC2XML is the second-order force constant of dfc2 reforming, I will confirm whether this process is correct, and its result is the thermal conductivity after reforming. In addition, I mentioned in your answer to other people's questions that dfc2 needs to specify a temperature, and this temperature you pointed out should be equal to Tscph, does that mean that I should use TMIN=TMAX when calculating Tscph, or can I set the range as I like? Finally, thank you for answering my questions.

[ttadano]

As you mentioned, we do not currently have a handy script to compute/output renormalized third-order IFCs.

does that mean that I should use TMIN=TMAX when calculating Tscph,

When the FC2XML is an scph-corrected IFC2 at 800 K, TMIN and TMAX should be set to 800 as you did.

[liudadie12138]

@ttadano
Thanks for your answer, so I can think that the above process for calculating the thermal conductivity after reforming is correct,
Also, I have one last question, I should do opt after fitting the fourth-order force constant, cv, this is my opt file
&general
PREFIX = ct221
MODE = optimize
NAT = 160; NKD = 3
KD = Na Nb O
/

&interaction
NORDER = 3 # 1: harmonic, 2: cubic, ..
NBODY = 2 3 3
/

&cutoff
- None 10.0 9.0
/

&optimize
LMODEL = ols
NDATA = 40
DFSET = DFSET_AIMD_random
FC2XML = NaNbO221.xml

CV = 0 # switch off CV
L1_RATIO = 1.0 # LASSO
L1_ALPHA = 4.64159e-05 # Add this line.
CV_MINALPHA = 1.0e-8
CV_MAXALPHA = 0.01
CV_NALPHA = 30

STANDARDIZE = 1
CONV_TOL = 1.0e-9 # was CONV_TOL = 1.0e-8 in CV

/

&cell
1.889726
10.9833793640 0.0000000000 0.0000000000
0.0000000000 11.1092567444 0.0000000000
0.0000000000 0.0000000000 15.4741468430
/
My previous LMODEL=enet fit residuals were 30%, which is not ideal
I tried to modify LMODEL = ols, it is still being calculated, if this fitting effect is not good, which parameters can I adjust to improve the fitting effect,
Then I would like to ask how to determine the cutoff here, the recommended value you gave in the tutorial is none 15 9, is it possible that this value can be explicitly obtained from the log like the second-order force constant, or is it an empirical value that needs to be tested?

[ttadano]

Can you make a separate discussion topic for a new question?

[vivekgusain09]

Dear all,
I am writing this query regarding the strange results i am getting when i am using the second order force constants (harmonic as well as temperature normalized) obtained using alamode. When computing the phonon dispersion in alamode and in phonopy seperately it gives the correct results and there are no imaginary frequencies. However, the moment I am using the second order force constants obtained using alamode to computed phonon dispersions and lattice thermal conductivty in phonopy and ShengBTE respectively, i am getting wrong phonon dispersion with imaginary frequencies and unphysically low values of lattice thermal conductivity.

The procedure i am following is as follows:

obtained the second order force constants using alamode in .xml format.
Converted the .xml file into shengbte/phonopy format using the convert_fc2.py script that has been provided with the alamode in the tools directory.
Command used to extract the force constants: python3 convert_fc2.py -i s222.xml -o FORCE_CONSTANTS --full

After comparing the converted force constants with the force constants obtained using phonopy, i can see the values are different for the same pair of atoms. I computed the lattice thermal conductivity and phonons with phonopy obtained second order force constants and they are giving the correct and physically meaningfull results, similar to the results for lattice thermal conductivity obtained

using the alamode code. I think there is some problem with the conversion which i am not able to identify.

I kindly request for your help and guidance in this regard as it is extremely important for my research. Any response is highly appreciated.
Thank you for your time and consideration.

I have uploaded the phonon dispersions and force constant files here for the reference

converted_FORCE_CONSTANTS.txt
phonopy_FORCE_CONSTANTS.txt

[sabofavian829-eng]

@LiuHaoXuan0403 Hello friend, I have the same question as you regarding the renormalized second-order force constant. I now have correction values at different temperatures. When calculating the thermal conductivity at different temperatures, can I fix a corrected second-order moment (force constant) to compute the thermal conductivity?Thank you very much, and I hope you can answer this question.