Changeset 1075 for trunk/GSASIImath.py


Ignore:
Timestamp:
Oct 2, 2013 3:42:53 PM (8 years ago)
Author:
vondreele
Message:

implement new March-Dollase term in MC/SA

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1070 r1075  
    26392639                        parmDict[pfx+name] = A
    26402640       
    2641     def mcsaCalc(values,refList,rcov,ifInv,allFF,RBdata,varyList,parmDict):
     2641    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
    26422642        ''' Compute structure factors for all h,k,l for phase
    26432643        input:
     
    26592659        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
    26602660        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
    2661         MDaxis = parmDict['0:MDaxis']
    2662         Gmat = parmDict['Gmat']
    26632661        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
    26642662        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
     
    26682666            Bterm = refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
    26692667            refList[5] = Aterm+Bterm
    2670         #apply MD correction here
     2668        if len(cosTable):        #apply MD correction
     2669            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
    26712670        sumFcsq = np.sum(refList[5])
    26722671        scale = parmDict['sumFosq']/sumFcsq
     
    26822681    generalData = data['General']
    26832682    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
    2684     Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
     2683    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
    26852684    SGData = generalData['SGData']
    2686     SGM = [SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))]
    2687     SGT = [SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))]
     2685    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
     2686    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
     2687    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
    26882688    fixAtoms = data['Atoms']                       #if any
    26892689    cx,ct,cs = generalData['AtomPtrs'][:3]
     
    27092709    upper = []
    27102710    lower = []
     2711    MDvec = np.zeros(3)
    27112712    for i,item in enumerate(MCSAObjs):
    27122713        mfx = str(i)+':'
     
    27142715        if item['Type'] == 'MD':
    27152716            getMDparms(item,mfx,parmDict,varyList)
     2717            MDvec = np.array(item['axis'])
    27162718        elif item['Type'] == 'Atom':
    27172719            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
     
    27302732    refs = []
    27312733    allFF = []
     2734    cosTable = []
    27322735    sumFosq = 0
    27332736    if 'PWDR' in reflName:
     
    27402743                refs.append([h,k,l,m,f*m,pos,sig])
    27412744                sumFosq += f*m
     2745                Heqv = np.inner(np.array([h,k,l]),SGMT)
     2746                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
    27422747        nRef = len(refs)
     2748        cosTable = np.array(cosTable)**2
    27432749        rcov = np.zeros((nRef,nRef))
    27442750        for iref,refI in enumerate(refs):
     
    27662772                refs.append([h,k,l,m,f*m,iref,0.])
    27672773                sumFosq += f*m
     2774                Heqv = np.inner(np.array([h,k,l]),SGMT)
     2775                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
     2776        cosTable = np.array(cosTable)**2
    27682777        nRef = len(refs)
    27692778        if covData['freshCOV'] and generalData['doPawley'] and MCSA.get('newDmin',True):
     
    27982807    x0 = [parmDict[val] for val in varyList]
    27992808    ifInv = SGData['SGInv']
    2800     results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,allFF,RBdata,varyList,parmDict),
     2809    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
    28012810        schedule=MCSA['Algorithm'], full_output=True,
    28022811        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
     
    28052814        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
    28062815        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
    2807     M = mcsaCalc(results[0],refs,rcov,ifInv,allFF,RBdata,varyList,parmDict)
     2816    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
    28082817#    for ref in refs.T:
    28092818#        print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6])
Note: See TracChangeset for help on using the changeset viewer.