Skip to content
Snippets Groups Projects
Commit f2b13c88 authored by alex mirone's avatar alex mirone
Browse files

finito calcolo secondo ordine ma mettere multiresoluzione

parent 570a77fa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ab2tds/make_TDS_tds2elHarvest_Clone.py
+ 179
− 23
View file @ f2b13c88
......@@ -109,17 +109,19 @@ except:
class LUT:
def __init__(self, qs, freqs, vects ):
def __init__(self, qs_reds, qs, freqs, vects ):
self.qs_reds = qs_reds
self.qs = qs
self.freqs = freqs
self.vects = vects
self.origin = qs.min( axis=0 )
self.ends = qs.max( axis=0 )
self.origin = qs_reds.min( axis=0 )
self.ends = qs_reds.max( axis=0 )
self.steps = []
for i in range(3):
xs = numpy.array(qs[:,i]).sort()
xs = numpy.array(qs_reds[:,i]).sort()
dists = xs[1:]-xs[:-1]
dists = dists[dists>1.0e-3]
self.steps.append( dists.min() )
self.steps = numpy.array( self.steps )
......@@ -128,7 +130,7 @@ class LUT:
self.lut = numpy.zeros( shape )
self.lut[:]=-1
for i, q in enumerate(qs):
for i, q in enumerate(qs_reds):
w = ( q-self.origin )/self.steps
iw = numpy.round(w).astype("i")
if numpy.abs(w-iw).sum()>0.1:
......@@ -138,24 +140,52 @@ class LUT:
def getFV(self,q):
w = ( q-self.origin )/self.steps
if numpy.less( w,0 ).sum():
if numpy.less( w, -1.0e-3 ).sum():
return None
if numpy.less( self.ends + 1.0e-3 ,q ).sum():
return None
iw = numpy.round(w).astype("i")
if numpy.less( self.ends + self.steps ,w ).sum():
return None
if numpy.abs(w-iw).sum()>0.1:
if numpy.abs(w-iw).sum() > 1.0e-3 :
raise Exception, "LUT is not working properly"
i = self.lut[tuple(iw)]
if i==-1:
return None
return self.qs[i], self.freqs[i], self.vects[i]
return self.qs_reds[i], self.qs[i], self.freqs[i], self.vects[i]
def lut_getFV(q,
lut_origin,
lut_steps,
lut_ends,
lut_lut,
lut_qs_reds, lut_qs, lut_freqs, lut_vects):
w = ( q-lut_origin )/lut_steps
if numpy.less( w, -1.0e-3 ).sum():
return None
if numpy.less( lut_ends + 1.0e-3 ,q ).sum():
return None
iw = numpy.round(w).astype("i")
if numpy.abs(w-iw).sum() > 1.0e-3 :
raise Exception, "LUT is not working properly"
i = lut_lut[tuple(iw)]
if i==-1:
return None
return lut_qs_reds[i], lut_qs[i], lut_freqs[i], lut_vects[i]
if(0):
Temperature=100
......@@ -201,8 +231,6 @@ fitSpots = [[-2,-1,-2],[-1,0,4],[0,1,-8] ,[2,0,-2], [-2,-1,-2], ]
"""
"""
Ctransf = None
"""
if Ctransf is None, no transformation is done. This is equivalent to do the transformation [ [1,0,0],[0,1,0],[0,0,1] ].
......@@ -218,6 +246,106 @@ Ctransf = None
"""
def secondorder(qnotshifted,red, qshifted, Qs, Qs_reds,
lut , massfactors, dwfacts_mass_kin_scatt,factor_forcoth,
lut_origin,
lut_steps,
lut_ends,
lut_lut,
lut_qs_reds, lut_qs, lut_freqs, lut_vects ):
add=0.0
for q1_red in Qs_reds:
q2_red = res-q1_red
if 1:
q1_data = lut.getFV(q1_red)
q2_data = lut.getFV(q2_red)
else:
q1_data = lut_getFV(q1_red,
lut_origin,
lut_steps,
lut_ends,
lut_lut,
lut_qs_reds, lut_qs, lut_freqs, lut_vects)
q2_data = lut_getFV(q2_red,
lut_origin,
lut_steps,
lut_ends,
lut_lut,
lut_qs_reds, lut_qs, lut_freqs, lut_vects)
if q1_data is None or q2_data is None:
continue
q1_red_p , q1, evals1, evects1 = q1_data
q2_red_p , q2, evals2, evects2 = q2_data
assert( (abs(q1_red_p-q1_red ) ).sum<1.0e-6 )
assert( (abs(q2_red_p-q2_red ) ).sum<1.0e-6 )
# originally the procedure beklow which has been copied,
# was getting vects obtained from numpy.linalg.eig
# which gives vectors as columns
evects1=numpy.transpose(evects1)
evects2=numpy.transpose(evects2)
evals1=(numpy.maximum(evals1, numpy.array([bottom_meV/(0.0001239852 *1000)], numpy.float32) ))
evals2=(numpy.maximum(evals2, numpy.array([bottom_meV/(0.0001239852 *1000)], numpy.float32) ))
qvect_s1 = numpy.tensordot( qshifted1 , evects1.reshape([-1,3, Nbranches]) , [[0], [1] ] )
qvect_s2 = numpy.tensordot( qshifted2 , evects1.reshape([-1,3, Nbranches]) , [[0], [1] ] )
Fjs = numpy.sum(dwfacts_mass_kin_scatt* (qvect_s1.T * qvect_s2.T ) * massfactors, axis=-1) # sommato sugli ioni e con un mass fctors aggiuntivo
##### MASSA ATOMICA 1,007 825 032 u (1,673 534·10-27 kg)
##### h plank 6.626070040(81)×10−34 J s
##### 4.135667662(25)×10−15 eV s
##### c = 29979245800 cm/s
hplanck = 6.626070041e-27 ## ergs seconds
clight = 2.9979245800e+10 ## cm/seconds
massaatomica = 1.673534e-24 /1.007825032 ## grammi
## deve equilibrare gli (2pi A-1)^2 di E.Q ei fattori 2( che va messo), e l'omega in cm-1 (senza pi) al denominatore
equilibrating_factor = ( hplanck/( 2* clight * massaatomica *(2*math.pi)*(2*math.pi) ) ) *1.0e16
#########
## il 2pi al quadrato perche c'e' h invece di hbar e gli omega sono dati in cm-1 senza 2pi.
## la velocita della luce perche' gli omega sono in cm-1. La massa atomica perche i massfactors
## sono formati con quella relativa solamente, infine il 2 a denominatore perche' nella formula c'e'
## Si ottengono cosi' dei cm**2 e si motiplica per 1.0e16 per trasformarli in A-1
Fjs = Fjs * equilibrating_factor
## va applicato prima del quadrato perche ci sono i due fattori sqrt(hbar/omega/m/2 ) E.Q
intensity_array = ( Fjs*(Fjs.conjugate())).real
exp_plus1 = numpy.exp( factor_forcoth * evals1)
exp_minus1 = numpy.exp(-factor_forcoth * evals1)
deno1 = (exp_plus1-exp_minus1)*evals1
totalstats1 = ( exp_plus1+exp_minus1)/deno1
exp_plus2 = numpy.exp( factor_forcoth * evals2)
exp_minus2 = numpy.exp(-factor_forcoth * evals2)
deno2 = (exp_plus2-exp_minus2)*evals2
totalstats2 = ( exp_plus2+exp_minus2)/deno2
totalstats = totalstats1*totalstats2
add += numpy.sum( intensity_array* totalstats ) * equilibrating_factor
return add
if(sys.argv[0][-12:]!="sphinx-build"):
s=open(sys.argv[2], "r").read()
......@@ -313,6 +441,7 @@ if(sys.argv[0][-12:]!="sphinx-build"):
if iprog%100 ==0:
print "iprog ", iprog
vals, vects = numpy.linalg.eigh(dm)
vals = numpy.sqrt( numpy.maximum( vals,0 ) )
frequencies.append(vals)
eigenvectors.append(vects.T)
......@@ -323,10 +452,10 @@ if(sys.argv[0][-12:]!="sphinx-build"):
print " maxBrill " , maxBrill
print " norme.max(), norme.min() ", norme.max(), norme.min()
reds = numpy.abs(numpy.tensordot( Qs, BVinv_simple, axes = [[-1],[-1]] ))
reds = numpy.max(reds, axis=-1)
Qs_reds = numpy.abs(numpy.tensordot( Qs, BVinv_simple, axes = [[-1],[-1]] ))
reds = numpy.max(Qs_reds, axis=-1)
lut = LUT(calculatedDatas.QsReduce,frequencies ,eigenvectors)
lut = LUT( Qs_reds, Qs, frequencies ,eigenvectors)
Qs_roi = Qs[ numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max ) *numpy.less(reds, 0.5) ]
frequencies_roi = frequencies[numpy.less(qradius_min, norme)*numpy.less(norme, qradius_max )*numpy.less(reds, 0.5) ]
......@@ -382,7 +511,6 @@ if(sys.argv[0][-12:]!="sphinx-build"):
print " MASSES "
print calculatedDatas.atomMasses
h5=h5py.File(sys.argv[3] , "r")
DWs_3X3 = h5["Codes_"+(str(APPLYTIMEREVERSAL))]["DWs"][str(Temperature)]["DWf_33"][:]
h5.close()
......@@ -402,7 +530,6 @@ if(sys.argv[0][-12:]!="sphinx-build"):
else:
dwfacts_mass_kin_scatt_perq = dwfacts_perq*massfactors*kinematicalFactors*scattFactors_q_site
# fake_int = scattFactors_q_site.sum(axis=-1)
print dwfacts_perq[0]
......@@ -417,6 +544,7 @@ if(sys.argv[0][-12:]!="sphinx-build"):
# 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)
resevals=[]
......@@ -456,12 +584,40 @@ if(sys.argv[0][-12:]!="sphinx-build"):
kinantistk= 1.0
deno = (exp_plus-exp_minus)*evals
totalstats = ( exp_plus*kinstokes+exp_minus*kinantistk)/deno
intensity_line[count] = numpy.sum( intensity_array* totalstats )
if DOSECONDORDER:
##### MASSA ATOMICA 1,007 825 032 u (1,673 534·10-27 kg)
##### h plank 6.626070040(81)×10−34 J s
##### 4.135667662(25)×10−15 eV s
##### c = 29979245800 cm/s
hplanck = 6.626070041e-27 ## ergs seconds
clight = 2.9979245800e+10 ## cm/seconds
massaatomica = 1.673534e-24 /1.007825032 ## grammi
## deve equilibrare gli (2pi A-1)^2 di E.Q ei fattori 2, e l'omega in cm-1 (senza pi) al denominatore
equilibrating_factor = ( hplanck/( 2* clight * massaatomica *(2*math.pi)*(2*math.pi) ) ) *1.0e16
#########
## il 2pi al quadrato perche c'e' h invece di hbar e gli omega sono dati in cm-1 senza 2pi.
## la velocita della luce perche' gli omega sono in cm-1. La massa atomica perche i massfactors
## sono formati con quella relativa solamente, infine il 2 a denominatore perche' nella formula c'e'
## Si ottengono cosi' dei cm**2 e si motiplica per 1.0e16 per trasformarli in A-1
intensity_line[count] = numpy.sum( intensity_array * totalstats * equilibrating_factor )
if DOSECONDORDER:
red = numpy.dot( BVinv , qnotshifted )
add = secondorder(qnotshifted,red, qshifted, Qs, calculatedDatas.QsReduced,
lut , massfactors, dwfacts_mass_kin_scatt,factor_forcoth )
add = secondorder(qnotshifted,red, qshifted, Qs, Qs_reds,
lut , massfactors, dwfacts_mass_kin_scatt,factor_forcoth,
lut.origin,
lut.steps,
lut.ends,
lut.lut,
lut.qs_reds, lut.qs, lut.freqs, lut.vects
)
intensity_line[count] += add
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment