Changeset 883


Ignore:
Timestamp:
Apr 9, 2013 9:22:09 AM (9 years ago)
Author:
vondreele
Message:

begin RB refinement mods & fixed so non-RB refinements still work

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIconstrGUI.py

    r881 r883  
    101101    rigidbodyDict = G2frame.PatternTree.GetItemPyData(   
    102102        G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Rigid bodies'))
    103     rbVary,rbDict,rbIds = G2str.GetRigidBodyModels(rigidbodyDict,Print=False)
     103    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
     104    rbVary,rbDict = G2str.GetRigidBodyModels(rigidbodyDict,Print=False)
    104105    AtomDict = dict([Phases[phase]['pId'],Phases[phase]['Atoms']] for phase in Phases)
    105106    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable = G2str.GetPhaseData(Phases,rbIds=rbIds,Print=False)
  • trunk/GSASIIstruct.py

    r882 r883  
    136136    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
    137137    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
    138     varyList = phaseVary+hapVary+histVary
     138    varyList = rbVary+phaseVary+hapVary+histVary
    139139    constrDict,fixedList = GetConstraints(GPXfile)
    140140    return G2mv.CheckConstraints(varyList,constrDict,fixedList)
     
    339339    return GPXback
    340340
    341 def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
     341def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
    342342    ''' Updates gpxfile from all histograms that are found in any phase
    343     and any phase that used a histogram
     343    and any phase that used a histogram. Also updates rigid body definitions.
    344344    input:
    345345        GPXfile = .gpx full file name
    346346        Histograms = dictionary of histograms as {name:data,...}
    347347        Phases = dictionary of phases that use histograms
     348        RigidBodies = dictionary of rigid bodies
    348349        CovData = dictionary of refined variables, varyList, & covariance matrix
    349350        makeBack = True if new backup of .gpx file is to be made; else use the last one made
     
    369370        elif datum[0] == 'Covariance':
    370371            data[0][1] = CovData
     372        elif datum[0] == 'Rigid bodies':
     373            data[0][1] = RigidBodies
    371374        try:
    372375            histogram = Histograms[datum[0]]
     
    571574    rbDict = {}
    572575    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    573     if len(rigidbodyDict['Vector']):
    574         for irb,item in enumerate(rigidbodyDict['RBIds']['Vector']):
     576    if len(rbIds['Vector']):
     577        for irb,item in enumerate(rbIds['Vector']):
    575578            if rigidbodyDict['Vector'][item]['useCount']:
    576579                RBmags = rigidbodyDict['Vector'][item]['VectMag']
     
    584587                    print >>pFile,'\nVector rigid body model:'
    585588                    PrintVecRBModel(rigidbodyDict['Vector'][item])
    586     if len(rigidbodyDict['Residue']):
    587         for item in rigidbodyDict['RBIds']['Residue']:
     589    if len(rbIds['Residue']):
     590        for item in rbIds['Residue']:
    588591            if rigidbodyDict['Residue'][item]['useCount']:
    589592                if Print:
     
    592595    return rbVary,rbDict
    593596   
    594 def ApplyRBModels(parmDict,Phases,rigidbodyDict):
     597def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
     598   
     599    def PrintRBVectandSig(VectRB,VectSig):
     600        print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
     601        namstr = '  names :'
     602        valstr = '  values:'
     603        sigstr = '  esds  :'
     604        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
     605            namstr += '%12s'%('Vect '+str(i))
     606            valstr += '%12.4f'%(val)
     607            if sig:
     608                sigstr += '%12.4f'%(sig)
     609            else:
     610                sigstr += 12*' '
     611        print >>pFile,namstr
     612        print >>pFile,valstr
     613        print >>pFile,sigstr       
     614       
     615    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
     616    if not RBIds['Vector']:
     617        return
     618    for irb,item in enumerate(RBIds['Vector']):
     619        if rigidbodyDict['Vector'][item]['useCount']:
     620            VectSig = []
     621            RBmags = rigidbodyDict['Vector'][item]['VectMag']
     622            for i,mag in enumerate(RBmags):
     623                name = '::RBV;'+str(irb)+':'+str(i)
     624                mag = parmDict[name]
     625                if name in sigDict:
     626                    VectSig.append(sigDict[name])
     627            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
     628       
     629def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
    595630    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
    596631    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
     
    598633    atxIds = ['Ax:','Ay:','Az:']
    599634    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
    600     RBIds = rigidbodyDict['RBIds']
     635    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
    601636    if not RBIds['Vector'] and not RBIds['Residue']:
    602637        return
    603638    VRBIds = RBIds['Vector']
    604639    RRBIds = RBIds['Residue']
    605     RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
     640    if Update:
     641        RBData = rigidbodyDict
     642    else:
     643        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
    606644    if RBIds['Vector']:                       # first update the vector magnitudes
    607645        VRBData = RBData['Vector']
     
    617655        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
    618656        pfx = str(Phase['pId'])+'::'
    619         RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
    620         for irb,rbId in enumerate(VRBIds):
    621             RBObj = RBModels['Vector'][irb]
     657        if Update:
     658            RBModels = Phase['RBModels']
     659        else:
     660            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
     661        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
     662            jrb = VRBIds.index(RBObj['RBId'])
     663            rbsx = str(irb)+':'+str(jrb)
     664            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
     665                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
     666            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
     667                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
     668            TLS = RBObj['ThermalMotion']
     669            if 'T' in TLS[0]:
     670                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
     671                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
     672            if 'L' in TLS[0]:
     673                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
     674                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
     675            if 'S' in TLS[0]:
     676                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
     677                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
     678            if 'U' in TLS[0]:
     679                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBVU:'+rbsx]
    622680            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
    623681            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
     
    631689                elif UIJ[0] == 'I':
    632690                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = Uij[j+1]
    633         for irb,rbId in enumerate(RRBIds):
    634             RBObj = RBModels['Residue'][irb]           
     691        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
     692            jrb = RRBIds.index(RBObj['RBId'])
     693            rbsx = str(irb)+':'+str(jrb)
     694            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
     695                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
     696            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
     697                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
     698            TLS = RBObj['ThermalMotion']
     699            if 'T' in TLS[0]:
     700                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
     701                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
     702            if 'L' in TLS[0]:
     703                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
     704                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
     705            if 'S' in TLS[0]:
     706                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
     707                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
     708            if 'U' in TLS[0]:
     709                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
     710            for itors,tors in enumerate(RBObj['Torsions']):
     711                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
    635712            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
    636713            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
     
    791868        print >>pFile,ptstr
    792869       
    793     def MakeRBParms(rbKey):
     870    def MakeRBParms(rbKey,phaseVary,phaseDict):
    794871        rbid = str(rbids.index(RB['RBId']))
    795872        pfxRB = pfx+'RB'+rbKey+'P'
     
    805882            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
    806883            phaseDict[name] = RB['Orient'][0][i]
    807             if RB['Orient'][1] == 'V' and i:
     884            if RB['Orient'][1] == 'AV' and i:
    808885                phaseVary += [name,]
    809886            elif RB['Orient'][1] == 'A' and not i:
    810887                phaseVary += [name,]
    811888           
    812     def MakeRBThermals(rbKey):
     889    def MakeRBThermals(rbKey,phaseVary,phaseDict):
    813890        rbid = str(rbids.index(RB['RBId']))
    814891        tlstr = ['11','22','33','12','13','23']
     
    830907        if 'S' in RB['ThermalMotion'][0]:
    831908            pfxRB = pfx+'RB'+rbKey+'S'
    832             for i in range(5):
     909            for i in range(8):
    833910                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
    834911                phaseDict[name] = RB['ThermalMotion'][1][i+12]
     
    841918                phaseVary += [name,]
    842919               
    843     def MakeRBTorsions(rbKey):
     920    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
    844921        rbid = str(rbids.index(RB['RBId']))
    845922        pfxRB = pfx+'RB'+rbKey+'Tr;'
     
    887964            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
    888965            for iRB,RB in enumerate(resRBData):
    889                 MakeRBParms('R')
    890                 MakeRBThermals('R')
    891                 MakeRBTorsions('R')
     966                MakeRBParms('R',phaseVary,phaseDict)
     967                MakeRBThermals('R',phaseVary,phaseDict)
     968                MakeRBTorsions('R',phaseVary,phaseDict)
    892969       
    893970        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
     
    895972            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
    896973            for iRB,RB in enumerate(vecRBData):
    897                 MakeRBParms('V')
    898                 MakeRBThermals('V')
     974                MakeRBParms('V',phaseVary,phaseDict)
     975                MakeRBThermals('V',phaseVary,phaseDict)
    899976                   
    900977        Natoms[pfx] = 0
     
    11751252    return cellSig           
    11761253   
    1177 def SetPhaseData(parmDict,sigDict,Phases,covData,RestraintDict=None,pFile=None):
     1254def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
    11781255   
    11791256    def PrintAtomsAndSig(General,Atoms,atomsSig):
     
    12241301                print >>pFile,sigstr
    12251302               
    1226     #def PrintRBObjectsAndSig()
     1303    def PrintRBObjPOAndSig(rbfx,rbsx):
     1304        namstr = '  names :'
     1305        valstr = '  values:'
     1306        sigstr = '  esds  :'
     1307        for i,px in enumerate(['Px:','Py:','Pz:']):
     1308            name = pfx+rbfx+px+rbsx
     1309            namstr += '%12s'%('Pos '+px[1])
     1310            valstr += '%12.5f'%(parmDict[name])
     1311            if name in sigDict:
     1312                sigstr += '%12.5f'%(sigDict[name])
     1313            else:
     1314                sigstr += 12*' '
     1315        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
     1316            name = pfx+rbfx+po+rbsx
     1317            namstr += '%12s'%('Ori '+po[1])
     1318            valstr += '%12.5f'%(parmDict[name])
     1319            if name in sigDict:
     1320                sigstr += '%12.5f'%(sigDict[name])
     1321            else:
     1322                sigstr += 12*' '
     1323        print >>pFile,namstr
     1324        print >>pFile,valstr
     1325        print >>pFile,sigstr
     1326       
     1327    def PrintRBObjTLSAndSig(rbfx,rbsx):
     1328        namstr = '  names :'
     1329        valstr = '  values:'
     1330        sigstr = '  esds  :'
     1331        for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
     1332            name = pfx+rbfx+pt+rbsx
     1333            namstr += '%12s'%(pt[:3])
     1334            valstr += '%12.5f'%(parmDict[name])
     1335            if name in sigDict:
     1336                sigstr += '%12.5f'%(sigDict[name])
     1337            else:
     1338                sigstr += 12*' '
     1339        print >>pFile,namstr
     1340        print >>pFile,valstr
     1341        print >>pFile,sigstr
     1342        namstr = '  names :'
     1343        valstr = '  values:'
     1344        sigstr = '  esds  :'
     1345        for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
     1346            name = pfx+rbfx+pt+rbsx
     1347            namstr += '%12s'%(pt[:3])
     1348            valstr += '%12.3f'%(parmDict[name])
     1349            if name in sigDict:
     1350                sigstr += '%12.3f'%(sigDict[name])
     1351            else:
     1352                sigstr += 12*' '
     1353        print >>pFile,namstr
     1354        print >>pFile,valstr
     1355        print >>pFile,sigstr
     1356        namstr = '  names :'
     1357        valstr = '  values:'
     1358        sigstr = '  esds  :'
     1359        for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
     1360            name = pfx+rbfx+pt+rbsx
     1361            namstr += '%12s'%(pt[:3])
     1362            valstr += '%12.3f'%(parmDict[name])
     1363            if name in sigDict:
     1364                sigstr += '%12.3f'%(sigDict[name])
     1365            else:
     1366                sigstr += 12*' '
     1367        print >>pFile,namstr
     1368        print >>pFile,valstr
     1369        print >>pFile,sigstr
     1370       
     1371    def PrintRBObjTorAndSig(rbsx):
     1372        namstr = '  names :'
     1373        valstr = '  values:'
     1374        sigstr = '  esds  :'       
    12271375               
    12281376    def PrintSHtextureAndSig(textureData,SHtextureSig):
     
    12571405        print >>pFile,ptstr
    12581406        print >>pFile,sigstr
    1259        
    12601407           
    12611408    print >>pFile,'\n Phases:'
     
    13191466                    refl[7] = 0
    13201467        else:
    1321             #new RBObj parms here, mod atoms, & PrintRBObjectsAndSig
     1468            VRBIds = RBIds['Vector']
     1469            RRBIds = RBIds['Residue']
     1470            RBModels = Phase['RBModels']
     1471            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
     1472                jrb = VRBIds.index(RBObj['RBId'])
     1473                rbsx = str(irb)+':'+str(jrb)
     1474                print >>pFile,' Vector rigid body '
     1475                PrintRBObjPOAndSig('RBV',rbsx)
     1476                PrintRBObjTLSAndSig('RBV',rbsx)
     1477            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
     1478                jrb = VRBIds.index(RBObj['RBId'])
     1479                rbsx = str(irb)+':'+str(jrb)
     1480                print >>pFile,' Residue rigid body '
     1481                PrintRBObjPOAndSig('RBR',rbsx)
     1482                PrintRBObjTLSAndSig('RBR',rbsx)
     1483                PrintRBObjTorAndSig(rbsx)
    13221484            atomsSig = {}
    13231485            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
     
    27092871        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
    27102872        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
     2873# or else do RB modification of dFdvDict here? Perhaps not...
    27112874    return dFdvDict
    27122875   
     
    36423805            refList = Histogram['Data']
    36433806            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
     3807# do RB modification of dFdvDict here; varylist will contain RB variables so dFdvDict needs corresponding entries
    36443808            dMdvh = np.zeros((len(varylist),len(refList)))
    36453809            for iref,ref in enumerate(refList):
     
    37323896            refList = Histogram['Data']
    37333897            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
     3898# do RB modification of dFdvDict here; varylist will contain RB variables so dFdvDict needs corresponding entries
    37343899            dMdvh = np.zeros((len(varylist),len(refList)))
    37353900            wdf = np.zeros(len(refList))
     
    39884153        Values2Dict(parmDict, varyList, result[0])
    39894154        G2mv.Dict2Map(parmDict,varyList)
    3990        
    39914155        Rvals['Nobs'] = Histograms['Nobs']
    39924156        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
     
    40374201    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
    40384202    G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
    4039     SetPhaseData(parmDict,sigDict,Phases,covData,restraintDict,printFile)
     4203    ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
     4204    SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
     4205    SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
    40404206    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
    40414207    SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
    4042     SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
     4208    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData)
    40434209    printFile.close()
    40444210    print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
     
    42244390            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
    42254391        # add the uncertainties into the esd dictionary (sigDict)
     4392        ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
     4393#??        SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
    42264394        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
    42274395        SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
    42284396        SeqResult[histogram] = covData
    4229         SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
     4397        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,covData,makeBack)
    42304398        makeBack = False
    42314399    SetSeqResult(GPXfile,Histograms,SeqResult)
Note: See TracChangeset for help on using the changeset viewer.