Source code for TDS_Simmetry


#/*##########################################################################
# 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