The ABINIT Users Mailing List ( CLOSED )

[verstraete@pcpm.ucl.ac.be]

Yes.

In the next version of abinit, a python script of mine will be included,
which looks up the irreps for a given space group, and uses them to
project the characters obtained in anaddb, giving the irreps for each
mode. It does not yet give Raman/IR, and you have to consider the LO TO
splitting yourself, as in anaddb it mixes irreps (a priori) and has to be
disentangled. This means 1) get the irreps for the unsplit modes, and 2)
see where the mode goes, to get the irreps for the LO/TO ones.

The script is attached with an additional python module. Both must be
present, and let me stress I offer no support for the moment. Please tell
me if it works or crashes, but the version is not official yet, and I
don't have time to tune it, test it, or debug it.

Matthieu

-- 
===================================================================
Matthieu Verstraete               mailto:verstraete@pcpm.ucl.ac.be 
PCPM, Boltzmann, pl. Croix du Sud, 1            tel: 010/ 47 86 81
B-1348 Louvain-la-Neuve  Belgium                fax: 010/ 47 34 52
#! /usr/bin/python

#
#    Copyright (C) 2003 ABINIT group
#
#    Written by Matthieu Verstraete in python (compatible v1.9).
#    This is free software, and you are welcome to redistribute it
#    under certain conditions (GNU General Public License,
#    see ~ABINIT/Infos/copyright or http://www.gnu.org/copyleft/gpl.txt).
#
#    ABINIT is a project of the Universite Catholique de Louvain,
#    Corning Inc. and other collaborators, see ~ABINIT/Infos/contributors.
#    Please read ~ABINIT/Infos/acknowledgments.htm for suggested
#    acknowledgments of the ABINIT effort.
#
#    For more information, see http://www.abinit.org .
#
#    Input is files from ABINIT: ground state or phonon runs, and from anaddb
#    for the Gamma point at least.
#
#    TODO:
#    * LO/TO splitting is not taken into account yet for the calculation
#      of the characters in anaddb so the characters are for the TO modes
#      only. The splitting can be identified afterwards.
#    * Add the primes and plus/minus to the irrep symbols according to the
#      2-axis and mirror symmetry character of each irrep.
#    * If the symmetries are one day output to the anaddb files, only
#      one input file will be needed (ie parse it for syms, instead of 
#      "abifilename"
#


#  reg expss
import re
import string
import sys
#import fileinput
import math
import xml.dom.minidom
from tokenize_file import *
import urllib


elemabbrev = 
['XX','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf','Db','Sg','Bh','Hs','Mt','Un']
ha2cmm1 = 219474.63

def get_sym_elem(chardom):
        """Retrieve symmetry elements from XML file"""
        isym=0
        sym_elem = []
        sym_tnons = []
        sym_elem_nodes = chardom.getElementsByTagName('genpos')
#       print " number of sym elements ", len(sym_elem_nodes)
        for sym_elem_node in sym_elem_nodes:
                #print sym_elem_node
                sym_elem.append([])
                sym_tnons.append([])
                rows = sym_elem_node.getElementsByTagName('row')
                for row in rows:
                        try:
                                data = row.childNodes[0].data
                                #print data
                                tmp = re.split ('[ \n\t\[\]]*',data)
                                tmp.pop()
                                sym_tnons[isym].append(tmp.pop())
                                tmp.pop(0)
                                #print tmp
                                sym_elem[isym].extend(map(string.atoi,tmp))
                        except AttributeError:
                                print "value error in get_sym_elem" 
                                raise
                                
                isym = isym+1
        
#       print sym_elem
        return sym_elem,sym_tnons

def get_irreps_chars(chardom,nsym):
        """Retrieve characters for sym elements from XML file"""
        irreps_node = chardom.getElementsByTagName('irreps')
        if (len(irreps_node) > 1):
                print "Error : more than one irreps field in XML"
                raise ValueError
        try :
                nirreps = irreps_node[0].getAttribute('num')
        except ValueError:
                print "no nirreps field in irreps tag"
                raise
#       print " number of irreps ", nirreps
        repr_nodes = irreps_node[0].getElementsByTagName('repr')
        irrep_character = []
        irrep_name = []
        iirrep=0
        for repr_node in repr_nodes:
                irrep_character.append([])
                #print repr_node
                irrep_name.append(repr_node.getAttribute('name'))
                dim = string.atoi(repr_node.getAttribute('dim'))
                for mat in repr_node.getElementsByTagName('mat'):
                        idim=0
                        character = 0.0
                        for row in mat.getElementsByTagName('row'):
                                try:
                                        data = row.childNodes[0].data
                                        #print data
                                        tmp = re.split ('[ 
\n\t\[\]\(\)\,]*',data)
                                        tmp.pop()
                                        tmp.pop(0)
                                        #print tmp
                                        char = 
int(float(tmp[2*idim])*math.cos(math.pi*float(tmp[2*idim+1])/180.0))
                                        #print char
                                        character = character + char
                                except AttributeError:
                                        print "value error in 
get_irreps_chars" 
                                        raise
                                
                                idim=idim+1
                        irrep_character[iirrep].append(character)
                iirrep=iirrep+1
        
#       print irrep_name,irrep_character
        return irrep_name,irrep_character

def get_abinit_syms(tokens):
        """Retrieve symmetry elements from abinit file tokenized previously"""
        
        ind_nsym = tokens.index('nsym')
        nsym = int (tokens[ind_nsym+2])
#       print "nsym = ", nsym
        ind_symrel = tokens.index('symrel')
        symrel = []
        ii = 0
        for isym in range(nsym):
                
symrel.append(map(int,tokens[ind_symrel+ii+1:ind_symrel+ii+10]))
                ii = ii+9
        ind_tnons = tokens.index('tnons')
        tnons = []
        ii = 0
        for itnons in range(nsym):
                tnons.append(map(float,tokens[ind_tnons+ii+1:ind_tnons+ii+3]))
                ii = ii+3
#       print "symrel = ", symrel
        return symrel

def get_abinit_group(tokens):
        """Retrieve symmetry group from abinit output file"""
        groupnum = 0
        try:
                ind_spgroup = tokens.index('spgroup')
                groupnum = tokens[ind_spgroup+1]
        except ValueError:
                ind_spgroup = tokens.index('space')
                for ii in range(10):
                        if tokens[ind_spgroup+1] == 'group':
                                while not 
re.compile('\(\#[0-9]*\)').match(tokens[ind_spgroup]):
                                        ind_spgroup = ind_spgroup+1
                                groupnum = 
int(re.sub('[\(\)\#;]*','',tokens[ind_spgroup]))
                                break
                        else:
                                tokens.pop(ind_spgroup)
#       print groupnum
        return groupnum
#               if not re.compile('degenerate').match(tokens[ind1]):

                        
def get_abinit_chars(tokens):
        """Retrieve characters and frequencies for modes from anaddb output 
file"""

        for ii in range(tokens.count('Phonon')):
                ind1 = tokens.index('Phonon')
                if tokens[ind1+1] != 'wavevector' or tokens[ind1+2] != 
'(reduced' :
                        tokens.pop(ind1)
                else:
                        if (float(tokens[ind1+5]) == 0.0 and \
                                float(tokens[ind1+6]) == 0.0 and \
                                float(tokens[ind1+7]) == 0.0 ) :
                                tokens.pop(ind1)
                        break
        abi_freq = []
        ind1=ind1+12
#       print tokens[ind1],tokens[ind1+1]
        while 
re.compile('[-]?[0-9].[0-9]*[-]?[Ee]?[0-9]*').match(tokens[ind1]):
                abi_freq.append(float(tokens[ind1]))
#               print abi_freq
                ind1=ind1+1

        for ii in range(tokens.count('Analysis')):
                ind1 = tokens.index('Analysis')
                if tokens[ind1+2] != 'degeneracies' or tokens[ind1+4] != 
'characters':
                        tokens.pop(ind1)
                else:
                        break
        while not re.compile('Symmetry').match(tokens[ind1]):
                ind1=ind1+1
        imode=0
        abi_characters = []
        while not re.compile('=[=]*').match(tokens[ind1]):
                abi_characters.append([])
                while not re.compile('[0-9]').match(tokens[ind1]):
                        ind1=ind1+1
                ind1=ind1+1
                #print "tokens now ", tokens[ind1:ind1+20]
                #print "abi_characters ", abi_characters
                degen = 0
                if not re.compile('degenerate').match(tokens[ind1]):
                        while 
re.compile('[-]?[0-9]+\.[0-9]+').match(tokens[ind1]):
                                charlist = int(float(tokens[ind1]))
                                abi_characters[imode].append(charlist)
                                ind1 = ind1+1
                else:
                        degen = 1
                        while not 
re.compile('[-]?[0-9]+\.[0-9]+').match(tokens[ind1]):
                                ind1=ind1+1
                        while 
re.compile('[-]?[0-9]+\.[0-9]+').match(tokens[ind1]):
                                charlist = int(float(tokens[ind1]))
                                abi_characters[imode].append(charlist)
                                ind1 = ind1+1
                        
                        
                if degen == 0:
                        imode=imode+1
                elif degen == 1:
                        abi_characters.append(abi_characters[imode])
                        imode=imode+2

#       print "abi_characters = ", abi_characters,abi_freq
        return abi_characters,abi_freq

        
def get_sym_corresp(abinit_sym_elem,bilbao_sym_elem):
        """determine correspondence between Bilbao and abinit symmetry 
lists"""
        abi_to_bilbao = []
        nsym = len(abinit_sym_elem)
        if nsym != len(bilbao_sym_elem):
                print "Error: different number of syms in abinit/bilbao"
                raise ValueError
        for iabi in range(nsym):
                for ibil in range(nsym):
                        found = 1
                        for ielem in range(9):
                                if abinit_sym_elem[iabi] != 
bilbao_sym_elem[ibil]:
                                        found=0
                                        break
                        if found == 1:
                                abi_to_bilbao.append(ibil)
                                break
        if len(abi_to_bilbao) != nsym:
                print "Error : not all correspondences found"
                raise ValueError
#       print abi_to_bilbao
        return abi_to_bilbao

def get_bilbao_symbols(bilbao_sym_elem,bilbao_irreps_chars):
        """Deteremine irreps names from symmetry considerations"""
        
        iident = find_ident(bilbao_sym_elem)
        iinv = -1
        try:
                iinv = find_inv(bilbao_sym_elem)
                has_inv = 1
        except ValueError:
                print "Inv not found"
                has_inv = 0
        irot,rot_order = find_princ_rotation(bilbao_sym_elem,has_inv,iinv)
#       print " irot,rot_order = ",irot,rot_order
        if has_inv==1:
                print "ident,inv,irot,rot_order = ", 
iident,iinv,irot,rot_order
        elif has_inv==0:
                print "ident,irot,rot_order = ", iident,irot,rot_order

        ia = 1
        ib = 1
        ie = 1
        it = 1

        bilbao_symbols=[]
        for iirrep in range(len(bilbao_irreps_chars)):
                if bilbao_irreps_chars[iirrep][iident] == 1:
                        if bilbao_irreps_chars[iirrep][irot] == 1:
                                irrepletter = 'A'
                                irrepnum = ia
                                ia = ia+1
                        elif bilbao_irreps_chars[iirrep][irot] == -1:
                                irrepletter = 'B'
                                irrepnum = ib
                                ib = ib+1
                        else:
                                print "couldn't classify the irrep into A or 
B"
                                raise ValueError
                elif bilbao_irreps_chars[iirrep][iident] == 2:
                        irrepletter = 'E'
                        ie = ie+1
                        irrepnum = ie
                elif bilbao_irreps_chars[iirrep][iident] == 3:
                        irrepletter = 'T'
                        it = it+1
                        irrepnum = ie
                else:
                        print "couldn't classify the irrep into A,B,E, or T"
                        raise ValueError
                
                geradeletter = ' '
                if has_inv==1:
                        if bilbao_irreps_chars[iirrep][iinv] > 0:
                                geradeletter = 'g'
                        elif bilbao_irreps_chars[iirrep][iinv] < 0:
                                geradeletter = 'u'
                        irrepnum = int((irrepnum+1.0)/2.0)
                        
#               print "%1c%1d%1c" % (irrepletter,irrepnum,geradeletter)
                bilbao_symbols.append("%1c%1d%1c" % 
(irrepletter,irrepnum,geradeletter))

        return bilbao_symbols

def find_ident(sym_elem):
        """find the index of the identity operation in the symmetries"""
        for isym in range(len(sym_elem)):
                if sym_elem[isym] == [1,0,0,0,1,0,0,0,1]:
                        return isym
        print "Error : ident not found"
        raise ValueError

def find_inv(sym_elem):
        """find the index of the inverse operation in the symmetries, if it 
exists"""
        for isym in range(len(sym_elem)):
                if sym_elem[isym] == [-1,0,0,0,-1,0,0,0,-1]:
                        return isym
#       print "Error : inv not found"
        raise ValueError

def find_princ_rotation(sym_elem,has_inv,iinv):
        """Find the principal rotation of the group, as the highest order 
member which doesn't give -I"""
        morder = 0
        msym = 0
        for isym in range(len(sym_elem)):
                is_Sop = 0
                #if has_inv==1:
                #       if 
                # suppose that order is at most 6 for elements of sym
                tmpsym = [1,0,0,0,1,0,0,0,1]
                for iorder in range(1,7):
#                       print "iorder= ", iorder
                        tmpsym = sym_elem_prod(tmpsym,sym_elem[isym])
                        if tmpsym == [-1,0,0,0,-1,0,0,0,-1]:
                                is_Sop = 1
                                break
                        if tmpsym == [1,0,0,0,1,0,0,0,1]:
                                break
                # if we found an S operation (axis+Inverse) dont count it
                #   should have a corresponding C operation
                #   the danger is that order(S) = 2*order(C) and we want the C
                #   to deterine the principal rotation
                if is_Sop == 1:
                        break
                if iorder > morder:
                        msym = isym
                        morder = iorder
        return msym,morder

def sym_elem_prod (s1,s2):
        """ calculate (non commuting) product of s2*s1"""
        sym = []
        for ii in range(3):
                for jj in range(3):
                        ss = 0
                        for kk in range(3):
                                ss = ss + s2[3*ii+kk]*s1[3*kk+jj]
                        sym.append(ss)
        return sym

def print_abi_irreps 
(abi_to_bilbao,abinit_chars,abinit_freq,irreps_names,irreps_chars):
        """Print out irreps for each abinit mode"""
        
        
        order = len(irreps_chars[0])
#       print "abinit_chars = ", abinit_chars
        for ichar in range(len(abinit_chars)):
                freqcmm1 = abinit_freq[ichar]*ha2cmm1
                print "Abinit mode # ", "%2d = " % ichar, " %15.2f cm-1 " % 
freqcmm1,
                for iirrep in range(len(irreps_chars)):
                        sum = 0
                        for isym in range(len(irreps_chars[iirrep])):
                                sum = sum + 
irreps_chars[iirrep][abi_to_bilbao[isym]]*abinit_chars[ichar][isym]
                        if sum != 0:
                                print int(sum/order), "x", 
irreps_names[iirrep],
                print
        

if __name__=='__main__':
#       for ielem in range(len(elemabbrev)):
#               print ielem, elemabbrev[ielem]

#
        if len(sys.argv) < 2:
                print "Error : need two filename arguments : abinit output 
file and anaddb output file"
        else:
#               charfilename = sys.argv[1]
#               groupnum = int(sys.argv[1])
                abifilename = sys.argv[1]
                anaddbfilename = sys.argv[2]

#
#  get group and symmetry operations from abinit output file
#
        abitokens=tokenize_file(abifilename)
        #print abitokens
        
        groupnum = get_abinit_group(abitokens)
        
        abinit_sym_elem = get_abinit_syms(abitokens)
        
#
#  get characters of different modes from anaddb output file. Get the 
information
#    only for the Gamma point.
#
        anaddbtokens=tokenize_file_keep_comments(anaddbfilename)
        #print anaddbtokens

        abinit_chars,abinit_freq = get_abinit_chars(anaddbtokens)
        
#
#   open url and request xml input for characters and symops of 
#     point group "groupnum"
#
        print "http://www.cryst.ehu.es"+\
                "/cgi-bin/cryst/xml/nph-get_doc?p=irreps&g=%d&k=0,0,0,GM,p\n" 
% groupnum
        try:
                fhan = urllib.urlopen("http://www.cryst.ehu.es"+\
                        
"/cgi-bin/cryst/xml/nph-get_doc?p=irreps&g=%d&k=0,0,0,GM,p\n" % groupnum)
        except IOError:
                print "Error opening http://www.cryst.ehu.es";
                return -1
        
        bilbao_xml_lines = fhan.readlines()
        XML = ''
        for line in bilbao_xml_lines:
                XML = XML + string.strip(line)
        

#
#   parse the XML input obtained from www.cryst.ehu.es Bilbao server
#    use the minimalistic minidom. No dtd support in minidom, although 
#    one exists on the server.
#
        chardom = xml.dom.minidom.parseString(XML)
        #print chardom.documentElement.tagName
        if (chardom.documentElement.hasAttributes()):
                print chardom.documentElement.attributes
#       else:
#               print "No attributes here!"


#
#  extract symmetry elements from Bilbao xml,
#
        bilbao_sym_elem,bilbao_sym_tnons = get_sym_elem(chardom)

#
#  extract irreps from Bilbao xml,
#
        
bilbao_irreps_names,bilbao_irreps_chars=get_irreps_chars(chardom,len(bilbao_sym_elem))

#
#  get correspondance between abinit and Bilbao symmetry elements
#
        abi_to_bilbao = get_sym_corresp(abinit_sym_elem,bilbao_sym_elem)

#
#  generate names of irreps A1u ... from the characters of each irrep for -I, 
main rotation axis...
#
        bilbao_irreps_names = 
get_bilbao_symbols(bilbao_sym_elem,bilbao_irreps_chars)

#
#  print out abinit frequencies and irreps for each phonon mode.
#
        print_abi_irreps 
(abi_to_bilbao,abinit_chars,abinit_freq,bilbao_irreps_names,bilbao_irreps_chars)
        

#! /usr/bin/python

#
#    Copyright (C) 2003 ABINIT group
#
#    Written by Matthieu Verstraete in python (compatible v1.9).
#    This is free software, and you are welcome to redistribute it
#    under certain conditions (GNU General Public License,
#    see ~ABINIT/Infos/copyright or http://www.gnu.org/copyleft/gpl.txt).
#
#    ABINIT is a project of the Universite Catholique de Louvain,
#    Corning Inc. and other collaborators, see ~ABINIT/Infos/contributors.
#    Please read ~ABINIT/Infos/acknowledgments.htm for suggested
#    acknowledgments of the ABINIT effort.
#
#    For more information, see http://www.abinit.org .
#
#  This module makes a file into a list of tokens, conserving
#   or eliminating comments of the type "#  bla bla   \n"
#  Probably redundant with other, more powerful python modules, but hey.
#


def tokenize_file(infilename):
        """Take input file name and return list of tokens, eliminating 
comments"""
#  reg expss
        import re
        import sys
        import fileinput
        import math

#
#  open and parse input 
#

        try:
                fp = open (infilename, 'r')
        except IOError:
                print "Error opening file"
                raise

        lines = fp.readlines ()

#
#  put all tokens into tokens and remove comments
#
        tokens = []
        for line in lines:
                tmp = re.split ('[ \t\n]*',line)
#               print "tmp = ", tmp
                for tok in tmp:
                        if (tok != ''):
                                if (re.compile('[#!][.]*').match(tok)):
                                        break
                                tokens.append(tok)
#       print "tokens = ", tokens

        fp.close()

        return tokens

def tokenize_file_keep_comments(infilename):
        """Take input file name and return list of tokens, keep comments"""
#  reg expss
        import re
        import sys
        import fileinput
        import math

#
#  open and parse input 
#

        try:
                fp = open (infilename, 'r')
        except IOError:
                print "Error opening file"
                raise

        lines = fp.readlines ()

#
#  put all tokens into tokens and remove comments
#
        tokens = []
        for line in lines:
                tmp = re.split ('[ \t\n]*',line)
#               print "tmp = ", tmp
                for tok in tmp:
                        if (tok != ''):
                                tokens.append(tok)
#       print "tokens = ", tokens

        fp.close()

        return tokens