#/*##########################################################################
# Copyright (C) 2011-2014 European Synchrotron Radiation Facility
#
# Ab2tds
# European Synchrotron Radiation Facility, Grenoble,France
#
# Ab2tds is developed at
# the ESRF by the Scientific Software staff.
# Principal author for Ab2tds: Alessandro Mirone, Bjorn Wehinger
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 2 of the License, or (at your option)
# any later version.
#
# Ab2tds is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# Ab2tds; if not, write to the Free Software Foundation, Inc., 59 Temple Place,
# Suite 330, Boston, MA 02111-1307, USA.
#
# Ab2tds follows the dual licensing model of Trolltech's Qt and Riverbank's PyQt
# and cannot be used as a free plugin for a non-free program.
#
# Please contact the ESRF industrial unit (industry@esrf.fr) if this license
# is a problem for you.
#############################################################################*/
import math
import string
import sys
import time
import numpy
import hashlib
import os
from yaml import load, dump
try:
from yaml import CLoader as Loader, CDumper as Dumper
except ImportError:
from yaml import Loader, Dumper
USE_HDF5=1
SMALLQ=1.0e-4
if USE_HDF5:
import h5py
else:
import cPickle
def ciclicity(a,b,c):
if( a==b or a==c or c==b):
return 0
elif( (b-a)*(c-b)*(a-c) <0 ):
return 1
else:
return -1
def CellVolume( cellVects ):
AntiSymm= numpy.array([ [ [ ciclicity(i,j,k) for k in range(0,3) ] for j in range (0,3) ] for i in range(0,3) ])
cellVolume=numpy.dot( numpy.dot(AntiSymm, cellVects[2]),cellVects[1] )
cellVolume=numpy.dot(cellVects[0],cellVolume)
return cellVolume
def BrillVects(cellVects):
AntiSymm= numpy.array([ [ [ ciclicity(i,j,k) for k in range(0,3) ] for j in range (0,3) ] for i in range(0,3) ])
cellVolume= CellVolume(cellVects)
Brillvectors=numpy.zeros([3,3], numpy.float32)
Brillvectors[0]=2 * numpy.pi * numpy.dot( numpy.dot(AntiSymm, cellVects[2]),cellVects[1] )/ cellVolume
Brillvectors[1]=2 * numpy.pi * numpy.dot( numpy.dot(AntiSymm, cellVects[0]),cellVects[2] )/ cellVolume
Brillvectors[2]=2 * numpy.pi * numpy.dot( numpy.dot(AntiSymm, cellVects[1]),cellVects[0] )/ cellVolume
return Brillvectors
class CalcDatas:
@staticmethod
def get_CalcDatas_Object_From_File(filename, file_extra=None,
energy_scaling=1):
tipo = "Castep"
try:
class params:
s=open(sys.argv[2], "r").read()
exec(s)
tipo = params.CODE
except:
pass
returndata= ReadCastep(filename, file_extra,
energy_scaling=energy_scaling, tipo=tipo)
return returndata
def Reduce(self, Qs):
# res=numpy.dot(Qs,self.cellVects )/(2 * numpy.pi)
bv = self.GetBrillVects()
res = numpy.dot( Qs, numpy.linalg.inv (bv ) )
return res
def GetCellVolume(self):
cellVolume = CellVolume( self.cellVects )
return cellVolume
def GetBrillVects(self, CTransf=None):
if CTransf is None:
Brillvectors= BrillVects(self.cellVects)
else:
cellVects = numpy.dot(CTransf, self.cellVects)
Brillvectors= BrillVects(cellVects)
return Brillvectors
def GetClippedQs(self, QsReduced=None):
Brillvectors = self.GetBrillVects()
if QsReduced is None:
QsReduced=self.QsReduced
else:
QsReduced =numpy.array(QsReduced, copy=True )
pass
QsReduced[:] = ( (QsReduced +0.5 ) % 1) -0.5 # 0.9 1.4 0.4 -0.1
## positivita = numpy.less( QsReduced , 0.0)
## if numpy.sum(numpy.sum(positivita) ):
## raise Exception, " procedura of reduction of Qs inside the Brillouin zone by now works only with positive coefficients "
Qs_Brill = (numpy.dot( QsReduced , Brillvectors ) )
mods_ref= numpy.sum( Qs_Brill* Qs_Brill, axis=-1 )
factors_list= [ [i,j,k] for i in range(-1,2) for j in range(-1,2) for k in range(-1,2)]
a = numpy.zeros( Qs_Brill.shape, Qs_Brill.dtype)
Qs_new= numpy.zeros( Qs_Brill.shape, Qs_Brill.dtype)
a[:]=Qs_Brill
for (i,j,k) in factors_list:
Qs_new[:] = Qs_Brill - (i* Brillvectors[0] + j*Brillvectors[1] + k*Brillvectors[2] )
mods = numpy.sum( Qs_new * Qs_new, axis=-1 )
mask = numpy.less(mods, mods_ref)
notmask= numpy.less(mask , 1)
a[mask] = Qs_new[mask]
a[notmask ] = a[ notmask]
mods_ref[mask] = mods[mask]
return a
def GetNonClippedQs(self, QsReduced=None):
Brillvectors = self.GetBrillVects()
if QsReduced is None:
QsReduced=self.QsReduced
else:
pass
Qs_Brill = (numpy.dot( QsReduced , Brillvectors ) )
return Qs_Brill
def __init__(self, NofIons=None, NofBranches=None, NofWaves=None,
cellVects=None, atomReducedPositions=None, atomNames=None, atomMasses=None,
frequencies=None, eigenvectors=None, weights=None, atomAbsolutePositions=None,
QsReduced=None, namemd5=None):
self.gvectsdump=None
if NofIons is None:
return
self.NofIons= NofIons
self.NofBranches= NofBranches
self.NofWaves= NofWaves
self.cellVects= cellVects
self.atomReducedPositions= atomReducedPositions
self.atomNames= atomNames
self.atomMasses= atomMasses
self.frequencies= frequencies
self.eigenvectors= eigenvectors
self.weights= weights
self.QsReduced=QsReduced
self.namemd5=namemd5
self.gvectsdump=None
if atomAbsolutePositions is None:
self.atomAbsolutePositions=[]
#############################################################"""
# The positions are given in the cell vectos basis.
# Now we transform the position in xyz absolute space.
for k in range(0, len(atomReducedPositions )):
self.atomAbsolutePositions.append(numpy.dot(atomReducedPositions[k], cellVects ))
self.atomAbsolutePositions =numpy.array(self.atomAbsolutePositions)
assert( NofIons == len( self.atomAbsolutePositions ) )
assert( NofBranches == eigenvectors.shape[1] )
assert( NofBranches == eigenvectors.shape[2] )
assert( NofWaves == eigenvectors.shape[0] )
assert( NofWaves == frequencies.shape[0] )
assert( NofBranches == frequencies.shape[1] )
assert( NofIons == len( atomMasses ) )
def Get_md5postfix(self):
lfilename=string.split(self.namemd5,".")
postpend=lfilename[-1]
if "md5_code=" in postpend and len(postpend[9:])==8*4:
return postpend
else:
raise Exception, " not able to retrieve the md5 postfix"
def Get_filename(self):
lfilename=string.split(self.namemd5,".")
postpend=lfilename[-1]
if "md5_code=" in postpend and len(postpend[9:])==8*4:
return self.namemd5[:-(8*4+1+9)]
else:
return self.namemd5
def dump_on_file(self, nomefile):
h=h5py.File(nomefile,'w')
h['/NofIons'] = self.NofIons
h['/NofBranches'] = self.NofBranches
h['/NofWaves'] = self.NofWaves
h['/cellVects'] = self.cellVects
h['/atomReducedPositions'] = self.atomReducedPositions
h['/atomAbsolutePositions'] = self.atomAbsolutePositions
h['/atomNames'] = self.atomNames
h['/atomMasses'] = self.atomMasses
h['/frequencies'] = self.frequencies
h['/eigenvectors'] = self.eigenvectors
h['/weights'] = self.weights
h['/Qs'] = self.QsReduced
h.flush()
h.close()
def read_from_file(self, nomefile):
self.namemd5=nomefile
h=h5py.File(nomefile,'r+')
self. NofIons= ( h['/NofIons']).value
self.NofBranches=(h['/NofBranches'] ).value
self.NofWaves=(h['/NofWaves'] ).value
self.cellVects=h['/cellVects'] [:]
self.atomReducedPositions=h['/atomReducedPositions'] [:]
self.atomAbsolutePositions=h['/atomAbsolutePositions'] [:]
self.atomNames=h['/atomNames'] [:]
self.atomMasses=h['/atomMasses'] [:]
self.frequencies=h['/frequencies'] [:]
self.eigenvectors=h['/eigenvectors'] [:]
self.weights=h['/weights'] [:]
self.QsReduced = h['/Qs'][:]
h.flush()
h.close()
def dump_extra_on_file(self, nomefile):
h5=h5py.File(nomefile,'r+')
res={"BornCharges": self.BornCharges, "Polariz":self.Polariz }
for keyd in res.keys():
nuova_key= keyd
if(nuova_key in h5):
del h5[nuova_key]
h5[nuova_key]= res[keyd]
h5.flush()
h5.close()
def read_extra_from_file(self, nomefile):
self.namemd5=nomefile
h=h5py.File(nomefile,'r+')
try:
self.BornCharges = ( h['/BornCharges']).value
self.Polariz = ( h['/Polariz'] ).value
h.flush()
h.close()
return 1
except:
print " DID NOT FIND BornCharges Polariz in pre-existing HDF5 file "
h.flush()
h.close()
return 0
def GetGvectsDump(self) :
if self.gvectsdump is not None:
return self.gvectsdump
Brillvectors = self.GetBrillVects()
gvects=[]
for i in range(-2,3):
for j in range(-2,3):
# i=j=0
for k in range(-2,3):
dum = i*Brillvectors[0]+j*Brillvectors[1]+k*Brillvectors[2]
gvects.append(dum)
gs = [ numpy.sum( b*b) for b in Brillvectors ]
G = min(gs)
dum= 30.0/G
# self.gvectsdump = numpy.array(gvects), dum
self.gvectsdump = numpy.array([[0,0,0.0]]), dum
return self.gvectsdump
def GetNonAnalytic( self, q ):
if( not hasattr(self,"BornCharges") ):
if len(q.shape)==1:
return numpy.zeros([ 3*len(self.atomNames), 3*len(self.atomNames)],"f")
else:
nqs=len(q)
return numpy.zeros([nqs, 3*len(self.atomNames), 3*len(self.atomNames)],"f")
Gvects, dumpfact = self.GetGvectsDump() ##
# pfs = numpy.tensordot( Gvects ,self.atomAbsolutePositions,[[-1],[-1]])
# pfs = numpy.exp( -1.0j*pfs)
if(len(q.shape)==1) :
q= Gvects + q
elif(len(q.shape)==2):
q= Gvects[:,None,:] + q
elif(len(q.shape)==3):
q= Gvects[:,None,None,:] + q
elif(len(q.shape)==4):
q= Gvects[:,None,None,None,:] + q
else:
raise Exception, "problem : shape unknown "
facts = numpy.sum(q*q,axis=-1)*dumpfact ##
facts=numpy.exp(-facts) ##
pfs = numpy.tensordot( q ,self.atomAbsolutePositions,[[-1],[-1]])
pfs = numpy.exp( -1.0j*pfs)
P = self.Polariz / abs(self.GetCellVolume()) # adimensional
dielectric = 4*math.pi * P
dielectric[0,0]+=1.0
dielectric[1,1]+=1.0
dielectric[2,2]+=1.0
lenshape=len(q.shape)
if lenshape==1:
qq = q[None,:]*q[:,None]
elif lenshape==2:
qq = q[:, None,:]*q[:, :,None]
elif lenshape==3:
qq = q[:,:,None,:]*q[:,:, :,None]
elif lenshape==4:
qq = q[:,:,:,None,:]*q[:,:,:, :,None]
elif lenshape==5:
qq = q[:,:,:,:,None,:]*q[:,:,:,:, :,None]
else:
raise Exception, "problem : shape unknown "
## tmp = numpy.tensordot( q , dielectric , axes=[[-1],[0]] )
deno = numpy.tensordot( qq , dielectric , axes=[[-2,-1],[0,1 ]] )
Z = numpy.reshape(self.BornCharges, [-1,3,3])
# print Z
Z=numpy.transpose( numpy.transpose( Z)/ numpy.sqrt(self.atomMasses ) )
avdie = ( dielectric[0,0]+dielectric[1,1] +dielectric[2,2])/3.0
ZZself = Z[:,:,None,:]* Z[:,None,:,:]
bohr = 0.52917721092
# for i in range(len(ZZself)):
# print numpy.sum(ZZself[i], axis=-1)
# print "########### "
# raise Exception , " OK "
ZZself = numpy.sum( ZZself, axis=-1) *(4.0*math.pi/3.0*(1.0/2.0/math.pi/(dumpfact/2.0))**(3.0/2.0) /avdie )
# Z = numpy.reshape( Z , [-1,3])
Zq = numpy.tensordot( q , Z, axes=[[-1],[-1]] )
# if len(Zq.shape)>3:
# Zq= numpy.array(numpy.swapaxes( numpy.transpose( numpy.transpose(numpy.swapaxes ( Zq ,1,-2 ))* numpy.transpose(pfs) ) ,1,-2 ) )
# else:
# Zq=numpy.array(numpy.transpose( numpy.transpose(Zq)* numpy.transpose(pfs) ))
# numpy.transpose( numpy.transpose(Zq)* numpy.transpose(pfs) )
Zq=numpy.reshape( Zq , tuple( list(Zq.shape[:-2]) +[Zq.shape[-2]* Zq.shape[-1]]))
lenshape=len(Zq.shape)
if lenshape==1:
nume = Zq[None,:]*(Zq[:,None].conj())
elif lenshape==2:
nume = Zq[:,None,:]*(Zq[:,:,None].conj())
elif lenshape==3:
nume = Zq[:,:,None,:]*(Zq[:,:,:,None].conj())
elif lenshape==4:
nume = Zq[:,:,:,None,:]*(Zq[:,:,:,:,None].conj())
elif lenshape==5:
nume = Zq[:,:,:,:,None,:]*(Zq[:,:,:,:,:,None].conj())
else:
raise Exception, "problem : shape unknown "
res = 4*math.pi / abs(self.GetCellVolume()) * numpy.transpose (( numpy.transpose(facts)* numpy.transpose(nume))/numpy.transpose(deno)) ##
res=numpy.sum(res,axis=0) ##
# for i in range(len(self.atomNames)):
# if lenshape==2:
# res[i*3:(i+1)*3, i*3:(i+1)*3] = res[i*3:(i+1)*3, i*3:(i+1)*3] - ZZself[i]
# elif lenshape==3:
# res[:,i*3:(i+1)*3, i*3:(i+1)*3] = res[:,i*3:(i+1)*3, i*3:(i+1)*3] - ZZself[i]
# elif lenshape==4:
# res[:,:,i*3:(i+1)*3, i*3:(i+1)*3] = res[:,:,i*3:(i+1)*3, i*3:(i+1)*3] - ZZself[i]
# elif lenshape==5:
# res[:,:,:,i*3:(i+1)*3, i*3:(i+1)*3] = res[:,:,:,i*3:(i+1)*3, i*3:(i+1)*3] - ZZself[i]
# else:
# raise Exception, "problem : shape unknown "
res = res /4 / math.pi/ math.pi
bohr = 0.52917721092
res = res *( bohr ) # pour la distance qui s'en va avec e2
hartree = 4.3597482e-11
nuclearunitmass = 1.6605402e-24
c_light = 2.997925e+10
res=res*hartree/nuclearunitmass/c_light/c_light *1.0e16 # le dernier pour L2 au denominateur
return res
[docs]def ReadCastep(filename, filename_extra=None,
energy_scaling=1, tipo="Castep"):
# check if filename is already pospended by a md5code
lfilename=string.split(filename,".")
postpend=lfilename[-1]
res=None
resExtra=None
## try to read preparsed data for phonon from a hdf5 file (res)
if "md5_code=" in postpend and len(postpend[9:])==8*4:
namemd5=filename
else:
stime=time.time()
print "-----reading the file ", filename
s=open(filename,"r").read()
etime = time.time()
print " Reading Took ",etime - stime, "seconds"
m=hashlib.md5()
m.update(s)
postfix=m.hexdigest()
namemd5=filename+".md5_code="+postfix
if os.path.exists(namemd5):
print " ()()()() Pre-serialised object for this file exists already in file ", namemd5
s=None
print " ()() !! Now reloading the structure "
stime=time.time()
if(USE_HDF5):
print " ()() !! using HDF5"
res=CalcDatas()
res.read_from_file(namemd5)
else:
print " ()() !! using cPickle"
f=open(namemd5,"r")
res=cPickle.load(f)
etime = time.time()
print " Reloading Took ",etime - stime, "seconds"
# end of try from hdf5 for phonons
if res is None : # start parsing phonon file
if tipo=="Castep":
stime=time.time()
print "----- splitting the file into lines ", filename
sl=string.split(s,"\n")
s=None
etime = time.time()
print " splitting Took ",etime - stime, "seconds"
linecount=0
while( "Number of ions" not in sl[linecount]):
linecount+=1
NofIons= string.atoi(string.split(sl[linecount])[-1])
while( "Number of branches" not in sl[linecount]):
linecount+=1
NofBranches= string.atoi(string.split(sl[linecount])[-1])
while( "Number of wavevectors" not in sl[linecount]):
linecount+=1
NofWaves= string.atoi(string.split(sl[linecount])[-1])
print "\tNofIons\t\tNofBranches\tNofWaves"
print "\t", NofIons, "\t\t",NofBranches, "\t\t",NofWaves
while( "Unit cell vectors (A)" not in sl[linecount]):
linecount+=1
cellVects=[]
for icellvect in range(3):
linecount+=1
cellVects.append(map(string.atof, string.split(sl[linecount])))
cellVects=numpy.array(cellVects)
print "\tCell vectors "
for v in cellVects:
print "\t%e20.7 %e20.7 %e20.7"%tuple(v.tolist() )
while( "Fractional Co-ordinates" not in sl[linecount]):
linecount+=1
atomReducedPositions=[]
atomNames =[]
atomMasses =[]
for i in range(NofIons):
linecount+=1
items=string.split(sl[linecount])
atomReducedPositions.append(map(string.atof, items[1:4]))
atomNames.append(items[4])
atomMasses.append(string.atof(items[5]))
atomMasses=numpy.array(atomMasses)
qs= numpy.zeros([NofWaves,3],"f")
weights = numpy.zeros([NofWaves],"f")
frequencies = numpy.zeros([NofWaves,NofBranches],"f")
eigenvectors = numpy.zeros([NofWaves,NofBranches, NofBranches],"F")
linecount+=1
stime = time.time()
print "Read waves\t ", 0, "/", NofWaves,
oldQn=-1
while(oldQn<NofWaves-1): # CASTEP counting starts from 1
#for iwave in range(NofWaves):
iwave=oldQn
if (iwave and iwave%100 ==0):
etime = time.time()
print "\rRead waves\t ", iwave, "/", NofWaves, "in ", etime-stime, "seconds. Total is approx.",NofWaves*(etime-stime)/iwave,
sys.stdout.flush()
linecount+=1
items=string.split(sl[linecount], "=")
items=string.split(items[1])
newQn=string.atoi(items[0])-1 # CASTEP counting starts from 1
if(newQn==oldQn):
print " PRESENCE QPOINTS WITH MULTEPLICITY>1 . LINE = ", sl[linecount]
oldQn=newQn
qs[newQn]=map(string.atof, items[1:4])
if numpy.sum( qs[newQn]*qs[newQn] )==0.0 :
if(len(items)>8):
direction = numpy.array(map(string.atof, items[5:8]))
bvs = BrillVects(cellVects)
directionQ = numpy.dot(direction , bvs )
direction = SMALLQ * direction / numpy.sqrt( numpy.sum( directionQ*directionQ ) )
print " FOUND A Gamma POINT, replacing it with direction ", direction
qs[newQn]= direction/10
else:
qs[newQn]= [0.0,0.0,1.0e-6]
weights[newQn]=string.atof(items[4])
for imode in range(NofBranches):
linecount+=1
items=string.split(sl[linecount])
frequencies[newQn,imode] = string.atof(items[1])*energy_scaling
linecount+=2
for imode in range(NofBranches):
for idegree in range(0,NofBranches,3):
linecount+=1
items=map(string.atof,string.split(sl[linecount])[2:8])
eigenvectors[newQn,imode, idegree:idegree+3].real=items[0:8:2]
eigenvectors[newQn,imode, idegree:idegree+3].imag=items[1:8:2]
elif tipo in ["Vasp","phonopy"]:
stime=time.time()
print "----- reding yaml file ", filename
file = open(filename,"r")
yamlData = load(file, Loader=Loader)
print " reading Took ",etime - stime, "seconds"
NofIons= yamlData["natom"]
NofBranches = 3 * NofIons
NofWaves= len(yamlData["phonon"])
print "\tNofIons\t\tNofBranches\tNofWaves"
print "\t", NofIons, "\t\t",NofBranches, "\t\t",NofWaves
bvs = yamlData["reciprocal_lattice"]
cellVects =BrillVects(bvs)/2.0/math.pi
print "\tCell vectors "
for v in cellVects:
print "\t%e20.7 %e20.7 %e20.7"%tuple(v.tolist() )
atomReducedPositions = yamlData["atomReducedPositions"]
atomNames = yamlData["atomNames"]
atomMasses = numpy.array(yamlData["atomMasses"])
qs= numpy.zeros([NofWaves,3],"f")
weights = numpy.zeros([NofWaves],"f")
frequencies = numpy.zeros([NofWaves,NofBranches],"f")
eigenvectors = numpy.zeros([NofWaves,NofBranches, NofBranches],"F")
Phonons = yamlData["phonon"]
for iwav in range(NofWaves):
phonon = Phonons[iwav]
qs[iwav] =phonon["q-position"]
if numpy.sum( qs[iwav]*qs[iwav] )==0.0 :
qs[iwav]= [0.0,0.0,1.0e-6]
weights[iwav]= phonon["weight"]
for imode in range(NofBranches):
frequencies[iwav,imode] = phonon["band"][imode]["frequency"]*energy_scaling*33.35641
for imode in range(NofBranches):
vvv = numpy.array( phonon["band"][imode]["eigenvector"])
for idegree in range(0,NofBranches,3):
eigenvectors[iwav,imode, idegree:idegree+3].real = vvv[idegree/3][:,0]
eigenvectors[iwav,imode, idegree:idegree+3].imag = +vvv[idegree/3][:,1]
else:
raise Exception, " Unknown phonon file type: not Castep neither Vasp, phonopy"
# eigenvectors=numpy.array(numpy.swapaxes(eigenvectors, 1,2))
res=CalcDatas( NofIons= NofIons, NofBranches= NofBranches, NofWaves= NofWaves,
cellVects= cellVects, atomReducedPositions= atomReducedPositions,
atomNames= atomNames, atomMasses= atomMasses,
frequencies= frequencies, eigenvectors= eigenvectors,
weights= weights, QsReduced=qs,
namemd5 = namemd5)
qs = res.GetNonClippedQs()
print " GetNonClippedQs "
print qs
assert( len(qs) == len (eigenvectors) )
for i in range(len(qs)):
if tipo in ["Vasp","phonopy"] :
q=qs[i]
facts = numpy.tensordot( q , res.atomAbsolutePositions , axes= [(0),(1)])
facts=(numpy.exp(+1.0j*facts)).astype( numpy.complex64 )
facts = facts[:,None]*(numpy.array([1,1,1]).astype( numpy.complex64 ))
facts.shape = ( facts.shape[0]*3 , )
eigenvectors[i] = eigenvectors[i] *facts
# print qs[i]
# print eigenvectors[i]
print "----- Saving the readen structure to file ", namemd5
stime=time.time()
if(USE_HDF5):
print " ()() !! using HDF5"
res.dump_on_file(namemd5)
else:
print " ()() !! using cPickle"
f=open(namemd5,"w")
cPickle.dump(res,f)
f=None
etime = time.time()
print " Saving Took ",etime - stime, "seconds"
# end of parsing phonon file
## try to read preparsed data for Born charges and Polarisability from a hdf5 file (resExtra)
Extra_has_been_read = False
if filename_extra is None :
namemd5=namemd5 # the info might eventually be here. Such name has been set before.
else:
stime=time.time()
print "-----reading the file ", filename_extra
s=open( filename_extra,"r").read()
etime = time.time()
print " Reading Took ",etime - stime, "seconds"
# Read what has to be preprocessed. Result will then be (over)written to namemd5
if filename_extra is None and os.path.exists(namemd5):
print " ()()()() Looking for Born Charges and P in prestored file ", namemd5
s=None
print " ()() !! Now reloading the structure "
stime=time.time()
if(USE_HDF5):
print " ()() !! using HDF5"
Extra_has_been_read= res.read_extra_from_file(namemd5)
else:
raise Exception , " only hdf5 is supported "
etime = time.time()
print " Reloading Took ",etime - stime, "seconds"
# end of try from hdf5 for phonons
if not Extra_has_been_read and filename_extra is not None: # start parsing extra
print "----- reading Extra BornCharges and Polariz from ", filename_extra
stime=time.time()
if tipo == "Castep":
print "----- reading Extra BornCharges and Polariz from ", filename_extra
sl=string.split(s,"\n")
s=None
etime = time.time()
print " splitting Took ",etime - stime, "seconds"
linecount=0
while( "Polarisabilities" not in sl[linecount]):
linecount+=1
linecount+=2 # to jump over line :Optical (f->infinity) Static (f=0)
# and following
Polariz=[]
for ipol in range(3):
linecount+=1
Polariz.append(map(string.atof, string.split(sl[linecount])[:3] ))
Polariz=numpy.array(Polariz) # Polarizability is a symmetric tensor. No worry about transposition
# The Born effective charges are laid out with the columns representing
# the X,Y,Z electric field directions and the rows the X,Y,Z displacement directions.
while( "Born Effective Charges" not in sl[linecount]):
linecount+=1
linecount+=1
Borns=[]
for idisp in range(3* res.NofIons):
linecount+=1
if idisp%3 ==0:
Borns.append(map(string.atof, string.split(sl[linecount])[2:] ))
else:
Borns.append(map(string.atof, string.split(sl[linecount]) ))
elif tipo in ["Vasp","phonopy"] :
yamlData = load(s, Loader=Loader)
print " reading Took ",etime - stime, "seconds"
dielectric =numpy.array(yamlData["Dielectric constant"])
dielectric[0,0] -=1.0
dielectric[1,1] -=1.0
dielectric[2,2] -=1.0
Polariz = dielectric * abs(res.GetCellVolume()) / 4.0 / math.pi
Borns = yamlData["Born effective charges"]
res.BornCharges= numpy.array(Borns )
res.Polariz = numpy.array(Polariz)
print "----- Saving the readen Extra to file ", namemd5
res.dump_extra_on_file(namemd5)
return res