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

[zhangyuwen984]

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

[zhangyuwen984]

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.

[Lamed-git]

FULL_FORCE_CONSTANTS = .TRUE. WRITE_FORCE_CONSTANTS = .TRUE. Note: the old FORCE_CONSTANTS will be overwrited
Please refer to gangliu-github, use the flags in phonopy toconvert to FULL_FORCE_CONSTANTS.

[sc696]

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.

This code works perfectly, however, this code generates the FORCE_CONSTANT in compact form, which cannot be used in Shengbte. In order to use in Shengbte one need the FORCE_CONSTANT in full form. Can you kindly provide a script which can convert the alamode.xml file into phonopy FORCE_CONSTANT in full form.

[ttadano]

Hi,

I've added convert_fc2.py in the tools directory.
This is an improved version of the script shown in this thread above.

To convert the ALAMODE XML FC2 into FORCE_CONSTANTS of phonopy in the full form, please issue

python convert_fc2.py -i alamode_fc.xml -o FORCE_CONSTANTS -full

[sc696]

Hi,

Thanks for your help. However, I did not find the file "convert_fc2.py" in the tools directory of the most recent version of v.1.5.0. Can you kindly help me sorting this problem.

[ttadano]

The script is in the develop branch at https://github.com/ttadano/alamode/blob/develop/tools/convert_fc2.py.