Changeset 1062


Ignore:
Timestamp:
Sep 20, 2013 3:28:58 PM (8 years ago)
Author:
vondreele
Message:

a better MC/SA? Still checking it.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1061 r1062  
    25302530        return np.array(Mdata)
    25312531       
    2532     def GetAtomTX(RBdata,parmDict):
     2532    def GetAtomT(RBdata,parmDict):
     2533        'Needs a doc string'
     2534        atNo = parmDict['atNo']
     2535        nfixAt = parmDict['nfixAt']
     2536        Tdata = atNo*[' ',]
     2537        for iatm in range(nfixAt):
     2538            parm = ':'+str(iatm)+':Atype'
     2539            if parm in parmDict:
     2540                Tdata[iatm] = aTypes.index(parmDict[parm])
     2541        iatm = nfixAt
     2542        for iObj in range(parmDict['nObj']):
     2543            pfx = str(iObj)+':'
     2544            if parmDict[pfx+'Type'] in ['Vector','Residue']:
     2545                if parmDict[pfx+'Type'] == 'Vector':
     2546                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
     2547                    nAtm = len(RBRes['rbVect'][0])
     2548                else:       #Residue
     2549                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
     2550                    nAtm = len(RBRes['rbXYZ'])
     2551                for i in range(nAtm):
     2552                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
     2553                    iatm += 1
     2554            elif parmDict[pfx+'Type'] == 'Atom':
     2555                atNo = parmDict[pfx+'atNo']
     2556                parm = pfx+'Atype'              #remove extra ':'
     2557                if parm in parmDict:
     2558                    Tdata[atNo] = aTypes.index(parmDict[parm])
     2559                iatm += 1
     2560            else:
     2561                continue        #skips March Dollase
     2562        return Tdata
     2563       
     2564    def GetAtomX(RBdata,parmDict):
    25332565        'Needs a doc string'
    25342566        Bmat = parmDict['Bmat']
    25352567        atNo = parmDict['atNo']
    25362568        nfixAt = parmDict['nfixAt']
    2537         Tdata = atNo*[' ',]
    25382569        Xdata = np.zeros((3,atNo))
    2539         keys = {':Atype':Tdata,':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
     2570        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
    25402571        for iatm in range(nfixAt):
    25412572            for key in keys:
    25422573                parm = ':'+str(iatm)+key
    25432574                if parm in parmDict:
    2544                     if key == ':Atype':
    2545                         keys[key][iatm] = aTypes.index(parmDict[parm])
    2546                     else:
    2547                         keys[key][iatm] = parmDict[parm]
     2575                    keys[key][iatm] = parmDict[parm]
    25482576        iatm = nfixAt
    25492577        for iObj in range(parmDict['nObj']):
     
    25672595                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
    25682596                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
    2569                 Xdata.T[iatm:iatm+len(Cart)] = np.inner(prodQVQ(Qori,Cart),Bmat)+Pos
    2570                 for i in range(len(Cart)):
    2571                     Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
    2572                     iatm += 1
     2597                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat.T,prodQVQ(Qori,Cart)).T+Pos
     2598                iatm += len(Cart)
    25732599            elif parmDict[pfx+'Type'] == 'Atom':
    25742600                atNo = parmDict[pfx+'atNo']
     
    25762602                    parm = pfx+key[1:]              #remove extra ':'
    25772603                    if parm in parmDict:
    2578                         if key == ':Atype':
    2579                             keys[key][atNo] = aTypes.index(parmDict[parm])
    2580                         else:
    2581                             keys[key][atNo] = parmDict[parm]
     2604                        keys[key][atNo] = parmDict[parm]
    25822605                iatm += 1
    25832606            else:
    25842607                continue        #skips March Dollase
    2585         return Tdata,Xdata.T
     2608        return Xdata.T
    25862609       
    25872610    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
     
    26262649        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
    26272650        t0 = time.time()
    2628         Xdata = GetAtomTX(RBdata,parmDict)[1]                   #get new atom coords from RB
     2651        Xdata = GetAtomX(RBdata,parmDict)                   #get new atom coords from RB
    26292652        tsum += (time.time()-t0)
    26302653        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
     
    26442667        refList[6] = refList[4]-refList[5]
    26452668        M = np.inner(refList[6],np.inner(rcov,refList[6]))
    2646 #        print M,parmDict['sumFosq'],np.sum(refList[6]**2),np.sum(refList[4]**2)
    2647 #        print np.sum(refList[6]**2)/np.sum(refList[4]**2)
    26482669        return M/np.sum(refList[4]**2)
    26492670
     
    26942715    parmDict['nObj'] = len(MCSAObjs)
    26952716    aTypes = list(aTypes)
    2696     Tdata,Xdata = GetAtomTX(RBdata,parmDict)
     2717    Tdata = GetAtomT(RBdata,parmDict)
     2718    Xdata = GetAtomX(RBdata,parmDict)
    26972719    Mdata = GetAtomM(Xdata,SGData)
    26982720    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
Note: See TracChangeset for help on using the changeset viewer.