The ABINIT Users Mailing List ( CLOSED )

[Matthieu Verstraete <mjv500@york.ac.uk>]

Further update:

The Bilbao web site gives conventional basis symops, so I added the double 
check to get_irreps: if the reduced symops don't work, tries the cartesian 
ones (attached).

If you publish anything using this script or the Bilbao web site, please 
cite:

M. I. Aroyo, A. Kirov, C. Capillas, J. M. Perez-Mato & H. Wondratschek.
"Bilbao Crystallographic Server II: Representations of crystallographic
point groups and space groups". Acta Cryst. (2006), A62, 115-128

If you want Raman and IR active modes, they have an online program for 
this too: http://www.cryst.ehu.es/rep/sam.html

It is forbidden to find bugs in my code for the next month!

Matthieu

--
================================================================
Dr. Matthieu Verstraete                 mailto:mjv500@york.ac.uk
Dept. of Physics, University of York,     tel: +44 1904 43 22 08
Heslington, YO10 5DD York, United Kingdom fax: +44 1904 43 22 14#! /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
        if (len(irreps_node) < 1):
                print "Error : no 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
                                        #print "tmp[2*idim],tmp[2*idim+1] ", 
tmp[2*idim],tmp[2*idim+1]
                                        char = 
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_names and irrep_character"
#       for irrepnum in range(len(irrep_name)):
#               print irrep_name[irrepnum],irrep_character[irrepnum]

        return irrep_name,irrep_character

def get_abinit_basis_vectors(tokens):
        """Retrieve real and reciprocal space lattice vectors from abinit 
file tokenized previously"""

        ind_rgprim = tokens.index('Real(R)+Recip(G)')
        rprim = []
        rprim.append(map(float,tokens[ind_rgprim+8:ind_rgprim+11]))
        rprim.append(map(float,tokens[ind_rgprim+16:ind_rgprim+19]))
        rprim.append(map(float,tokens[ind_rgprim+24:ind_rgprim+27]))

        gprim = []
        gprim.append(map(float,tokens[ind_rgprim+12:ind_rgprim+15]))
        gprim.append(map(float,tokens[ind_rgprim+20:ind_rgprim+23]))
        gprim.append(map(float,tokens[ind_rgprim+28:ind_rgprim+31]))
        #print "rprim, gprim = ", rprim, gprim
        return rprim,gprim

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
        tnons = []
        try:
                ind_tnons = tokens.index('tnons')
                ii = 0
                for itnons in range(nsym):
                        
tnons.append(map(float,tokens[ind_tnons+ii+1:ind_tnons+ii+3]))
                        ii = ii+3
        except ValueError:
                for itnons in range(nsym):
                        tnons.append((0.0,0.0,0.0))
                
#       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 = int(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"""

#  find the correct occurence of Phonon which is followed by the gamma point 
freqs
        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 gamma point
                        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]):
        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
                imode_file = int(tokens[ind1])-1
                if (imode_file != imode):
                        print "Probably error : modes don't match ", imode, 
imode_file
                ind1=ind1+1
                #print "tokens now ", tokens[ind1:ind1+20]
                #print "abi_characters ", abi_characters
                degen = 1
                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:
#  Find degeneracy of present mode
                        while not  re.compile('[0-9]+').match(tokens[ind1]):
                                ind1=ind1+1
                        if re.compile('to').match(tokens[ind1+1]):
                                ind1=ind1+2
                        jmode = int(tokens[ind1])
                        degen = jmode-imode
                        
                        ind1=ind1+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
                        
                        
# Duplicate characters for degenerate modes
                for ii in range(degen-1):
                        abi_characters.append(abi_characters[imode])
                imode=imode+degen

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

def abinit_syms_to_cart(abinit_sym_elem,rprim,gprim):
        """Change abinit symops to cartesian coordinates"""

        abinit_sym_elem_cart = []
        for isym in range(len(abinit_sym_elem)):
                abinit_sym_elem_cart.append([0, 0, 0, 0, 0, 0, 0, 0, 0])
                tmpmat  = [[0., 0., 0.],[0., 0., 0.],[0., 0., 0.]]
                tmpmat2 = [[0., 0., 0.],[0., 0., 0.],[0., 0., 0.]]
                for irow in range(3):
                        for icol in range(3):
                                for iin in range(3):
                                        tmpmat[irow][icol] += 
abinit_sym_elem[isym][irow*3+iin]*rprim[iin][icol]
#                                       tmpmat[irow][icol] += 
abinit_sym_elem[isym][irow+iin*3]*gprim[icol][iin]

                for irow in range(3):
                        for icol in range(3):
                                for iin in range(3):
                                        tmpmat2[irow][icol] += 
gprim[irow][iin]*tmpmat[iin][icol]
#                                       tmpmat2[irow][icol] += 
rprim[iin][irow]*tmpmat[iin][icol]

#
#  round off: should be all 0 1 or -1
#
                for irow in range(3):
                        for icol in range(3):
                                abinit_sym_elem_cart[isym][irow*3+icol] = 
round(tmpmat2[irow][icol])
        return abinit_sym_elem_cart

        
def get_sym_corresp(abinit_sym_elem,abinit_sym_elem_cart,bilbao_sym_elem):
        """determine correspondence between Bilbao and abinit symmetry 
lists"""
        abi_to_bilbao = []
        nsym = len(abinit_sym_elem_cart)
#       print " nsyms = ", nsym, len(bilbao_sym_elem)

#       print "abinit_sym_elem_cart = "
#       for symelem in abinit_sym_elem_cart:
#               print 
(symelem[0]+1)*100000000+(symelem[1]+1)*10000000+(symelem[2]+1)*1000000+(symelem[3]+1)*100000+(symelem[4]+1)*10000+(symelem[5]+1)*1000+(symelem[6]+1)*100+(symelem[7]+1)*10+(symelem[8]+1),
 symelem
#       print "bilbao_sym_elem = "
#       for symelem in bilbao_sym_elem:
#               print 
(symelem[0]+1)*100000000+(symelem[1]+1)*10000000+(symelem[2]+1)*1000000+(symelem[3]+1)*100000+(symelem[4]+1)*10000+(symelem[5]+1)*1000+(symelem[6]+1)*100+(symelem[7]+1)*10+(symelem[8]+1),
 symelem
        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 found != 1:
                        print "Error: didnt find symop ", iabi, " = ", 
abinit_sym_elem[iabi], " in reduced sym list"

#  if not all sym correspondences were found, we may be in a 
#    conventional basis instead of a reduced primitive one:
        if len(abi_to_bilbao) != nsym:
                print "Not all correspondences found for reduced syms. Try 
cartesian"
                print "If your system is not cubic or hexagonal, probably an 
error."
                abi_to_bilbao = []
                for iabi in range(nsym):
                        for ibil in range(nsym):
                                found = 1
                                for ielem in range(9):
                                        if abinit_sym_elem_cart[iabi] != 
bilbao_sym_elem[ibil]:
                                                found=0
                                                break
                                if found == 1:
                                        abi_to_bilbao.append(ibil)
                                        break
                        if found != 1:
                                print "Error: didnt find symop ", iabi, " = 
", abinit_sym_elem_cart[iabi], " in cartesian sym list"
                
                
        if len(abi_to_bilbao) != nsym:
                print "Error : not all correspondences found"
                print "abi_to_bilbao = "
                print abi_to_bilbao
                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

        print "List of symbols for Bilbao irreps:"
        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'
                        irrepnum = ie
                        ie = ie+1
                elif bilbao_irreps_chars[iirrep][iident] == 3:
                        irrepletter = 'T'
                        irrepnum = it
                        it = it+1
                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
                ichar_eff=ichar+1
                print "Abinit mode # ", "%2d = " % ichar_eff, " %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 math.fabs(sum) > 0.001:
                                #print sum
                                print int(sum/order), "x", 
irreps_names[iirrep], '(Irrep num. ', iirrep+1, ')',
                print
        

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

#
#       print "sys.argv = ", sys.argv
        if len(sys.argv) < 2:
                print "Error : need two filename arguments : abinit output 
file and anaddb output file"
                raise ValueError
        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)
        print "Space group", groupnum
        
        abinit_sym_elem = get_abinit_syms(abitokens)
        
        
        rprim,gprim = get_abinit_basis_vectors(abitokens)
        
        abinit_sym_elem_cart = 
abinit_syms_to_cart(abinit_sym_elem,rprim,gprim)

#
#  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&gor=2\n" % 
groupnum)
        except IOError:
                print "Error opening http://www.cryst.ehu.es";
                raise
        
        bilbao_xml_lines = fhan.readlines()
        XML = ''
        for line in bilbao_xml_lines:
                XML = XML + string.strip(line)
                #print bilbao_xml_lines
        

#
#   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,abinit_sym_elem_cart,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)