Changeset 391 for trunk/GSASIIstruct.py


Ignore:
Timestamp:
Oct 13, 2011 10:48:51 AM (10 years ago)
Author:
vondreele
Message:

set default atom frac=1 for cif files
store covMatrix not normalized by diagonal
make esds on cell parameters

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstruct.py

    r390 r391  
    473473    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables
    474474   
    475 def SetPhaseData(parmDict,sigDict,Phases):
     475def getVCov(varyNames,varyList,covMatrix):
     476    vcov = np.zeros((len(varyNames),len(varyNames)))
     477    for i1,name1 in enumerate(varyNames):
     478        for i2,name2 in enumerate(varyNames):
     479            try:
     480                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
     481            except ValueError:
     482                vcov[i1][i2] = 0.0
     483    return vcov
     484   
     485def getCellEsd(pfx,SGData,A,covData):
     486    dpr = 180./np.pi
     487    rVsq = G2lat.calc_rVsq(A)
     488    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
     489    cell = np.array(G2lat.Gmat2cell(g))   #real cell
     490    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
     491    scos = cosd(cellst[3:6])
     492    ssin = sind(cellst[3:6])
     493    scot = scos/ssin
     494    rcos = cosd(cell[3:6])
     495    rsin = sind(cell[3:6])
     496    rcot = rcos/rsin
     497    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
     498    varyList = covData['varyList']
     499    covMatrix = covData['covMatrix']
     500    vcov = getVCov(RMnames,varyList,covMatrix)
     501    Ax = np.array(A)
     502    Ax[3:] /= 2.
     503    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
     504        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
     505    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
     506    Vol = 1/np.sqrt(rVsq)
     507    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
     508    R123 = Ax[0]*Ax[1]*Ax[2]
     509    dsasdg = np.zeros((3,6))
     510    dadg = np.zeros((6,6))
     511    for i0 in range(3):         #0  1   2
     512        i1 = (i0+1)%3           #1  2   0
     513        i2 = (i1+1)%3           #2  0   1
     514        i3 = 5-i2               #3  5   4
     515        i4 = 5-i1               #4  3   5
     516        i5 = 5-i0               #5  4   3
     517        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
     518        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
     519        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
     520        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
     521        denom = np.sqrt(denmsq)
     522        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
     523        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
     524        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
     525        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
     526        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
     527        dadg[i5][i5] = -Ax[i0]/denom
     528    for i0 in range(3):
     529        i1 = (i0+1)%3
     530        i2 = (i1+1)%3
     531        i3 = 5-i2
     532        for ij in range(6):
     533            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
     534            if ij == i0:
     535                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
     536            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
     537    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
     538    var = np.diag(sigMat)
     539    CS = np.where(var>0.,np.sqrt(var),0.)
     540    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
     541    return cellSig           
     542   
     543def SetPhaseData(parmDict,sigDict,Phases,covData):
    476544   
    477545    def cellFill(pfx,SGData,parmDict,sigDict):
     
    572640        if cell[0]:
    573641            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
     642            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
    574643            print ' Reciprocal metric tensor: '
    575644            ptfmt = "%15.9f"
     
    590659            cell[1:7] = G2lat.A2cell(A)
    591660            cell[7] = G2lat.calc_V(A)
    592             print ' New unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]), \
    593                 ' c = %.5f'%(cell[3]),' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]), \
    594                 ' gamma = %.3f'%(cell[6]),' volume = %.3f'%(cell[7])
     661            print ' New unit cell:'
     662            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
     663            names = ['a','b','c','alpha','beta','gamma','Volume']
     664            namstr = '  names :'
     665            ptstr =  '  values:'
     666            sigstr = '  esds  :'
     667            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
     668                namstr += '%12s'%(name)
     669                ptstr += fmt%(a)
     670                if siga:
     671                    sigstr += fmt%(siga)
     672                else:
     673                    sigstr += 12*' '
     674            print namstr
     675            print ptstr
     676            print sigstr
     677           
    595678        if 'Pawley' in Phase['General']['Type']:
    596679            pawleyRef = Phase['Pawley ref']
     
    17621845            ymb = np.array(y-yb)
    17631846            ycmb = np.array(yc-yb)
    1764             ratio = np.where(ycmb>0.,ymb/ycmb,1.0)           
     1847            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
    17651848            xB = np.searchsorted(x,Limits[0])
    17661849            xF = np.searchsorted(x,Limits[1])
     
    17891872                            iFin2 = np.searchsorted(x,pos2+fmax)
    17901873                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
    1791                         refl[8] = np.sum(yp[iBeg:iFin2]*ratio[iBeg:iFin2])/(refl[13]*(1.+kRatio))
     1874                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
    17921875                    elif 'T' in calcControls[hfx+'histType']:
    17931876                        print 'TOF Undefined at present'
     
    20832166            covMatrix = result[1]*GOF
    20842167            sig = np.sqrt(np.diag(covMatrix))
    2085             xvar = np.outer(sig,np.ones_like(sig))
    2086             cov = np.divide(np.divide(covMatrix,xvar),xvar.T)
    20872168            if np.any(np.isnan(sig)):
    20882169                print '*** Least squares aborted - some invalid esds possible ***'
     
    21062187    GetFobsSq(Histograms,Phases,parmDict,calcControls)
    21072188    sigDict = dict(zip(varyList,sig))
    2108     SetPhaseData(parmDict,sigDict,Phases)
     2189    covData = {'variables':result[0],'varyList':varyList,'covMatrix':covMatrix}
     2190    SetPhaseData(parmDict,sigDict,Phases,covData)
    21092191    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
    21102192    SetHistogramData(parmDict,sigDict,Histograms)
    2111     covData = {'variables':result[0],'varyList':varyList,'covariance':cov}
    21122193    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
    21132194#for testing purposes!!!
Note: See TracChangeset for help on using the changeset viewer.