#/*##########################################################################
# 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 numpy
import copy
import numpy.linalg
from math import cos,sin
import math
import os
import time
import sys
import h5py
import TDS_Reading
import sys
import time
from TDS_Reading import SMALLQ
try:
from mpi4py import MPI
myrank = MPI.COMM_WORLD.Get_rank()
nprocs = MPI.COMM_WORLD.Get_size()
procnm = MPI.Get_processor_name()
comm = MPI.COMM_WORLD
print "MPI LOADED , nprocs = ", nprocs
except:
nprocs=1
myrank = 0
print "MPI NOT LOADED "
####################################################
# Returns 1 ,-1 or zero
# depending on the ciclicity of the a,b,c list.
# Utility function.
# If a,b,c can be reordered in a monotonic arrangement
# by an even number of permutations retuns 1
# If the number is odd returns -1
# If two numbers are equal retorns zero
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
###################################
# Utility function.
# Return simply the element-wise difference
# of two tuples : a-b
def difftuple(a,b):
res=()
for i in range(0,len(a)):
res= res+ (a[i]-b[i],)
return res
######################################""
# utility function.
# Return the intersection of two lists.
# Uses a dictionary to go faster
#
def intersection(list1, list2):
int_dict = {}
list1_dict = {}
for e in list1:
list1_dict[e] = 1
for e in list2:
if list1_dict.has_key(e):
int_dict[e] = 1
return int_dict.keys()
###########################################################
# This function test the equality between two OP_cella
# objects. The OP_cella is defined further
def OP_isDifferent(cella, cellb, verbose=0):
#################################
# First find the linear combination of axis that gives the shift.
# The cell shifts are decomposede on cella.vectors basis
A = numpy.transpose(cella.cellvectors)
A = numpy.linalg.inv (A )
coeffa = numpy.dot(A,cella.shift)
coeffb = numpy.dot(A,cellb.shift)
############################################################
# calculates the relative shift of cellb relative to cella
# in units of the cellb.vectors
cellb.unitsshift=numpy.zeros(3,numpy.float32)
for i in range(0,3):
cellb.unitsshift[i]= (( coeffb [i] + 10.0 ) % 1.0 ) -0.
# the differences are modulus 1. (periodicity)
if verbose:
print "coeffa ", coeffa
print "coeffb ", coeffb
print "cella.cellvectors ", cella.cellvectors
print "cellb.cellvectors ", cellb.cellvectors
unitsshift=numpy.zeros(3,numpy.float32)
for i in range(0,3):
unitsshift[i]= (((coeffb-coeffa)[i]+10.5) %1.0) -0.5
if(numpy.sum(abs(unitsshift))>0.1):
return 1 # they are different
diff=abs(cella.cellvectors-cellb.cellvectors)
diff=numpy.sum(numpy.sum(diff))
if(diff>0.1):
return 1
# if we arrive here, it is really the same cella
# So they are not different
return 0
# end of OP_isDifferent(cella, cellb)
############################################################
# class OP_cella
# This class stores things like atom names, their positions,
# cell axis, shift ......
# it has methods for rotation, reflection,
# to calculate Fourier transform, for finding the way one can rotate a set
# of atoms to another set of atoms, for traking interactions and so on......
class OP_cella:
#########################################################"
# Constructor of OP_cella.
# Here one passes the relevant parameter (see default list)
def __init__(self,
cellvectors=numpy.array( [ (1.0,0.0,0.0),(0.0,1.0,0.0), (0.0,0.0,1.0) ], "f" ),
AtomNames_long=None,
PositionsList_flat=None
):
if cellvectors is not None: # this happens when deepcopying
AtomNames={}
Names=[]
self.PositionsList_flat = numpy.array(PositionsList_flat)
for name in AtomNames_long:
if not name in Names:
Names.append(name)
AtomNames[name]=[]
for name, v in zip ( AtomNames_long, PositionsList_flat ):
AtomNames[name].append(v)
PositionsList=[]
for name in Names:
PositionsList.append(AtomNames[name])
AtomNames=Names
self.cellvectors=cellvectors
self.AtomNames =AtomNames
self.PositionsList =(PositionsList)
self.atomlist=[]
for i in range(0, len(AtomNames)):
for j in range(0, len(PositionsList[i])):
self.atomlist.append( (AtomNames[i],j ) )
# always called even if deepcopy
self.shift=numpy.array([0.0,0.0,0.0], "f")
self.unitsshift=numpy.zeros(3,numpy.float32)
self.reflection=1
####################################################
# OP_cella.Reflection(self)
# Reflection reflects all the position on the cell origin.
# The shift is left untouched because it is the shift of the origin itself.
def Reflection(self):
self.cellvectors=(-1.0)*self.cellvectors
self.reflection=-self.reflection
for i in range(0, len(self.PositionsList)):
self.PositionsList[i]=(-1.0)*self.PositionsList[i]
self.PositionsList_flat=(-1.0)*self.PositionsList_flat
##########################################################
# OP_cella.getK(self, N=2)
# Returns an array of 3D vectors, that are spanned by
# integer multiples of cell's Brillouin vectors.
# integers run form -N to N.
# (the vectors are already multiplied by imaginary unit
# in order to be ready for fourier transforming)
def getK(self, N=2):
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, self.cellvectors[2]),self.cellvectors[1] )
cellVolume=numpy.dot(self.cellvectors[0],cellVolume)
Brillvectors=numpy.zeros([3,3], numpy.float32)
Brillvectors[0]=2 * numpy.pi * numpy.dot( numpy.dot(AntiSymm, self.cellvectors[2]),self.cellvectors[1] )/ cellVolume
Brillvectors[1]=2 * numpy.pi * numpy.dot( numpy.dot(AntiSymm, self.cellvectors[0]),self.cellvectors[2] )/ cellVolume
Brillvectors[2]=2 * numpy.pi * numpy.dot( numpy.dot(AntiSymm, self.cellvectors[1]),self.cellvectors[0] )/ cellVolume
K=[ numpy.dot( [ i,j,k], Brillvectors ) for k in range(-N,N+1) for j in range (-N,N+1) for i in range(-N,N+1) ]
K=numpy.array(K)*complex(0,1.0)
self.Brillvectors=Brillvectors
self.volume= cellVolume
return K
def get_Correspondance_Table(self, cella, Rot):
numat=0
for positions in self.PositionsList:
for pos in positions:
numat=numat+1
totnumat=numat
res=numpy.zeros([numat*3,numat*3], "d" )
numat=0
# for positions in self.PositionsList:
# for pos in positions:
for pos in self.PositionsList_flat:
numat2=0
ok=0
# for positions2 in cella.PositionsList:
# if ok:
# break
# for pos2 in positions2:
for pos2 in cella.PositionsList_flat :
diff = pos + self.shift -pos2
diff = numpy.dot( cella.Brillvectors, diff )/(2*math.pi)
diff = numpy.sum(numpy.abs(((diff+100.00001)%1.0)-0.00001))
if diff<1.0e-5:
# parte con la rotazione e atterra su numat2ruotato
res[ numat2*3:(numat2+1)*3, numat*3:(numat+1)*3 ] = Rot
# res[ numat2*3:(numat2+1)*3, numat*3:(numat+1)*3 ] = self.Rotation.transpose()
ok=1
break
numat2=numat2+1
if numat2==totnumat:
raise " Correspondance not found "
numat=numat+1
return res
###########################################################
# OP_cella.getFourier(self, K)
# return an array (for each atom kind(element))
# of arrays (for each vector in K) of fourier component
# Dirac delta are put on the atoms locations.
def getFourier(self, K):
res=[]
for positions in self.PositionsList:
resAtom=numpy.zeros(len(K))
for pos in positions:
pos=pos+self.shift
resAtom=resAtom+numpy.exp(numpy.dot( K , pos ) )
res.append(resAtom)
return numpy.array(res)
################################################################################
# OP_cella.Transform(self,Az1, Ax, Az2 ,DPOS )
# Rotates by the three Eulerian angles Az1, Ax, Az2 around
# the x and z axis through the origin and shifts the origin by DPOS
def Transform(self,Az1, Ax, Az2 ,DPOS ):
rotAz1=numpy.array([ [cos(Az1) , -sin(Az1) , 0 ],
[sin(Az1) , cos(Az1) , 0 ],
[ 0 , 0 , 1 ]
]
)
rotAx=numpy.array([ [ 1 , 0 , 0 ],
[ 0 , cos(Ax) , -sin(Ax) ],
[ 0 , sin(Ax) ,cos(Ax) ]
]
)
rotAz2=numpy.array([ [cos(Az2) , -sin(Az2) , 0 ],
[sin(Az2) , cos(Az2) , 0 ],
[ 0 , 0 , 1 ]
]
)
rot=numpy.dot( rotAz2, numpy.dot(rotAx,rotAz1) )
self.TransformByMatrix( rot ,DPOS )
################################################################################
# OP_cella.TransformByMatrix(self, rot ,DPOS )
# Rotates by the rot matrix and shifts the origin by DPOS
def TransformByMatrix(self, rot ,DPOS ):
####################################################
# Please to notice that as the vectors
# are stored row by row in self.cellvectors
# we have to rotate a row at once.
# ... we cannot do numpy.dot(rot, self.cellvectors)
for i in range(0, len(self.cellvectors)):
self.cellvectors[i]=numpy.dot(rot, self.cellvectors[i])
for ka in range(0, len(self.AtomNames)):
L=self.PositionsList[ka]
for i in range(0,len(L)):
L[i]=numpy.dot(rot, L[i])
L=self.PositionsList_flat
for i in range(0,len(L)):
L[i]=numpy.dot(rot, L[i])
# stores the most recent rotation parameters
self.shift=DPOS
self.Rotation=rot
############################################################
# OP_cella.computeDistances(self, maxdist)
# It builds up a Big multiple dictionary, namely
#
# self.distances[Aa][Bb][kA][kB][(DX,DY,DZ)]
#
# where Aa and kA are the atom name ( for example "Si",..)
# and the position in the array of vectors in cell.PositionsList
# starting from 0. Same thing for Bb and kB.
# (DX,DY,DZ) are integer numbers specifying the cell containing (Bb,kB)
# relatively to (Aa,kA).
# DX,DY,DZ run from -N to N
def computeDistances(self, maxdist, N=2):
self.distances={}
for iAa in range(0,len(self.AtomNames)):
Aa=self.AtomNames[iAa]
D1=self.distances[Aa]={}
for iBb in range(0,len(self.AtomNames)):
Bb=self.AtomNames[iBb]
D2=D1[Bb]={}
for kA in range(0, len(self.PositionsList[iAa])):
D3=D2[kA]={}
P1= self.PositionsList[iAa][kA]
for kB in range(0, len(self.PositionsList[iBb])):
D4=D3[kB]={}
Pb_base=self.PositionsList[iBb][kB]
for BX in range(-N,N+1):
for BY in range(-N,N+1):
for BZ in range(-N,N+1):
P2= Pb_base +numpy.dot( [BX,BY,BZ] , self.cellvectors )
dist=math.sqrt( numpy.sum( (P2-P1)*(P2-P1) ) )
if(dist<maxdist):
D4[(BX,BY,BZ)]=[ dist , P1+self.shift,P2+self.shift,[] ]
##################################################################################
# OP_cella. FindTetragons(self, Latoms, CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)])
# This is the key routine to find good candidates to symmetry group.
# Latoms is a list of the kind [ (Aa,kA) ,(Bb,kB)....]
# ( i.e. couples formed by atom-name and position in the position list.
# Latoms must be formed of four atoms defining a non-degenerate tetraedron.
# A check is permormed on the non-degeneracy.
# The funtions finds all the possible equivalent tetraedrons (which have the same
# set of distances, and the same atom kinds)
# The function return the list of all these tetraedrons
def FindTetragons(self, Latoms, CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)], cellspan=3):
Positions= self.FindPositions(Latoms, CellsCoo=CellsCoo )
listindexes= [ self.AtomNames.index(Latoms[i][0]) for i in range(0,len(Latoms))]
vectors = Positions[1:4]-Positions[0]
#####
# check for non-degeneracy
det= numpy.linalg.det(vectors)
if(abs(det)<0.00001):
print " ERROR : DETERMINANT IS ", det, " PROBABLE DEGENERACY "
os.exit()
# CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0) ]
CooList=[(i,j,k) for i in range(-cellspan,cellspan+1) for j in range(-cellspan,cellspan+1) for k in range(-cellspan,cellspan+1) ]
(constraintNames,constraintTripod,constraintCircle)=self.FindConstraints(Latoms, CellsCoo=CellsCoo)
cN=constraintNames
res=[(Latoms,CellsCoo )]
for iA0 in range(0, len(self.PositionsList[listindexes[0]]) ):
Coo0=(0,0,0)
for iA1 in range(0, len(self.PositionsList[listindexes[1]]) ):
for Coo1 in CooList:
d= self.distances[cN[0]][cN[1]][iA0 ][iA1 ][difftuple(Coo1, Coo0)][0]
if ( abs(d-constraintTripod[0])>0.0001 ): continue
for iA2 in range(0, len(self.PositionsList[listindexes[2]]) ):
for Coo2 in CooList:
d= self.distances[cN[0]][cN[2]][iA0 ][iA2 ][difftuple(Coo2, Coo0)][0]
if ( abs(d-constraintTripod[1])>0.0001 ): continue
for iA3 in range(0, len(self.PositionsList[listindexes[3]]) ):
for Coo3 in CooList:
d= self.distances[cN[0]][cN[3]][iA0 ][iA3 ][difftuple(Coo3, Coo0)][0]
if ( abs(d-constraintTripod[2])>0.0001 ): continue
LatomsNew=[ (constraintNames[i] ,iA ) for i,iA in [(0,iA0),(1,iA1),(2,iA2),(3,iA3)] ]
(NamesNew,TripodNew,CircleNew)=self.FindConstraints(LatomsNew, [Coo0,Coo1,Coo2,Coo3] )
if( numpy.sum( abs(TripodNew-constraintTripod)+
abs(CircleNew-constraintCircle) )<0.001):
res.append((LatomsNew, [Coo0,Coo1,Coo2,Coo3] ))
PositionsNew= self.FindPositions( LatomsNew, [Coo0,Coo1,Coo2,Coo3] )
vectorsNew = PositionsNew[1:4]-PositionsNew[0]
return res
###############################################################################
# The next two funtion are auxilary funtions for
# Findtetragons
def FindPositions(self, Latoms, CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)]):
Positions=numpy.zeros([4,3], numpy.float32)
for i in range(0,4):
Positions[i]=self.PositionsList[self.AtomNames.index(Latoms[i][0])][ Latoms[i][1] ]
Positions[i]=Positions[i]+numpy.dot(CellsCoo[i],self.cellvectors)
return Positions
def FindConstraints(self, Latoms, CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)] ):
constraintNames=[Latoms[i][0] for i in range(0,4) ]
La=Latoms
constraintTripod=numpy.array([ self.distances[La[0][0]][La[i][0]][La[0][1] ][La[i][1] ][CellsCoo[i]][0]
for i in range(1,4) ] )
constraintCircle=numpy.array([ self.distances[La[i][0]][La[( i)%3 +1 ][0]][La[i][1] ]
[La[(i)%3 +1][1] ][difftuple(CellsCoo[i%3+1], CellsCoo[i])][0]
for i in range(1,4) ])
return (constraintNames,constraintTripod,constraintCircle)
def FindPosition(self, Latoms, CellCoo=(0,0,0)):
Position=self.PositionsList[self.AtomNames.index(Latoms[0])][ Latoms[1] ]
Position=Position +numpy.dot(CellCoo,self.cellvectors)
return Position
###############################################################
# OP_cella.FindEquivalences(self, CellsList)
# This function should receive a list of cells (CellList)
# that are superpose as they are, atom by atom to self.
# (program stops if this condition is not satisfied)
# It sets self.equivalences containing a list of lists containing set of equivalent atoms in the format [(Aa,kA),...],
# for each cell atom. It is exactly a multiplication table for the symmetry group
#
# self.grouplist containing the names given to the groups (like Cu_G1 )
# self.Groups a list of groups. Each group contains the equivalent atoms in the format (Aa,kA)
# self.DGroups self.DGroups[(Aa,kA)] gives the name (a name like Cu_G1 ) to which atom (Aa,kA) belongs
def FindEquivalences(self, CellsList):
equivalences={}
for i in range(0, len(self.AtomNames)):
name=self.AtomNames[i]
positions =self.PositionsList[i]
for iA in range(0, len(positions)):
EqAtoms=[]
equivalences[ (name,iA) ]= EqAtoms
for cell in CellsList:
eqa=cell.FindEquivalentOf(self.PositionsList[i][iA]+self.shift, name )
if(eqa==None):
print "BIG Problem"
print " Found no equivalent atom in a cella that should be equivalent"
os.exit()
EqAtoms.append(eqa)
self.equivalences= equivalences
Groups=[]
grouplist=[]
DGroups={}
for i in range(0, len(self.AtomNames)):
name=self.AtomNames[i]
NGr=[]
for iA in range(0, len(self.PositionsList[i] )):
Na= (name,iA)
ngr=0
for group in NGr:
if( Na in group):
break
if( intersection(self.equivalences[Na], group) !=[] ):
group.append(Na)
DGroups[Na]="%s_G%d" %(name,ngr)
break
ngr=ngr+1
else:
NGr.append([Na])
DGroups[Na]="%s_G%d" %(name,ngr)
grouplist.append("%s_G%d" %(name,ngr))
Groups=Groups+NGr
self.grouplist=grouplist
self.Groups=Groups
self.DGroups=DGroups
####################################################
# utility funtion used by FindEquivalences
def FindEquivalentOf(self, position, name ):
A = numpy.transpose(self.cellvectors)
A = numpy.linalg.inv (A )
positions =self.PositionsList[self.AtomNames.index(name)]
for iA in range(0, len(positions)):
pos=positions[iA]+self.shift
diff = numpy.dot(A, position-pos )
for i in range(0,3):
diff[i]= ((diff[i]+0.5) %1.0) -0.5
if(numpy.sum(abs(diff))<0.001):
return((name,iA))
return None
###################################################
# OP_cella.printEquivalences(self)
# prints in a more readeble format
# self.equivalences
# self.DGroups
def printEquivalences(self):
print "TABLE D'EQUIVALENCES"
for i in range(0, len(self.AtomNames)):
name=self.AtomNames[i]
positions=self.PositionsList[i]
for iA in range(0, len(positions)):
print "%2s%d -----> " %(name,iA),
EqAtoms=self.equivalences[ (name,iA) ]
count=0
for at in EqAtoms:
print " %2s%d " % at,
count=count+1
if(count%8 ==0):
print ""
print " ",
print ""
print ""
for group in self.Groups:
prefix=" "
print self.DGroups[group[0]]," : ",
for Na in group:
print "%s" % prefix, "%2s%d " % Na ,
prefix=" == "
print "\n"
################################################
# OP_cella.computeHistograms(self, DR=0.1)
# for all the couples of atom groups g1,g2, creates a multiple dictionary
# histo[g1][g2]
# giving a list of items in the format [ [ r1 ,r2 , Ninteractions)....
# where Ninteractions are the number of interactions between an atom of g1
# and one of g2 at a distance between r1 and r2, r1-r2 being the maximum interval smaller
# than DR
def computeHistograms(self, DR=0.1):
distances={}
where={}
histo={}
for g1 in self.grouplist:
distances[g1]={}
where[g1]={}
histo[g1]={}
for g2 in self.grouplist:
distances[g1][g2]=[]
where[g1][g2]={}
histo[g1][g2]=[]
for iAa in range(0,len(self.AtomNames)):
Aa=self.AtomNames[iAa]
for iBb in range(0,len(self.AtomNames)):
Bb=self.AtomNames[iBb]
for kA in range(0, len(self.PositionsList[iAa])):
for kB in range(0, len(self.PositionsList[iBb])):
g1=self.DGroups[(Aa,kA)]
g2=self.DGroups[(Bb,kB)]
lista=distances[g1][g2]
lw= where[g1][g2]
dictio=self.distances[Aa][Bb][kA][kB]
for Coo in dictio.keys():
lista.append(dictio[Coo][0])
lw[dictio[Coo][0]]=("like %s%d in 0 and %s%d in %d %d %d" % (Aa,kA,Bb,kB,Coo[0],Coo[1],Coo[2]))
for g1 in self.grouplist:
histo[g1]={}
for g2 in self.grouplist:
histo[g1][g2]=[]
distances[g1][g2] = numpy.sort(distances[g1][g2])
lw= where[g1][g2]
l= distances[g1][g2]
Rinf=-100
hi=histo[g1][g2]
tok=None
for r in l:
if(r<0.0001): continue
if(r-Rinf>DR):
if(tok!=None): hi.append(tok)
Rinf=r
tok=[Rinf,r,1,lw[ r ]]
else:
tok[1]=r
tok[2]=tok[2]+1
hi.append(tok)
self.histo=histo
for g1 in self.grouplist:
for g2 in self.grouplist:
for h in histo[g1][g2]:
h[2]=h[2]/2
####################################
# print self.histo
#
def printHisto(self,n):
for i in range(0,len(self.grouplist)):
for j in range(i ,len(self.grouplist)):
g1=self.grouplist[i]
g2=self.grouplist[j]
print g1, "---",g2," ",
h=self.histo[g1][g2]
count=0
for t in h:
print " %d in [%4f,%4f];" %(t[2],t[0],t[1]),
count+=1
if(count==n): break
print " "
def printParametersInput(self,filename, n=2):
fd=open(filename,"w")
fd.write("# This first part is just an initialisation\n")
fd.write("# Edit below \n")
fd.write("M_={}\nZ_={}\nY_={}\nK_={}\nVWs_={}\nLJs_={}\nBMs_={}\n")
fd.write("BK_L={}\nSM_L={}\nBK_T={}\nSM_T={}\nTens_P={}\nTens_S={}\nTens_T={}\nTens_YT={}\nangleBonds={}\nFunction_YT=None\nSCREENING_FUNCTION=None\n")
for i in range(0,len(self.grouplist)):
g1=self.grouplist[i]
fd.write("BK_L['%s']={}\nSM_L['%s']={}\nBK_T['%s']={}\nSM_T['%s']={}\nTens_P['%s']={}\nTens_S['%s']={}\nTens_T['%s']={}\nTens_YT['%s']={}\n"%( (g1,)*8 ))
for j in range(0 ,len(self.grouplist)):
g2=self.grouplist[j]
fd.write("SM_L['%s']['%s']={}\nSM_T['%s']['%s']={}\n"%( (g1,g2)*2 ))
fd.write("# HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH\n")
fd.write("# Edit the followin part\n")
fd.write("# Masses : M \n")
for i in range(0,len(self.grouplist)):
fd.write("M_['%s']=0.0\n" % self.grouplist[i])
fd.write("# Coulomb interaction : Z \n")
for i in range(0,len(self.grouplist)):
fd.write("Z_['%s']=0.0\n" % self.grouplist[i])
fd.write("\n# Shells Model : shell charge \n")
for i in range(0,len(self.grouplist)):
fd.write("Y_['%s']=0.0\n" % self.grouplist[i])
fd.write("\n# Shells Model : polarizability \n")
for i in range(0,len(self.grouplist)):
fd.write("K_['%s']=1000000.0\n" % self.grouplist[i])
fd.write("\n# Van der Walls : sigma \n")
fd.write("VW_V0=0.0\n")
for i in range(0,len(self.grouplist)):
fd.write("VWs_['%s']=1.0\n" % self.grouplist[i])
fd.write("\n# Lenard Jones : sigma \n")
fd.write("LJ_V0=0.0\n")
for i in range(0,len(self.grouplist)):
fd.write("LJs_['%s']=1.0\n" % self.grouplist[i])
fd.write("\n# Born-Mayer : sigma \n")
fd.write("BM_V0=0.0\n")
for i in range(0,len(self.grouplist)):
fd.write("BMs_['%s']=1.0\n" % self.grouplist[i])
for i in range(0,len(self.grouplist)):
g1=self.grouplist[i]
for j in range(i ,len(self.grouplist)):
g2=self.grouplist[j]
h=self.histo[g1][g2]
fd.write( "# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX \n" )
fd.write( "# Interactions between '%s' and '%s' \n" % (g1,g2) )
count=0
for t in h:
if(count==n): break
fd.write( "# Shell N%d going from %f to %f %s \n" % (count,t[0],t[1],t[3]))
count+=1
fd.write("BK_L['%s']['%s']=[]\n" %(g1,g2) )
fd.write("BK_T['%s']['%s']=[]\n" %(g1,g2) )
fd.write("SM_L['%s']['%s']['SS']=[]\n" %(g1,g2) )
fd.write("SM_T['%s']['%s']['SS']=[]\n" %(g1,g2) )
fd.write("SM_L['%s']['%s']['SC']=[]\n" %(g1,g2) )
fd.write("SM_T['%s']['%s']['SC']=[]\n" %(g1,g2) )
fd.write("SM_L['%s']['%s']['CS']=[]\n" %(g1,g2) )
fd.write("SM_T['%s']['%s']['CS']=[]\n" %(g1,g2) )
if(g1!=g2):
fd.write("# just for symmetry \n")
fd.write("SM_L['%s']['%s']=SM_L['%s']['%s']\n" %(g2,g1,g1,g2) )
fd.write("SM_T['%s']['%s']=SM_T['%s']['%s']\n" %(g2,g1,g1,g2) )
fd.write("BK_L['%s']['%s']=BK_L['%s']['%s']\n" %(g2,g1,g1,g2) )
fd.write("BK_T['%s']['%s']=BK_T['%s']['%s']\n" %(g2,g1,g1,g2) )
for letter in ["P","T","S","YT"]:
for i in range(0,len(self.grouplist)):
g1=self.grouplist[i]
for j in range(i ,len(self.grouplist)):
g2=self.grouplist[j]
h=self.histo[g1][g2]
fd.write( "# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX \n" )
fd.write( "# Tensorial Interactions between '%s' and '%s' \n" % (g1,g2) )
count=0
for t in h:
if(count==n): break
fd.write( "# Shell N%d going from %f to %f %s \n" % (count,t[0],t[1],t[3]))
count+=1
fd.write("Tens_%s['%s']['%s']=[]\n" %(letter, g1,g2) )
if(g1!=g2):
fd.write("# just for symmetry \n")
fd.write("Tens_%s['%s']['%s']=Tens_%s['%s']['%s']\n" %(letter,g2,g1,letter,g1,g2) )
def __deepcopy__(self, memo):
res=OP_cella(cellvectors=None)
res.cellvectors=numpy.array(self.cellvectors , copy=1 )
res.AtomNames =copy.copy(self.AtomNames )
res.PositionsList = copy.copy(self.PositionsList)
for i in range(0, len(res.PositionsList)):
res.PositionsList[i]=numpy.array(res.PositionsList[i], copy=1 )
res.PositionsList_flat = numpy.array(self.PositionsList_flat, copy=1 )
res.shift=self.shift
return res
def OP_TellAbout(newCell, short="YES"):
# findout the type of rotation
if(isinstance(newCell, OP_cella)):
A=newCell.Rotation *complex(1.0,0.0)
else:
A=newCell*complex(1.0,0.0)
(eval, evect)=numpy.linalg.eig(A)
if(short=="NO"):
print "Rotazione ="
print A
print " evect = "
print evect
print " eval="
print eval
# eigenvalue 1 corresponds to rotation axis
# Find rotation axis
for i in range(0,3):
if (abs(1-eval[i])<0.0001):
rotation_axis = evect[i]
# Now we are ready to calculate the rotation
# angle. We chose a different eigenvalue/vector
# than the fixed axis. The argument of the eigenvalue
# gives us the rotation angle. Its sign is calculated
# considering the eigenvector : the real and imaginary part
# are two orthogonal vectors. Their vector product times
# the fixed axis gives the sign.
othervector = evect[ :, (i+1) %3 ]
otherval = eval [(i+1) %3 ]
angle = math.atan2(otherval.imag, otherval.real )
v1 = othervector.real*math.sqrt(2.0)
v2 = othervector.imag*math.sqrt(2.0)
triedron = numpy.linalg.det(numpy.array([rotation_axis,v1,v2]))
if(triedron.real>0):
angle=-angle
break
else:
print " something wrong finding fixed point of rotation "
sys.exit()
s=str( rotation_axis.real)
s=s+" "*(40-len(s))
s=s+str( 180*angle/math.pi )
s=s+" "*(60-len(s))
if(isinstance(newCell, OP_cella)):
# print " Shift = ", newCell.shift
s=s+str( newCell.unitsshift)
if(newCell.reflection==-1):
s=s+ " YES "
else:
s=s+ " NO "
return s
def Talk_about_Distance_variation(evects,one_over_sqrtM,cella, Neigen):
""" Here evect should be already multiplied by the right phase
"""
evects =evects*one_over_sqrtM
for g1 in cella.grouplist :
for g2 in cella.grouplist :
max_merit=-1.0e-50
max_blabla=""
for a1 in cella.Groups[cella.grouplist.index(g1)]:
for a2 in cella.Groups[cella.grouplist.index(g2)]:
if(a1==a2):
continue
index_of_a1_incellalist=cella.atomlist.index(a1)
index_of_a2_incellalist=cella.atomlist.index(a2)
vect3D_a1= evects[Neigen][ index_of_a1_incellalist*3: index_of_a1_incellalist*3+3]
vect3D_a2= evects[Neigen][ index_of_a2_incellalist*3: index_of_a2_incellalist*3+3]
for Coo in [(0,0,0) , (1,1,1), (0,1,1), (1,0,1), (1,1,0) , (1,0,0),(0,1,0), (0,0,1), ]:
DR=cella.FindPosition(a2, CellCoo=Coo ) - cella.FindPosition(a1, CellCoo= (0,0,0))
dDR_real=numpy.sum( DR*(( vect3D_a2- vect3D_a1 ).real) )
dDR_imag=numpy.sum( DR*(( vect3D_a2- vect3D_a1 ).imag) )
dDR = math.sqrt(dDR_real*dDR_real+dDR_imag*dDR_imag)
dist= math.sqrt( numpy.sum(DR*DR))
merit = dDR / numpy.sum(DR*DR)
if(abs(merit)>max_merit):
max_merit=abs(merit)
max_blabla=" dDR=%e for a1=(%s,%d) and a2=(%s,%d) in (%d,%d,%d) at a distance of %e" % ( (dDR,)+a1+a2+Coo+(dist,))
print "for g1=%s and g2=%s the blablabla is %s\n" %( g1,g2, max_blabla)
def simmetrizzatransl(translation2symmetrize, RotList, equivalences , AtomNames):
result=copy.deepcopy(translation2symmetrize)
for i in range(0, len(AtomNames)):
result[i]=result[i]*0
for i in range(0, len(AtomNames)):
name=AtomNames[i]
shifts = translation2symmetrize[i]
for iA in range(0, len(shifts)):
shift = shifts[iA]
eqs = equivalences[ (name, iA) ]
assert(len(RotList)==len(eqs) )
print eqs
# la tavola delle equivalenze e viziosa. Se P sta nella cella non
# girata, e P' in quella girata da rot : P est equivalente a P'
# se ci cade sopra
for irot in range(len(RotList)):
eqpos = eqs[ irot][1]
newshift = numpy.dot( RotList[irot],shifts[eqpos] )
result[i][iA] += newshift
for i in range(0, len(AtomNames)):
result[i]=result[i]/len(RotList)
return result
def fillDicWithGroup(res, group):
def setta(a,b,d=res):
res[a]=b[:]
group.visititems(setta)
def FindSimmetries(cella, filename= None,
md5postfix=None, overwrite= False,
key = None ):
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r')
if( key in h5 and not overwrite):
res={}
fillDicWithGroup(res,h5[key] )
h5.close()
return res
h5.close()
h5=None
print " ---------- Simmetry study :creating the dictionary with distances "
stime=time.time()
cella.computeDistances(200.0,N=2)
print " ---------- Simmetry study : dictionary created in ", time.time()-stime, " seconds "
##################################################################################
# OP_cella. FindTetragons(self, Latoms, CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)])
# This is the key routine to find good candidates to symmetry group.
# Latoms is a list of the kind [ (Aa,kA) ,(Bb,kB)....]
# ( i.e. couples formed by atom-name and position in the position list.
# Latoms must be formed of four atoms defining a non-degenerate tetraedron.
# A check is permormed on the non-degeneracy.
# The funtions finds all the possible equivalent tetraedrons (which have the same
# set of distances, and the same atom kinds)
# The function return the list of all these tetraedrons
#
# Every tetraedron is represented by a list of four atoms and a list of cell coordinates.... like
# [ [ [(Aa,kA),(Aa2,kA2),(Aa3,kA3), (Aa4,kA4)], [(0,0,0),(0,0,0),(0,0,0),(0,0,0)] ]
# [ .............................................. ]
# ]
#LTetra = cella.FindTetragons(
#[('Cu',0),('O',0),('La',2),('O',1)]
#)
SeekedTetragon=[(cella.AtomNames[0],0),(cella.AtomNames[0],0),(cella.AtomNames[0],0),(cella.AtomNames[0],0), ]
CellsCoo=[(0,0,0),(1,0,0),(0,1,0),(0,0,1)]
print " ---------- Simmetry study :giveen a tetraedron find the equivalent ones "
stime=time.time()
LTetra = cella.FindTetragons(SeekedTetragon, CellsCoo, cellspan=1)
print " ---------- Simmetry study : found the equivalent ones in ", time.time()-stime, " seconds "
#########################################################################################
# OP_cella.FindPositions(self, Latoms, CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)])
# returns an array of four vectors.
# The nth vector gives the position of nth atom of Latoms with the nth cell-coordinates
# of CellsCoo.
# Here, LTetra[0] which is the first tetraedron of the list, correspond by construction
# to identity
Positions= cella.FindPositions(LTetra[0][0],LTetra[0][1])
Vectors = Positions[1:4]-Positions[0]
##########################################################
# OP_cella.getK(self, N=2)
# Returns an array of 3D vectors, that are spanned by
# integer multiples of cell's Brillouin vectors.
# integers run form -N to N.
# (the vectors are already multiplied by imaginary unit
# in order to be ready for fourier transforming)
K = cella.getK()
Fourier = cella.getFourier(K)
##############################################################
# in CellsList we put transformed Cells that superpose
# perfectly to the original one through simmetry operations
# that are looked for in LTetra
#
CellsList=[]
##############################################################
# in RotList we put rotations for the transformed Cells that superpose
# perfectly to the original one through simmetry operations
# that are looked for in LTetra
#
RotList=[]
ShiftList=[]
##############################
# different rotations
DiffRotList=[]
#################################################"
# Start the loop on candidates to symmetry group
#
## print """
## *** Now printing all the simmetry operation of the crystal. ***
## The operation is described by a rotation centered on the origin, followed
## eventually by a reflection ( centered on the origin ). All this followed
## by a translation ( shift ).
## ROTATION AXIS | rotation angle | Reflection | SHIFT( in cell vects units )
## ----------------------------------------------------------------------------
## """
## who is going where
wgw_list=[]
for i in range(0, len(LTetra)):
NewPositions= cella.FindPositions(LTetra[i][0],LTetra[i][1])
NewVectors = NewPositions[1:4]-NewPositions[0]
Rot= numpy.dot(numpy.transpose(NewVectors), numpy.linalg.inv(numpy.transpose(Vectors)) )
shift=NewPositions[0]-numpy.dot(Rot,Positions[0])
if(numpy.linalg.det(Rot)<0):
Reflection=1
Rot=-Rot
else:
Reflection=0
newCell=copy.deepcopy(cella)
if(Reflection):
newCell.Reflection()
newCell.TransformByMatrix( Rot ,shift )
newFourier=newCell.getFourier(K)
diffFourier=newFourier-Fourier
diffFourier= (diffFourier* numpy.conjugate(diffFourier)).real
error=numpy.sum(numpy.sum(diffFourier))
if(error<0.0001):
isnew=1
for tok in CellsList:
if ( not OP_isDifferent(tok,newCell) ):
isnew=0
if(isnew):
## if(len(CellsList)):
## OP_isDifferent(CellsList[-1],newCell, verbose=1)
s=OP_TellAbout(Rot, short="YES")
s=s+ ["NO","YES"][Reflection]
s=s+" "*(70-len(s))
s=s+str( newCell.unitsshift )
CellsList.append(newCell)
newRot= Rot*(1-2*Reflection)
RotList.append(numpy.array(newRot))
ShiftList.append(numpy.array(newCell.unitsshift))
isnew=1
for tok in DiffRotList:
if ( sum(sum(abs(newRot-tok))) <0.001):
isnew=0
if isnew:
DiffRotList.append(newRot)
wgw_list.append( newCell.get_Correspondance_Table(cella , newRot))
# print s
res={}
res["RotList"]= numpy.array(RotList )
res["ShiftList"]= numpy.array(ShiftList )
res["DiffRotList"]= numpy.array(DiffRotList)
res["wgw"]= numpy.array(wgw_list)
## res["CellsList"] = numpy.array(CellsList )
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
for keyd in res.keys():
nuova_key= key+"/"+keyd
if(nuova_key in h5):
del h5[nuova_key]
h5[nuova_key]= res[keyd]
h5.flush()
h5.close()
return res
########################################################
# Now CellsList contains all the symmetric versions
# of the original cells
if __name__=="__main__":
res= TDS_Reading.ReadCastep(sys.argv[1])
cella=OP_cella(cellvectors=res.cellVects ,
AtomNames_long=res.atomNames,
PositionsList_flat=res.atomAbsolutePositions )
simmetries_dict = FindSimmetries(cella)
## return value of FindSimmetries
## res["RotList"]=RotList
## res["ShiftList"]=ShiftList
## res["diffRotList"]=diffRotList
## res["CellsList"] = CellsList
print " OBTAINED SIMMETRY TRANSFORMATIONS calling FindSimmetries "
print "NOW PRINTING THEM "
print """
*** Now printing all the simmetry operation of the crystal. ***
The operation is described by a rotation centered on the origin, followed
eventually by a reflection ( centered on the origin ). All this followed
by a translation ( shift ).
ROTATION AXIS | rotation angle | Reflection | SHIFT( in cell vects units )
----------------------------------------------------------------------------
"""
for Rot, shift in zip(simmetries_dict["RotList"], simmetries_dict["ShiftList"] ):
if numpy.linalg.det(Rot)<0.0 :
Reflection=1
Rot = -Rot
else:
Reflection=0
s=OP_TellAbout(Rot, short="YES")
s=s+ ["NO","YES"][Reflection]
s=s+" "*(70-len(s))
s=s+str( shift )
print s
print """
*** Now printing all the simmetry operation of the crystal. ***
The operation is described by a rotation centered on the origin, followed
eventually by a reflection ( centered on the origin ).
ROTATION AXIS | rotation angle | Reflection
-----------------------------------------------------------------
"""
for Rot in simmetries_dict["DiffRotList"]:
if numpy.linalg.det(Rot)<0.0 :
Reflection=1
Rot = -Rot
else:
Reflection=0
s=OP_TellAbout(Rot, short="YES")
s=s+ ["NO","YES"][Reflection]
s=s+" "*(70-len(s))
print s
def printSimmetries( RotList , ShiftList , QQreduced=None, Brills =None , Bravais= None ):
print " OBTAINED SIMMETRY TRANSFORMATIONS calling FindSimmetries "
print "NOW PRINTING THEM "
print """
*** Now printing all the simmetry operation of the crystal. ***
The operation is described by a rotation centered on the origin, followed
eventually by a reflection ( centered on the origin ). All this followed
by a translation ( shift ).
ROTATION AXIS | rotation angle | Reflection | SHIFT( in cell vects units )
---------------------------------------------------------------------------- """
if QQreduced is not None:
QQ=numpy.dot( QQreduced, Brills)
else:
QQ=None
for Rot, shift in zip(RotList, ShiftList ):
if QQ is not None:
QQp = numpy.dot(Rot,QQ) -QQ
scals = numpy.dot(Bravais, QQp )
scals =abs((scals+0.0001) % (2*math.pi))
scals = numpy.sum(numpy.less(0.001, scals))
if scals >0 : continue
if numpy.linalg.det(Rot)<0.0 :
Reflection=1
Rot = -Rot
else:
Reflection=0
s=OP_TellAbout(Rot, short="YES")
s=s+ ["NO","YES"][Reflection]
s=s+" "*(70-len(s))
s=s+str( shift )
print s
def printRotations( Rotations):
print """
*** Now printing all the different rotations that can be applied to Qs vectors. ***
The operation is described by a rotation centered on the origin, followed
eventually by a reflection ( centered on the origin ).
ROTATION AXIS | rotation angle | Reflection
-----------------------------------------------------------------"""
for Rot in Rotations:
if numpy.linalg.det(Rot)<0.0 :
Reflection=1
Rot = -Rot
else:
Reflection=0
s=OP_TellAbout(Rot, short="YES")
s=s+ ["NO","YES"][Reflection]
s=s+" "*(70-len(s))
print s
def GetCodes ( calculatedDatas, simmetries_dict , filename= None,
md5postfix=None,
overwrite= False,
key = None ,
APPLYTIMEREVERSAL=1
):
Qs = calculatedDatas.GetClippedQs()
print " Qs.shape ", Qs.shape
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r')
if( key in h5 and not overwrite):
return (h5[key+"/totQs"][:],h5[key+"/totIrot"][:], h5[key+"/totItime"][:],
h5[key+"/totBVifacts"][:], h5[key+"/steps"][:], h5[key+"/offsets"][:])
h5.close()
h5=None
reducedQs = calculatedDatas.Reduce(Qs)
offsets, steps = getOffsetSteps(reducedQs)
print "########offsets######### "
print offsets
print "----steps --------- "
print steps
print "XXXXXXXXXXXXXXXXXXXXXXX "
lengths=numpy.array([1,1,1.0])+ 6*steps
ilengths = (lengths/steps).astype(numpy.int32)
codesSymm=numpy.zeros(ilengths.tolist(), numpy.int64)
codesNq =numpy.zeros(ilengths.tolist(), numpy.int64)
codesNq[:]=-1
centerposition = ilengths/2
wheretogo = (reducedQs-offsets)/steps
iwheretogo = (1000.0+wheretogo+0.5).astype(numpy.int32)-1000
if max(abs( iwheretogo- wheretogo ).sum( axis=1)) >1.0e-3 :
raise Exception, " bad distribution of points "
codesSymm[ ((iwheretogo+centerposition ).transpose()).tolist() ] = 1
codesNq [ ((iwheretogo+centerposition ).transpose() ).tolist() ] = numpy.arange( len(iwheretogo) )
baricenter = numpy.sum(iwheretogo , axis=0 )*1.0/( len(iwheretogo))
print " baricenter " , baricenter+centerposition
diffs = (abs(iwheretogo-baricenter))
print " diffs "
print diffs
diffs = (abs(iwheretogo-baricenter)).sum(axis=1)
mindiff=min(diffs)
posmin=(numpy.nonzero( diffs==mindiff ))[0][0]
bari = iwheretogo[posmin] + centerposition
for b1,b2,b3 in (iwheretogo+ centerposition):
print b1,b2,b3
print " " , codesSymm[ b1+1,b2,b3] , codesSymm[ b1,b2+1,b3], codesSymm[ b1,b2,b3+1]
if( codesSymm[ b1+1,b2,b3]==1 and codesSymm[ b1,b2+1,b3]==1 and codesSymm[ b1,b2,b3+1]==1 ) :
print " it looks like a cubic lattice "
qgridtype="cubic"
break
elif ( codesSymm[ b1+2,b2,b3]==1 and codesSymm[ b1,b2+2,b3]==1 and codesSymm[ b1,b2,b3+2]==1 and
codesSymm[ b1+1,b2,b3]==0 and codesSymm[ b1,b2+1,b3]==0 and codesSymm[ b1,b2,b3+1]==0 and
codesSymm[ b1+1,b2+1,b3+1]==1 ):
print " it looks like a bcc lattice "
qgridtype="bcc"
break
else:
raise Exception," this qgridtype is still unknown "
###############################33
irot=0
symm_list =simmetries_dict["DiffRotList"]
codesSymm[:] =0
codesNq[:] =0
if APPLYTIMEREVERSAL:
segni=[-1,1]
else:
segni=[1]
Nsymm = len(symm_list)
totQs = numpy.zeros( [ Nsymm * len(Qs ) *len(segni) ,3], numpy.float32 )
totIrot = numpy.ones( [ Nsymm * len(Qs ) *len(segni) ] , numpy.int32 )
totItime = numpy.ones( [ Nsymm * len(Qs ) *len(segni) ] , numpy.int32 )
totBVifacts = numpy.zeros( [ Nsymm * len(Qs ) *len(segni) ,3], numpy.int32 )
isymm=0
BV=calculatedDatas.GetBrillVects()
for BVifacts in [(0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,1,0),(1,0,1),(0,1,1),(1,1,1)][:1] :
SHIFT = numpy.dot( BVifacts, BV )
for Rot in symm_list[::-1] : # per finire con l'identita
for segno in segni :
irot=irot+1
repQs = numpy.dot( Rot, Qs.T ).T
totQs [ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = repQs
totIrot [ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = irot
totItime[ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = segno
totBVifacts[ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = BVifacts
isymm=isymm+1
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
if(key in h5):
del h5[key]
print h5.keys()
h5 .create_group(key)
h5[key+"/totQs"] = totQs
h5[key+"/totIrot"] = totIrot
h5[key+"/totItime"]= totItime
h5[key+"/totBVifacts"]= totBVifacts
h5[key+"/steps"]=steps
h5[key+"/offsets"]= offsets
h5.flush()
h5.close()
h5=None
return totQs,totIrot, totItime, totBVifacts, steps, offsets
## reducedRepQs = calculatedDatas.Reduce(repQs) *segno
## x={}
## y={}
## z={}
## for tok in reducedQs:
## x[tok[0]-offsets[0]]=1
## y[tok[1]-offsets[1]]=1
## z[tok[2]-offsets[2]]=1
## ritorna qui
## ## print x.keys()
## ## print steps
## ## print y.keys()
## ## print steps
## ## print z.keys()
## wheretogo = (reducedRepQs-offsets)/steps
## iwheretogo = (1000.0+wheretogo+0.5).astype(numpy.int32)-1000
## ## for tok, tok1 in zip(wheretogo,iwheretogo ) :
## ## print tok, tok1
## ## ## if max(abs( tok- tok1 ) ) >1.0e-3 :
## ## ## raise Exception, " bad distribution of points "
## if max(abs( iwheretogo- wheretogo ).sum( axis=1)) >1.0e-3 :
## raise Exception, " bad distribution of points "
## codesSymm[((iwheretogo+centerposition ).transpose()).tolist()]=irot *segno
## codesNq[ ((iwheretogo+centerposition ).transpose()).tolist() ]=numpy.arange( len(iwheretogo) )
## result={}
## result["codesSymm"] = codesSymm
## result["codesNq"] = codesNq
## result["offsets"] = offsets
## result["steps"] = steps
## result["centerposition"] =centerposition
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
h5[key]=result
h5["qgridtype"]=qgridtype
h5.flush()
h5.close()
h5=None
return result, qgridtype
def compose_dm ( evals, orig_evects , wgw , positions, qshifts, acustic=False ):
shape_evects= orig_evects.shape
res=numpy.zeros( shape_evects , numpy.complex64)
if qshifts is not None:
facts = numpy.tensordot(qshifts , positions, axes= [(1),(1)])
facts=(numpy.exp(-1.0j*facts)).astype( numpy.complex64 )
vshape=orig_evects.shape
orig_evects.shape=tuple( list(vshape[:-1] ) + [ vshape[-1]/3, 3 ] )
vects= orig_evects * facts[:,None,:,None]
vects.shape = vshape
orig_evects.shape = vshape
else:
vects = orig_evects
vects = numpy.tensordot(vects , wgw , axes=([2], [ 1]) )
if not acustic:
N = shape_evects[1]
else:
N = 3
for i in range(N):
v = numpy.array(vects[:,i, :])
v= v* (evals[:,i,None])
add = v[:,:,None] * ((v.conjugate())[:,None,:])
numpy.add(res,add , res)
return res
# N = shape_evects[1]
# for i in range(N):
# print i
# ev , v = evals[:,i], orig_evects[:,i]
# v=(numpy.dot(wgw,v.transpose())*ev).transpose()
# add = v[:,:,None] * ((v.conjugate())[:,None,:])
# numpy.add(res,add , res)
# return res
def GetDM_fromFFT_shifted(Four_Rvects, Four_DMs, calculatedDatas, subN1, subN2, iN3, subN3_global ) :
nf1 = Four_DMs.shape[0]/2
nf2 = Four_DMs.shape[1]/2
nf3 = Four_DMs.shape[2]/2
subN3 = Four_DMs.shape[2]
nyq1 = subN1/2
nyq2 = subN2/2
nyq3 = subN3/2
paddedComps = numpy.zeros([subN1, subN2, subN3],"F")
nbranches = Four_DMs.shape[-1]
interp_DMs = numpy.zeros([subN1, subN2, 1,nbranches,nbranches ],"F")
for ib1 in range(nbranches):
for ib2 in range(nbranches):
paddedComps[ nyq1-nf1:nyq1+nf1+1 , nyq2-nf2:nyq2+nf2+1 , nyq3-nf3:nyq3+nf3+1 ] = Four_DMs[:,:,:, ib1, ib2]*(subN1*subN2*subN3)
shifted_paddedComps = numpy.fft.ifftshift( paddedComps ) * (numpy.exp( numpy.fft.fftfreq(subN3)*(subN3*iN3*1.0/subN3_global )*(2j)*math.pi )[None,None,:])
interp_DMs[:,:,:, ib1, ib2 ]=numpy.fft.ifftn( shifted_paddedComps )[:,:,0:1]
BV=calculatedDatas.GetBrillVects()
Qs_n1 = numpy.array( numpy.arange(subN1) [:, None]* BV[0]/subN1, "f")
Qs_n2 = numpy.array( numpy.arange(subN2) [:, None]* BV[1]/subN2, "f")
Qs_n3 = numpy.array( numpy.arange(iN3,iN3+1)[:, None]* BV[2]/subN3_global, "f")
Qs_nnn = Qs_n1[:,None,None,: ] + Qs_n2[None,:,None,: ] + Qs_n3[None,None,:,: ]
Qs_nnn[0,0,0]=Qs_nnn[0,0,0]+numpy.array([0.0001,0,0],"f")
if hasattr(calculatedDatas,"BornCharges"):
interp_DMs = interp_DMs + calculatedDatas.GetNonAnalytic( Qs_nnn )
return interp_DMs
def GetDM_fromFFT(Four_Rvects, Four_DMs, calculatedDatas, subN1, subN2, subN3 ) :
paddedComps = numpy.zeros([subN1, subN2, subN3],"F")
nbranches = Four_DMs.shape[-1]
interp_DMs = numpy.zeros([subN1, subN2, subN3,nbranches,nbranches ],"F")
nyq1 = subN1/2
nyq2 = subN2/2
nyq3 = subN3/2
# print " Four_DMs.shape " , Four_DMs.shape
# print " subN1, subN2, subN3 " ,subN1, subN2, subN3
nf1 = Four_DMs.shape[0]/2
nf2 = Four_DMs.shape[1]/2
nf3 = Four_DMs.shape[2]/2
# print paddedComps[ nyq1-nf1:nyq1+nf1+1 , nyq2-nf2:nyq2+nf2+1 , nyq3-nf3:nyq3+nf3+1 ].shape
# print Four_DMs[:,:,:, 0, 0].shape
for ib1 in range(nbranches):
for ib2 in range(nbranches):
# print " frequenza centrale QUI ", ( Four_DMs[nf1,nf2,nf3, ib1, ib2] )
# print " frequenza centrale QUI ", numpy.fft.ifftshift( Four_DMs[:,:,:, ib1, ib2] )[0,0,0]
paddedComps[ nyq1-nf1:nyq1+nf1+1 , nyq2-nf2:nyq2+nf2+1 , nyq3-nf3:nyq3+nf3+1 ] = Four_DMs[:,:,:, ib1, ib2]*(subN1*subN2*subN3)
shifted_paddedComps = numpy.fft.ifftshift( paddedComps )
# print " frequenza centrale QUA ", shifted_paddedComps [0,0,0]
interp_DMs[:,:,:, ib1, ib2 ]=numpy.fft.ifftn( shifted_paddedComps )
#### NO interp_DMs[:,:,:, ib2, ib1 ]=(interp_DMs[:,:,:, ib1, ib2 ]).conjugate()
# print " ritorno shape " , interp_DMs.shape
BV=calculatedDatas.GetBrillVects()
Qs_n1 = numpy.array( numpy.arange(subN1)[:, None]* BV[0]/subN1, "f")
Qs_n2 = numpy.array( numpy.arange(subN2)[:, None]* BV[1]/subN2, "f")
Qs_n3 = numpy.array( numpy.arange(subN3)[:, None]* BV[2]/subN3, "f")
Qs_nnn = Qs_n1[:,None,None,: ] + Qs_n2[None,:,None,: ] + Qs_n3[None,None,:,: ]
Qs_nnn[0,0,0]=Qs_nnn[0,0,0]+numpy.array([0.0001,0,0],"f")
interp_DMs = interp_DMs + calculatedDatas.GetNonAnalytic( Qs_nnn )
return interp_DMs
# def GetDM_fromF(q, Four_Rvects, Four_DMs ):
# res=0
# sr = Four_Rvects.shape
# sm = Four_DMs .shape
# Four_Rvects.shape = (sr[0]*sr[1]*sr[2], 3 )
# Four_DMs .shape = (sm[0]*sm[1]*sm[2], sm[3], sm[4] )
# stime = time.time()
# Ntot = len(Four_DMs)
# icount=0
# for r, dm in zip(Four_Rvects, Four_DMs):
# icount+=1
# if ( 0 and icount and icount%100 ==0):
# etime = time.time()
# print "\rFourier transform DM n# \t ", icount, "/", Ntot, "in ", etime-stime, "seconds. Total is approx.",Ntot*(etime-stime)/icount,
# sys.stdout.flush()
# fact= numpy.dot( r ,q )
# fact=(numpy.exp(+1.0j*fact)).astype( numpy.complex64 )
# res=res+dm*fact
# Four_Rvects.shape = sr
# Four_DMs .shape = sm
# return res
def GetDM_fromF(q, Four_Rvects, Four_DMs , calculatedDatas, direction=None ):
res=0
sr = Four_Rvects.shape
sm = Four_DMs .shape
nyq0,nyq1,nyq2 = Four_Rvects.shape[:3]
nyq0,nyq1,nyq2 =nyq0/2,nyq1/2,nyq2/2
if( numpy.sum(q*q)==0.0):
if direction is not None:
q=q+direction
else:
q=q+numpy.array([0.00001,0,0],"f")
alphas = numpy.array([numpy.dot(r-Four_Rvects[0,0,0],q) for r in [Four_Rvects[1,0,0],Four_Rvects[0,1,0],Four_Rvects[0,0,1]]])
alphas = (numpy.exp(+1.0j*alphas)).astype( numpy.complex64 )
xfacts = numpy.zeros( Four_DMs.shape[0] , numpy.complex64)
yfacts = numpy.zeros( Four_DMs.shape[1] , numpy.complex64)
zfacts = numpy.zeros( Four_DMs.shape[2] , numpy.complex64)
xfacts[:] = alphas[0]
yfacts[:] = alphas[1]
zfacts[:] = alphas[2]
xfacts = numpy.cumprod(xfacts)
yfacts = numpy.cumprod(yfacts)
zfacts = numpy.cumprod(zfacts)
xfacts /= xfacts[nyq0]
yfacts /= yfacts[nyq1]
zfacts /= zfacts[nyq2]
res = numpy.tensordot( Four_DMs , xfacts , axes=[[0],[0]] )
res = numpy.tensordot( res , yfacts , axes=[[0],[0]] )
res = numpy.tensordot( res , zfacts , axes=[[0],[0]] )
corr = calculatedDatas.GetNonAnalytic( q )
# print " res serait "
# print res
# print " corr serait "
# print corr
# print (("%e "*3)%tuple(q.tolist() )) + (("%e "*len(corr[0]))%tuple(corr[0].tolist() ))
res = res + corr
return res
[docs]def CalcDWatT(Temperature , calculatedDatas, filename= None, md5postfix=None, overwrite= False,
key = None, MAKING=0 , simmetries_dict=None):
r"""
The Debye-Waller coefficients are calculated atom by atom.
For each atom they consist in a 3X3 matrix :
.. math:: \left< u_{\alpha}(v) u_{\beta}(v) \right> = \frac{\hbar}{2 M_v } \sum_{ {\bf k},\lambda} \frac{e_{\alpha}(v;{\bf k},\lambda) e^{\star}_{\beta}(v;{\bf k},\lambda) } {\omega_{{\bf k} \lambda}} coth(\frac{\hbar \omega({\bf k},\lambda)}{2 K_{B} T}) weight_{{\bf k}}
where the sum of weigths is one (weigths are given by ab-initio griding).
The Debye-Waller factors are given in atomic units. Internally the routine takes
temperature in Kelvin and converts it to Hartree, while the eigenvetors e are dimensionles,
the frequencies calculated from the ab-initio dynamical matrix are given in units of 1/cm and are converted also to Hartree
"""
# cambia temperatura in Hartree
Temperature_Hartree=Temperature /11604.50520/27.211396132
# amu = 1822.8897312974866 electrons
# cm-1 ==> Hartree cm-1 = 4.5563723521675276e-06 Hartree
# Bohr = 0.529177249 Angstroems
factor_forcoth = 4.5563723521675276e-06*1.0/(2.0*Temperature_Hartree)
factor_for_omegam1 = 1.0/ (2*math.pi*4.5563723521675276e-06)
factor_for_atomicmass =1.0/4 * 1.0/1822.8897312974866
# ----------------------------------------------------
# to be applied on the summed result with masses also
global_factor = factor_for_atomicmass*factor_for_omegam1
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
if( MAKING and (key not in h5 or "DWs" not in h5[key].keys()
or str(Temperature) not in h5[key]["DWs"].keys()
or "DWf_33" not in h5[key]["DWs"][str(Temperature)].keys()
or overwrite) ):
h5.flush()
h5.close()
h5=None
h5=h5py.File(namemd5,'r+')
h5.require_group(key+"/DWs")
h5.require_group(key+"/DWs/"+str(Temperature) )
Nions = len(calculatedDatas.atomAbsolutePositions)
DWf_33 = numpy.zeros([Nions, 3, 3], numpy.float32 )
totQs = h5[key+"/totQs"][:]
totDMs = h5[key+"/totDMs"][:]
stime = time.time()
Ntot = len(totDMs)
icount=0
Nw = len(calculatedDatas.weights)
if ( Ntot%Nw !=0) :
raise Exception, "Ntot%Nw !=0"
sumw=0.0
suma=0
sumb=0
sumc=0
for dm, q in zip(totDMs, totQs):
weig = calculatedDatas.weights[icount % Nw]
sumw+=weig
icount+=1
if (icount and icount%100 ==0):
etime = time.time()
print "\rDeby Waller calculating DM n# \t ", icount, "/", Ntot, "in ", etime-stime, "seconds. Total is approx.",Ntot*(etime-stime)/icount,
sys.stdout.flush()
evals, evects = numpy.linalg.eigh(dm) # gli autovettori qui sono la parte periodica
evals = numpy.sqrt(abs(evals) )
countomega=0
for omega, v in zip(evals, evects.T):
countomega+=1
if( numpy.sum(q*q)> 2*SMALLQ*SMALLQ or countomega>3):
v = numpy.reshape(v, [-1,3])
add=( v.conjugate()[:,:,None] * v[:,None,:]).real
cthfact_over_omega = 1.0/math.tanh( factor_forcoth *omega )/omega *weig
numpy.multiply( add, cthfact_over_omega , add )
numpy.add( DWf_33 , add , DWf_33 )
# print " suma " , suma
# print " sumb " , sumb
# print " sumc " , sumc
# fine loop
numpy.multiply(DWf_33 , 1.0/sumw ,DWf_33 )
# qua la massa dovrebbe essere a denominatore invece
numpy.multiply(DWf_33.transpose() , global_factor / calculatedDatas.atomMasses ,DWf_33.transpose() )
print " SUMW " , sumw
try:
del h5[key]["DWs"][str(Temperature)]["DWf_33"]
print " YOU NEED TO REPACK LATER "
except:
pass
h5[key]["DWs"][str(Temperature)]["DWf_33"]=DWf_33
res={}
res["DWf_33"]=h5[key]["DWs"][str(Temperature)]["DWf_33"]
h5.flush()
h5.close()
return res
[docs]def GetAllReplicated ( calculatedDatas, simmetries_dict , filename= None, Nfour_interp=None,
md5postfix=None,
overwrite= False,
key = None ,
APPLYTIMEREVERSAL=1,
CALCULATEFOURIER=1,
CALCULATECOMPLEMENT=0,
MAKING=0,
MAKINGFOURIER=0,
ZBORN_info=None
):
""" This routine writes informations about the whole Brillouin zone
that has been obtained by replication. The output is written into a hdf5 file.
Some of the most important calculated data, for further calculations are :
* h5[key+"/totQs"] = totQs
* h5[key+"/totDMs"] = totDMs
where totQs is a list of Q vectors, totDMs is a list of dynamical matrices ( units??? )
The output file name is filename+"."+md5postfix, where filename is the ab-initio (CASTEP, phonopy) filename and md5postfix is its md5 sum
"""
Qs = calculatedDatas.GetClippedQs()
origQs = calculatedDatas.GetNonClippedQs()
shiftsToBz = Qs-origQs
# exp( -jorigQs.R) * exp(-j*shifttobz.R)
# la linea sotto perche' ab-initio da gli autovettori non periodici
### THIS FACTOR SHOULD BE REARRANGED BECAUSE
# THE DM IS REMULTIPLIED BY ANIOTHER PHASE BEFORE USE
# AND SUCH PHASE SHOUDL BE REGROUPED HERE
shiftsToBz = shiftsToBz+origQs
print " Qs.shape ", Qs.shape
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
if( MAKING and (key not in h5 or overwrite) ):
if nprocs>1:
raise "mpi not yet activate in Simmetrization "
reducedQs = calculatedDatas.Reduce(Qs)
###############################33
irot=0
symm_list =simmetries_dict["DiffRotList"] # [:1]
if APPLYTIMEREVERSAL:
segni=[-1,1]
else:
segni=[1]
Nsymm = len(symm_list)
totQs = numpy.zeros( [ Nsymm * len(Qs ) *len(segni) ,3], numpy.float32 )
totWeights = numpy.zeros( [ Nsymm * len(Qs ) *len(segni) ], numpy.float32 )
totIrot = numpy.ones( [ Nsymm * len(Qs ) *len(segni) ] , numpy.int32 )
totItime = numpy.ones( [ Nsymm * len(Qs ) *len(segni) ] , numpy.int32 )
totBVifacts = numpy.zeros( [ Nsymm * len(Qs ) *len(segni) ,3], numpy.int32 )
print "nofbranches ", calculatedDatas.NofBranches
totDMs = numpy.zeros( [ Nsymm * len(Qs ) *len(segni) , calculatedDatas.NofBranches,
calculatedDatas.NofBranches ], numpy.complex64)
isymm=0
BV=calculatedDatas.GetBrillVects()
## for BVifacts in [(0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,1,0),(1,0,1),(0,1,1),(1,1,1)][:1] :
orig_evects = calculatedDatas.eigenvectors
evals = calculatedDatas.frequencies
for BVifacts in [ ( 0 , 0 , 0 ) ] :
SHIFT = numpy.dot( BVifacts, BV )
irot=-1
for Rot in symm_list:
irot=irot+1
print "Replicating with rotation ", irot
wgw = simmetries_dict["wgw"][irot]
new_dm = compose_dm ( evals, orig_evects , wgw,
calculatedDatas.atomAbsolutePositions,
shiftsToBz )
for segno in segni :
repQs = numpy.dot( Rot, Qs.T ).T
totWeights [ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = calculatedDatas.weights
totQs [ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = repQs*segno
totIrot [ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = isymm+1
totItime[ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = segno
totBVifacts[ isymm * len(Qs ) : ( isymm+1 )*len(Qs ) ] = BVifacts
if segno==-1:
new_dm_segno = new_dm.conjugate()
else:
new_dm_segno = new_dm
totDMs[ isymm * len(Qs ) : ( isymm+1 )*len(Qs ), :,:]= new_dm_segno
isymm=isymm+1
offsets, steps = getOffsetSteps(reducedQs)
step = numpy.dot(steps , BV )
step = numpy.sqrt( numpy.sum( step*step, axis = -1 ) )
AtotQs=[]
AtotIrot = []
AtotItime =[]
AtotBVifacts=[]
AtotDMs=[]
dist1 = numpy.sqrt( numpy.sum( Qs * Qs, axis = -1 ) )
if CALCULATECOMPLEMENT:
for ix in [-1,0,1]:
for iy in [-1,0,1]:
for iz in [-1,0,1]:
if (ix,iy,iz)==(0,0,0):
continue
BVifacts = (ix, iy, iz)
qshift = numpy.dot( BVifacts, BV )
neiQs = Qs + qshift
dist2 = numpy.sqrt( numpy.sum( neiQs*neiQs, axis = -1 ) )
diff = numpy.abs(dist2-dist1)
mask = (diff<3*step)
piece = neiQs[ mask ]
# qshifted_evects = self.shift_evects(orig_evects, qshift) # si usera self.PositionsList_flat
if(len(piece)):
isymm=0
irot=-1
for Rot in symm_list:
irot=irot+1
wgw = simmetries_dict["wgw"][irot]
print "replicating complement with rotation ", irot
repQs = numpy.dot( Rot, piece.T ).T
new_dm = compose_dm ( evals[mask], orig_evects[mask] , wgw,
calculatedDatas.atomAbsolutePositions,
qshift+shiftsToBz [mask] )
for segno in segni :
AtotQs.extend(repQs*segno )
AtotIrot.extend( [isymm+1]* len( repQs) )
AtotItime.extend( [segno]*len(repQs ) )
AtotBVifacts.extend( [ BVifacts] *len(repQs ) )
if segno==-1:
new_dm_segno = new_dm.conjugate()
else:
new_dm_segno = new_dm
AtotDMs.extend( numpy.array(new_dm_segno))
isymm=isymm+1
h5.flush()
h5.close()
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
h5.require_group(key)
for chiave in h5[key].keys():
if chiave != "complemento":
del h5[key][chiave]
h5[key+"/totWeights"] = totWeights
h5[key+"/totQs"] = totQs
h5[key+"/totIrot"] = totIrot
h5[key+"/totItime"]= totItime
h5[key+"/totBVifacts"]= totBVifacts
h5[key+"/totDMs"]= totDMs
h5[key+"/steps"]= steps
if CALCULATECOMPLEMENT:
h5.require_group(key+"/complemento")
for chiave in h5[key+"/complemento"].keys():
del h5[key+"/complemento"][chiave]
h5[key+"/complemento/totQs"] = numpy.array(AtotQs )
h5[key+"/complemento/totIrot"] = numpy.array(AtotIrot )
h5[key+"/complemento/totItime"] = numpy.array(AtotItime )
h5[key+"/complemento/totBVifacts"]= numpy.array(AtotBVifacts)
h5[key+"/complemento/totDMs"]= numpy.array(AtotDMs)
res={}
res["totQs" ]= totQs = h5[key+"/totWeights"][:]
res["totQs" ]= totQs = h5[key+"/totQs"][:]
res["totDMs"]= totDMs =h5[key+"/totDMs"][:]
res["totIrot" ]= h5[key+"/totIrot"][:]
res["totItime" ]= h5[key+"/totItime"][:]
res["totBVifacts"]= h5[key+"/totBVifacts"][:]
if CALCULATEFOURIER:
if( MAKINGFOURIER and ( overwrite or "Fourier_Transformations" not in h5[key].keys() or (str(Nfour_interp) not in h5[key+"/Fourier_Transformations"].keys())) ):
if(myrank==0):
h5.flush()
h5.close()
h5=h5py.File(namemd5,'r+')
h5.require_group(key+"/Fourier_Transformations")
Four_Rvects = numpy.zeros( [ 2*Nfour_interp +1, 2*Nfour_interp +1, 2*Nfour_interp +1, 3 ] , numpy.float32)
cellvects = calculatedDatas.cellVects
for j in range(-Nfour_interp , Nfour_interp +1):
for k in range(-Nfour_interp , Nfour_interp +1):
for l in range(-Nfour_interp , Nfour_interp +1):
Four_Rvects [ j+Nfour_interp, k+Nfour_interp, l+Nfour_interp] = j*cellvects[0]+k*cellvects[1]+l*cellvects[2]
Four_DMs = numpy.zeros( [ 2*Nfour_interp+1, 2*Nfour_interp+1, 2*Nfour_interp+1,
calculatedDatas.NofBranches , calculatedDatas.NofBranches ] , numpy.complex64)
stime = time.time()
Ntot = len(totDMs)
icount=0
Nw = len(calculatedDatas.weights)
if ( Ntot%Nw !=0) :
raise Exception, "Ntot%Nw !=0"
sumw=0.0
chunck = max(len(totQs)/100,1)
inizio=-chunck
while(inizio+chunck<len(totQs)):
inizio= inizio+chunck
fine = inizio+chunck
if fine>len(totQs):
fine = len(totQs)
nonAs = calculatedDatas.GetNonAnalytic(totQs[inizio:fine])
for dm, q, nonA in zip(totDMs[inizio:fine], totQs[inizio:fine], nonAs):
# for dm, q in zip(totDMs, totQs):
weig = calculatedDatas.weights[icount % Nw]
sumw+=weig
icount+=1
if (icount and icount%100 ==0):
etime = time.time()
print "\rFourier transform DM n# \t ", icount, "/", Ntot, "in ", etime-stime, "seconds. Total is approx.",Ntot*(etime-stime)/icount,
sys.stdout.flush()
facts = numpy.tensordot( q , calculatedDatas.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 , )
dm = dm * facts.conjugate()
dm = numpy.swapaxes(dm,-1,-2) * facts
dm = numpy.swapaxes(dm,-1,-2)
dm = dm - nonA
rcount=0
for j in range(-Nfour_interp , Nfour_interp +1):
for k in range(-Nfour_interp , Nfour_interp +1):
for l in range(-Nfour_interp , Nfour_interp +1):
rcount=rcount+1
if rcount%nprocs!=myrank:
continue
four_rv = Four_Rvects [ j+Nfour_interp, k+Nfour_interp, l+Nfour_interp]
fact= numpy.dot( four_rv ,q )
fact=(numpy.exp(-1.0j*fact)).astype( numpy.complex64 )
Four_DMs [j+Nfour_interp, k+Nfour_interp, l+Nfour_interp]+= dm*(fact*weig)
if nprocs>1:
target = numpy.zeros( Four_DMs.shape , "F")
comm.Reduce( [ Four_DMs , MPI.COMPLEX8 ] , [ target , MPI.COMPLEX8 ], op=MPI.SUM, root=0)
Four_DMs[:]=target
# target = numpy.zeros( Four_DMs.shape , "f")
# comm.Reduce( [ Four_DMs.real , MPI.FLOAT ] , [ target , MPI.FLOAT ], op=MPI.SUM, root=0)
# Four_DMs.real = target
# comm.Reduce( [ Four_DMs.imag , MPI.FLOAT ] , [ target , MPI.FLOAT ], op=MPI.SUM, root=0)
# Four_DMs.imag = target
print ""
print " SUMW " , sumw
numpy.multiply( Four_DMs , numpy.array([1.0/ sumw], numpy.float32) , Four_DMs )
h5.require_group(key+"/Fourier_Transformations/"+str(Nfour_interp))
for chiave in h5[key+"/Fourier_Transformations/"+str(Nfour_interp)].keys():
del h5[key+"/Fourier_Transformations/"+str(Nfour_interp)][chiave]
if(myrank==0) :
h5[key+"/Fourier_Transformations/"+str(Nfour_interp)+"/Four_DMs"]= Four_DMs
h5[key+"/Fourier_Transformations/"+str(Nfour_interp)+"/Four_Rvects"]=Four_Rvects
if(nprocs>1):
h5.close()
comm.Barrier()
return res
res["Four_DMs"] =h5[key+"/Fourier_Transformations/"+str(Nfour_interp)+"/Four_DMs"][:]
res["Four_Rvects"]=h5[key+"/Fourier_Transformations/"+str(Nfour_interp)+"/Four_Rvects"][:]
res["steps"]= h5[key+"/steps"][:]
if CALCULATECOMPLEMENT:
res["complemento"]={}
print "overwrite ", overwrite
print "leggo ", key+"/complemento/totQs"
res["complemento"]["totQs" ]= h5[key+"/complemento/totQs"][:]
res["complemento"]["totIrot" ]= h5[key+"/complemento/totIrot"][:]
res["complemento"]["totItime" ]= h5[key+"/complemento/totItime"][:]
res["complemento"]["totBVifacts"]= h5[key+"/complemento/totBVifacts"][:]
res["complemento"]["totDMs"]= h5[key+"/complemento/totDMs"][:]
h5.flush()
h5.close()
h5=None
return res
def GetSimplyDMs ( calculatedDatas, filename= None,
md5postfix=None,
overwrite= False, acustic = False
):
"""
"""
Qs = calculatedDatas.GetClippedQs()
origQs = calculatedDatas.GetNonClippedQs()
shiftsToBz = Qs-origQs
# exp( -jorigQs.R) * exp(-j*shifttobz.R)
# la linea sotto perche' ab-initio da gli autovettori non periodici
### THIS FACTOR SHOULD BE REARRANGED BECAUSE
# THE DM IS REMULTIPLIED BY ANIOTHER PHASE BEFORE USE
# AND SUCH PHASE SHOUDL BE REGROUPED HERE
shiftsToBz = shiftsToBz+origQs
print " Qs.shape ", Qs.shape
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
if( ("simplydms" not in h5 or overwrite) ):
if nprocs>1:
raise "mpi not yet activate in Simmetrization "
h5.flush()
h5.close()
reducedQs = calculatedDatas.Reduce(Qs)
###############################33
nofbranches = calculatedDatas.NofBranches
print "nofbranches ", calculatedDatas.NofBranches
totDMs = numpy.zeros( [ len(Qs ) , calculatedDatas.NofBranches,
calculatedDatas.NofBranches ], numpy.complex64)
BV=calculatedDatas.GetBrillVects()
orig_evects = calculatedDatas.eigenvectors
evals = calculatedDatas.frequencies
wgw = numpy.eye( nofbranches, dtype = "d")
dms = compose_dm ( evals, orig_evects , wgw,
calculatedDatas.atomAbsolutePositions,
shiftsToBz , acustic = acustic )
namemd5=filename+"."+md5postfix
h5=h5py.File(namemd5,'r+')
h5.require_group("simplydms")
for chiave in h5["simplydms"].keys():
del h5["simplydms"][chiave]
h5[ "simplydms" + "/dms" ] = dms
h5[ "simplydms" + "/Qs" ] = Qs
h5.close()
h5=h5py.File(namemd5,'r+')
res={}
res["dms" ] = h5[ "simplydms" + "/dms" ][:]
res["Qs" ] = h5[ "simplydms" + "/Qs" ][:]
dms = h5[ "simplydms" + "/dms" ][:]
h5.flush()
h5.close()
orig_evects = calculatedDatas.eigenvectors
evals = calculatedDatas.frequencies
print " SHAPES "
print Qs.shape
print dms.shape
print evals.shape
print orig_evects.shape
res["evals" ] = evals
return res
def getOffsetSteps(Qs):
offsets=Qs[0]
steps=numpy.zeros(3,numpy.float32)
oldq=offsets
for q in Qs[1:]:
for i in range(3):
step= abs( q[i]-oldq[i] )
if step<1.0e-4:
continue
if(steps[i]==0) :
steps[i]=step
elif ( step < steps[i] ):
steps[i]=step
oldq=q
offsets=offsets%steps
return offsets, steps