Changeset 686

Timestamp:

Jul 11, 2012 4:01:20 PM (11 years ago)

Author:

vondreele

Message:

add restraints to GSASIIstruct.py - input & printing only
no calcs yet

Location:

trunk

Files:

• GSASII.py (modified) (1 diff)
• GSASIIstruct.py (modified) (14 diffs)

trunk/GSASII.py

@@ -2039 +2039 @@
         Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
         print Histograms.keys()
-        Natoms,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable = G2str.GetPhaseData(Phases,Print=False)        
+        Natoms,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable = G2str.GetPhaseData(Phases,RestDict=None,Print=False)        
         hapVary,hapDict,controlDict = G2str.GetHistogramPhaseData(Phases,Histograms,Print=False)
         histVary,histDict,controlDict = G2str.GetHistogramData(Histograms,Print=False)

trunk/GSASIIstruct.py

@@ -128 +128 @@
     if not Histograms:
         return 'Error: no diffraction data',''
-    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,Print=False)
+    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestDict=None,Print=False)
     hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
     histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
@@ -134 +134 @@
     constrDict,fixedList = GetConstraints(GPXfile)
     return G2mv.CheckConstraints(varyList,constrDict,fixedList)
+    
+def GetRestraints(GPXfile):
+    '''Read the restraints from the GPX file
+    '''
+    fl = open(GPXfile,'rb')
+    while True:
+        try:
+            data = cPickle.load(fl)
+        except EOFError:
+            break
+        datum = data[0]
+        if datum[0] == 'Restraints':
+            restDict = datum[1]
+    fl.close()
+    return restDict
     
 def GetPhaseNames(GPXfile):
@@ -560 +575 @@
 ################################################################################        
                     
-def GetPhaseData(PhaseData,Print=True):
+def GetPhaseData(PhaseData,RestDict=None,Print=True):
             
     def PrintFFtable(FFtable):
@@ -633 +648 @@
             ptstr += '%12.4f'%(SHcoeff[item]) 
         print ptlbls
-        print ptstr    
+        print ptstr
+        
+    def PrintRestraints(phaseRest):
+        if phaseRest:
+            print '\n Restraints:'
+            names = ['Bonds','Angles','Planes','Volumes']
+            for i,name in enumerate(['Bond','Angle','Plane','Chiral']):
+                itemRest = phaseRest[name]
+                if itemRest[names[i]]:
+                    print '\n  %30s %10.3f'%(name+' restraint weight factor',itemRest['wtFactor'])
+                    print '  atoms(symOp), calc, obs, sig: '
+                    for item in phaseRest[name][names[i]]:
+                        text = '   '
+                        for a,at in enumerate(item[0]):
+                            text += at+'+('+item[1][a]+') '
+                            if (a+1)%5 == 0:
+                                text += '\n   '
+                        print text,' %.3f %.3f %.3f'%(item[3],item[4],item[5])
         
     if Print: print ' Phases:'
@@ -741 +773 @@
                     '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
                 PrintTexture(textureData)
+                if name in RestDict:
+                    PrintRestraints(RestDict[name])
                     
         elif PawleyRef:
@@ -2915 +2949 @@
     parmdict.update(zip(varylist,values))
     G2mv.Dict2Map(parmdict,varylist)
-    Histograms,Phases = HistoPhases
+    Histograms,Phases,restDict = HistoPhases
     nvar = len(varylist)
     dMdv = np.empty(0)
@@ -2980 +3014 @@
     parmdict.update(zip(varylist,values))
     G2mv.Dict2Map(parmdict,varylist)
-    Histograms,Phases = HistoPhases
+    Histograms,Phases,restDict = HistoPhases
     nvar = len(varylist)
     Hess = np.empty(0)
@@ -3066 +3100 @@
     Values2Dict(parmdict, varylist, values)
     G2mv.Dict2Map(parmdict,varylist)
-    Histograms,Phases = HistoPhases
+    Histograms,Phases,restDict = HistoPhases
     M = np.empty(0)
     SumwYo = 0
@@ -3184 +3218 @@
     calcControls.update(Controls)            
     constrDict,fixedList = GetConstraints(GPXfile)
+    restDict = GetRestraints(GPXfile)
     Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
     if not Phases:
@@ -3193 +3228 @@
         print ' *** Refine aborted ***'
         raise Exception        
-    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
+    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restDict)
     calcControls['Natoms'] = Natoms
     calcControls['FFtables'] = FFtables
@@ -3237 +3272 @@
             result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
                 ftol=Ftol,col_deriv=True,factor=Factor,
-                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+                args=([Histograms,Phases,restDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
             ncyc = int(result[2]['nfev']/2)
         elif 'Hessian' in Controls['deriv type']:
             result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
-                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+                args=([Histograms,Phases,restDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
             ncyc = result[2]['num cyc']+1
             Rvals['lamMax'] = result[2]['lamMax']
         else:           #'numeric'
             result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
-                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+                args=([Histograms,Phases,restDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
             ncyc = int(result[2]['nfev']/len(varyList))
 #        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
@@ -3336 +3371 @@
     ShowControls(Controls)            
     constrDict,fixedList = GetConstraints(GPXfile)
+    restDict = GetRestraints(GPXfile)
     Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
     if not Phases:
@@ -3345 +3381 @@
         print ' *** Refine aborted ***'
         raise Exception
-    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
+    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restDict,False)
     if 'Seq Data' in Controls:
         histNames = Controls['Seq Data']
@@ -3423 +3459 @@
                 result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
                     ftol=Ftol,col_deriv=True,factor=Factor,
-                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+                    args=([Histo,Phases,restDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
                 ncyc = int(result[2]['nfev']/2)
             elif 'Hessian' in Controls['deriv type']:
                 result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
-                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+                    args=([Histo,Phases,restDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
                 ncyc = result[2]['num cyc']+1                           
             else:           #'numeric'
                 result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
-                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
+                    args=([Histo,Phases,restDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
                 ncyc = int(result[2]['nfev']/len(varyList))