| 1 | # -*- coding: utf-8 -*- |
|---|
| 2 | ''' |
|---|
| 3 | *GSASIIstrMain: main structure routine* |
|---|
| 4 | --------------------------------------- |
|---|
| 5 | |
|---|
| 6 | ''' |
|---|
| 7 | ########### SVN repository information ################### |
|---|
| 8 | # $Date: 2015-06-18 18:25:47 +0000 (Thu, 18 Jun 2015) $ |
|---|
| 9 | # $Author: vondreele $ |
|---|
| 10 | # $Revision: 1896 $ |
|---|
| 11 | # $URL: trunk/GSASIIstrMain.py $ |
|---|
| 12 | # $Id: GSASIIstrMain.py 1896 2015-06-18 18:25:47Z vondreele $ |
|---|
| 13 | ########### SVN repository information ################### |
|---|
| 14 | import sys |
|---|
| 15 | import os |
|---|
| 16 | import os.path as ospath |
|---|
| 17 | import time |
|---|
| 18 | import math |
|---|
| 19 | import copy |
|---|
| 20 | import random |
|---|
| 21 | import cPickle |
|---|
| 22 | import numpy as np |
|---|
| 23 | import numpy.ma as ma |
|---|
| 24 | import numpy.linalg as nl |
|---|
| 25 | import scipy.optimize as so |
|---|
| 26 | import GSASIIpath |
|---|
| 27 | GSASIIpath.SetVersionNumber("$Revision: 1896 $") |
|---|
| 28 | import GSASIIlattice as G2lat |
|---|
| 29 | import GSASIIspc as G2spc |
|---|
| 30 | import GSASIImapvars as G2mv |
|---|
| 31 | import GSASIImath as G2mth |
|---|
| 32 | import GSASIIstrIO as G2stIO |
|---|
| 33 | import GSASIIstrMath as G2stMth |
|---|
| 34 | import GSASIIobj as G2obj |
|---|
| 35 | |
|---|
| 36 | sind = lambda x: np.sin(x*np.pi/180.) |
|---|
| 37 | cosd = lambda x: np.cos(x*np.pi/180.) |
|---|
| 38 | tand = lambda x: np.tan(x*np.pi/180.) |
|---|
| 39 | asind = lambda x: 180.*np.arcsin(x)/np.pi |
|---|
| 40 | acosd = lambda x: 180.*np.arccos(x)/np.pi |
|---|
| 41 | atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
|---|
| 42 | |
|---|
| 43 | ateln2 = 8.0*math.log(2.0) |
|---|
| 44 | DEBUG = True |
|---|
| 45 | |
|---|
| 46 | def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList, |
|---|
| 47 | calcControls,pawleyLookup,ifPrint,printFile,dlg): |
|---|
| 48 | 'Core optimization routines, shared between SeqRefine and Refine' |
|---|
| 49 | # print 'current',varyList |
|---|
| 50 | # for item in parmDict: print item,parmDict[item] ######### show dict just before refinement |
|---|
| 51 | G2mv.Map2Dict(parmDict,varyList) |
|---|
| 52 | Rvals = {} |
|---|
| 53 | while True: |
|---|
| 54 | begin = time.time() |
|---|
| 55 | values = np.array(G2stMth.Dict2Values(parmDict, varyList)) |
|---|
| 56 | # test code to compute GOF and save for external repeat |
|---|
| 57 | #args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg) |
|---|
| 58 | #print '*** before fit chi**2',np.sum(G2stMth.errRefine(values,*args)**2) |
|---|
| 59 | #fl = open('beforeFit.cpickle','wb') |
|---|
| 60 | #import cPickle |
|---|
| 61 | #cPickle.dump(values,fl,1) |
|---|
| 62 | #cPickle.dump(args[:-1],fl,1) |
|---|
| 63 | #fl.close() |
|---|
| 64 | Ftol = Controls['min dM/M'] |
|---|
| 65 | Factor = Controls['shift factor'] |
|---|
| 66 | if 'Jacobian' in Controls['deriv type']: |
|---|
| 67 | result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True, |
|---|
| 68 | ftol=Ftol,col_deriv=True,factor=Factor, |
|---|
| 69 | args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
|---|
| 70 | ncyc = int(result[2]['nfev']/2) |
|---|
| 71 | elif 'Hessian' in Controls['deriv type']: |
|---|
| 72 | maxCyc = Controls['max cyc'] |
|---|
| 73 | result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,Print=ifPrint, |
|---|
| 74 | args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
|---|
| 75 | ncyc = result[2]['num cyc']+1 |
|---|
| 76 | Rvals['lamMax'] = result[2]['lamMax'] |
|---|
| 77 | else: #'numeric' |
|---|
| 78 | result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor, |
|---|
| 79 | args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
|---|
| 80 | ncyc = 1 |
|---|
| 81 | if len(varyList): |
|---|
| 82 | ncyc = int(result[2]['nfev']/len(varyList)) |
|---|
| 83 | # table = dict(zip(varyList,zip(values,result[0],(result[0]-values)))) |
|---|
| 84 | # for item in table: print item,table[item] #useful debug - are things shifting? |
|---|
| 85 | runtime = time.time()-begin |
|---|
| 86 | Rvals['converged'] = result[2].get('Converged') |
|---|
| 87 | Rvals['DelChi2'] = result[2].get('DelChi2',-1.) |
|---|
| 88 | Rvals['chisq'] = np.sum(result[2]['fvec']**2) |
|---|
| 89 | G2stMth.Values2Dict(parmDict, varyList, result[0]) |
|---|
| 90 | G2mv.Dict2Map(parmDict,varyList) |
|---|
| 91 | Rvals['Nobs'] = Histograms['Nobs'] |
|---|
| 92 | Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to % |
|---|
| 93 | Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList))) |
|---|
| 94 | print >>printFile,' Number of function calls:',result[2]['nfev'], \ |
|---|
| 95 | ' No. of observations: ',Histograms['Nobs'],' No. of parameters: ',len(varyList), \ |
|---|
| 96 | ' User rejected: ',Histograms['Nrej'],' Sp. gp. extinct: ',Histograms['Next'] |
|---|
| 97 | print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc) |
|---|
| 98 | print >>printFile,' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']) |
|---|
| 99 | IfOK = True |
|---|
| 100 | try: |
|---|
| 101 | covMatrix = result[1]*Rvals['GOF']**2 |
|---|
| 102 | sig = np.sqrt(np.diag(covMatrix)) |
|---|
| 103 | if np.any(np.isnan(sig)): |
|---|
| 104 | print '*** Least squares aborted - some invalid esds possible ***' |
|---|
| 105 | # table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig))) |
|---|
| 106 | # for item in table: print item,table[item] #useful debug - are things shifting? |
|---|
| 107 | break #refinement succeeded - finish up! |
|---|
| 108 | except TypeError,FloatingPointError: #result[1] is None on singular matrix |
|---|
| 109 | IfOK = False |
|---|
| 110 | if not len(varyList): |
|---|
| 111 | covMatrix = [] |
|---|
| 112 | sig = [] |
|---|
| 113 | break |
|---|
| 114 | print '**** Refinement failed - singular matrix ****' |
|---|
| 115 | if 'Hessian' in Controls['deriv type']: |
|---|
| 116 | num = len(varyList)-1 |
|---|
| 117 | for i,val in enumerate(np.flipud(result[2]['psing'])): |
|---|
| 118 | if val: |
|---|
| 119 | print 'Removing parameter: ',varyList[num-i] |
|---|
| 120 | del(varyList[num-i]) |
|---|
| 121 | else: |
|---|
| 122 | Ipvt = result[2]['ipvt'] |
|---|
| 123 | for i,ipvt in enumerate(Ipvt): |
|---|
| 124 | if not np.sum(result[2]['fjac'],axis=1)[i]: |
|---|
| 125 | print 'Removing parameter: ',varyList[ipvt-1] |
|---|
| 126 | del(varyList[ipvt-1]) |
|---|
| 127 | break |
|---|
| 128 | G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls) |
|---|
| 129 | return IfOK,Rvals,result,covMatrix,sig |
|---|
| 130 | |
|---|
| 131 | def Refine(GPXfile,dlg): |
|---|
| 132 | 'Global refinement -- refines to minimize against all histograms' |
|---|
| 133 | import pytexture as ptx |
|---|
| 134 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
|---|
| 135 | |
|---|
| 136 | printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w') |
|---|
| 137 | G2stIO.ShowBanner(printFile) |
|---|
| 138 | varyList = [] |
|---|
| 139 | parmDict = {} |
|---|
| 140 | G2mv.InitVars() |
|---|
| 141 | Controls = G2stIO.GetControls(GPXfile) |
|---|
| 142 | G2stIO.ShowControls(Controls,printFile) |
|---|
| 143 | calcControls = {} |
|---|
| 144 | calcControls.update(Controls) |
|---|
| 145 | constrDict,fixedList = G2stIO.GetConstraints(GPXfile) |
|---|
| 146 | restraintDict = G2stIO.GetRestraints(GPXfile) |
|---|
| 147 | Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile) |
|---|
| 148 | if not Phases: |
|---|
| 149 | print ' *** ERROR - you have no phases to refine! ***' |
|---|
| 150 | print ' *** Refine aborted ***' |
|---|
| 151 | return False,'No phases' |
|---|
| 152 | if not Histograms: |
|---|
| 153 | print ' *** ERROR - you have no data to refine with! ***' |
|---|
| 154 | print ' *** Refine aborted ***' |
|---|
| 155 | return False,'No data' |
|---|
| 156 | rigidbodyDict = G2stIO.GetRigidBodies(GPXfile) |
|---|
| 157 | rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) |
|---|
| 158 | rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile) |
|---|
| 159 | Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile) |
|---|
| 160 | calcControls['atomIndx'] = atomIndx |
|---|
| 161 | calcControls['Natoms'] = Natoms |
|---|
| 162 | calcControls['FFtables'] = FFtables |
|---|
| 163 | calcControls['BLtables'] = BLtables |
|---|
| 164 | calcControls['maxSSwave'] = maxSSwave |
|---|
| 165 | hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile) |
|---|
| 166 | TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary) |
|---|
| 167 | constrDict += TwConstr |
|---|
| 168 | fixedList += TwFixed |
|---|
| 169 | calcControls.update(controlDict) |
|---|
| 170 | histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile) |
|---|
| 171 | calcControls.update(controlDict) |
|---|
| 172 | varyList = rbVary+phaseVary+hapVary+histVary |
|---|
| 173 | parmDict.update(rbDict) |
|---|
| 174 | parmDict.update(phaseDict) |
|---|
| 175 | parmDict.update(hapDict) |
|---|
| 176 | parmDict.update(histDict) |
|---|
| 177 | G2stIO.GetFprime(calcControls,Histograms) |
|---|
| 178 | # do constraint processing |
|---|
| 179 | varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed |
|---|
| 180 | try: |
|---|
| 181 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
|---|
| 182 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict) |
|---|
| 183 | #print G2mv.VarRemapShow(varyList) |
|---|
| 184 | #print 'DependentVars',G2mv.GetDependentVars() |
|---|
| 185 | #print 'IndependentVars',G2mv.GetIndependentVars() |
|---|
| 186 | except: |
|---|
| 187 | print ' *** ERROR - your constraints are internally inconsistent ***' |
|---|
| 188 | #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList) |
|---|
| 189 | #print 'Errors',errmsg |
|---|
| 190 | #if warnmsg: print 'Warnings',warnmsg |
|---|
| 191 | return False,' Constraint error' |
|---|
| 192 | # print G2mv.VarRemapShow(varyList) |
|---|
| 193 | |
|---|
| 194 | ifPrint = True |
|---|
| 195 | print >>printFile,'\n Refinement results:' |
|---|
| 196 | print >>printFile,135*'-' |
|---|
| 197 | try: |
|---|
| 198 | IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict, |
|---|
| 199 | rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg) |
|---|
| 200 | sigDict = dict(zip(varyList,sig)) |
|---|
| 201 | newCellDict = G2stMth.GetNewCellParms(parmDict,varyList) |
|---|
| 202 | newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList) |
|---|
| 203 | covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals, |
|---|
| 204 | 'varyListStart':varyListStart, |
|---|
| 205 | 'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict, |
|---|
| 206 | 'newCellDict':newCellDict,'freshCOV':True} |
|---|
| 207 | # add the uncertainties into the esd dictionary (sigDict) |
|---|
| 208 | sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict)) |
|---|
| 209 | G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile) |
|---|
| 210 | G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True) |
|---|
| 211 | G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile) |
|---|
| 212 | G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile) |
|---|
| 213 | G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile) |
|---|
| 214 | G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile) |
|---|
| 215 | G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData) |
|---|
| 216 | printFile.close() |
|---|
| 217 | print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst' |
|---|
| 218 | print ' ***** Refinement successful *****' |
|---|
| 219 | except G2obj.G2Exception,Msg: |
|---|
| 220 | printFile.close() |
|---|
| 221 | return False,Msg.msg |
|---|
| 222 | |
|---|
| 223 | #for testing purposes!!! |
|---|
| 224 | if DEBUG: |
|---|
| 225 | #needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup |
|---|
| 226 | import cPickle |
|---|
| 227 | fl = open('testDeriv.dat','wb') |
|---|
| 228 | cPickle.dump(result[0],fl,1) |
|---|
| 229 | cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1) |
|---|
| 230 | cPickle.dump([constrDict,fixedList],fl,1) |
|---|
| 231 | cPickle.dump(parmDict,fl,1) |
|---|
| 232 | cPickle.dump(varyList,fl,1) |
|---|
| 233 | cPickle.dump(calcControls,fl,1) |
|---|
| 234 | cPickle.dump(pawleyLookup,fl,1) |
|---|
| 235 | fl.close() |
|---|
| 236 | if dlg: |
|---|
| 237 | return True,Rvals['Rwp'] |
|---|
| 238 | |
|---|
| 239 | def SeqRefine(GPXfile,dlg): |
|---|
| 240 | '''Perform a sequential refinement -- cycles through all selected histgrams, |
|---|
| 241 | one at a time |
|---|
| 242 | ''' |
|---|
| 243 | import pytexture as ptx |
|---|
| 244 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
|---|
| 245 | |
|---|
| 246 | printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w') |
|---|
| 247 | print 'Starting Sequential Refinement' |
|---|
| 248 | G2stIO.ShowBanner(printFile) |
|---|
| 249 | Controls = G2stIO.GetControls(GPXfile) |
|---|
| 250 | G2stIO.ShowControls(Controls,printFile,SeqRef=True) |
|---|
| 251 | restraintDict = G2stIO.GetRestraints(GPXfile) |
|---|
| 252 | Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile) |
|---|
| 253 | if not Phases: |
|---|
| 254 | print ' *** ERROR - you have no phases to refine! ***' |
|---|
| 255 | print ' *** Refine aborted ***' |
|---|
| 256 | return False,'No phases' |
|---|
| 257 | if not Histograms: |
|---|
| 258 | print ' *** ERROR - you have no data to refine with! ***' |
|---|
| 259 | print ' *** Refine aborted ***' |
|---|
| 260 | return False,'No data' |
|---|
| 261 | rigidbodyDict = G2stIO.GetRigidBodies(GPXfile) |
|---|
| 262 | rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]}) |
|---|
| 263 | rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile) |
|---|
| 264 | Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True) |
|---|
| 265 | for item in phaseVary: |
|---|
| 266 | if '::A0' in item: |
|---|
| 267 | print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****' |
|---|
| 268 | print '**** instead use the Dij parameters for each powder histogram ****' |
|---|
| 269 | return False,'Lattice parameter refinement error - see console message' |
|---|
| 270 | if '::C(' in item: |
|---|
| 271 | print '**** WARNING - phase texture parameters should not be refined in a sequential refinement ****' |
|---|
| 272 | print '**** instead use the C(L,N) parameters for each powder histogram ****' |
|---|
| 273 | return False,'Phase texture refinement error - see console message' |
|---|
| 274 | if 'Seq Data' in Controls: |
|---|
| 275 | histNames = Controls['Seq Data'] |
|---|
| 276 | else: |
|---|
| 277 | histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',]) |
|---|
| 278 | if 'Reverse Seq' in Controls: |
|---|
| 279 | if Controls['Reverse Seq']: |
|---|
| 280 | histNames.reverse() |
|---|
| 281 | SeqResult = {'histNames':histNames} |
|---|
| 282 | makeBack = True |
|---|
| 283 | Histo = {} |
|---|
| 284 | NewparmDict = {} |
|---|
| 285 | for ihst,histogram in enumerate(histNames): |
|---|
| 286 | print('Refining with '+str(histogram)) |
|---|
| 287 | ifPrint = False |
|---|
| 288 | if dlg: |
|---|
| 289 | dlg.SetTitle('Residual for histogram '+str(ihst)) |
|---|
| 290 | calcControls = {} |
|---|
| 291 | calcControls['atomIndx'] = atomIndx |
|---|
| 292 | calcControls['Natoms'] = Natoms |
|---|
| 293 | calcControls['FFtables'] = FFtables |
|---|
| 294 | calcControls['BLtables'] = BLtables |
|---|
| 295 | calcControls['maxSSwave'] = maxSSwave |
|---|
| 296 | Histo = {histogram:Histograms[histogram],} |
|---|
| 297 | hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False) |
|---|
| 298 | calcControls.update(controlDict) |
|---|
| 299 | histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False) |
|---|
| 300 | calcControls.update(controlDict) |
|---|
| 301 | varyList = rbVary+phaseVary+hapVary+histVary |
|---|
| 302 | if not ihst: |
|---|
| 303 | # save the initial vary list, but without histogram numbers on parameters |
|---|
| 304 | saveVaryList = varyList[:] |
|---|
| 305 | for i,item in enumerate(saveVaryList): |
|---|
| 306 | items = item.split(':') |
|---|
| 307 | if items[1]: |
|---|
| 308 | items[1] = '' |
|---|
| 309 | item = ':'.join(items) |
|---|
| 310 | saveVaryList[i] = item |
|---|
| 311 | SeqResult['varyList'] = saveVaryList |
|---|
| 312 | origvaryList = varyList[:] |
|---|
| 313 | parmDict = {} |
|---|
| 314 | parmDict.update(phaseDict) |
|---|
| 315 | parmDict.update(hapDict) |
|---|
| 316 | parmDict.update(histDict) |
|---|
| 317 | if Controls['Copy2Next']: |
|---|
| 318 | parmDict.update(NewparmDict) |
|---|
| 319 | G2stIO.GetFprime(calcControls,Histo) |
|---|
| 320 | # do constraint processing |
|---|
| 321 | #reload(G2mv) # debug |
|---|
| 322 | G2mv.InitVars() |
|---|
| 323 | constrDict,fixedList = G2stIO.GetConstraints(GPXfile) |
|---|
| 324 | varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed |
|---|
| 325 | try: |
|---|
| 326 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
|---|
| 327 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst) |
|---|
| 328 | constraintInfo = (groups,parmlist,constrDict,fixedList,ihst) |
|---|
| 329 | except: |
|---|
| 330 | print ' *** ERROR - your constraints are internally inconsistent ***' |
|---|
| 331 | #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList) |
|---|
| 332 | #print 'Errors',errmsg |
|---|
| 333 | #if warnmsg: print 'Warnings',warnmsg |
|---|
| 334 | return False,' Constraint error' |
|---|
| 335 | #print G2mv.VarRemapShow(varyList) |
|---|
| 336 | if not ihst: |
|---|
| 337 | # first histogram to refine against |
|---|
| 338 | firstVaryList = [] |
|---|
| 339 | for item in varyList: |
|---|
| 340 | items = item.split(':') |
|---|
| 341 | if items[1]: |
|---|
| 342 | items[1] = '' |
|---|
| 343 | item = ':'.join(items) |
|---|
| 344 | firstVaryList.append(item) |
|---|
| 345 | newVaryList = firstVaryList |
|---|
| 346 | else: |
|---|
| 347 | newVaryList = [] |
|---|
| 348 | for item in varyList: |
|---|
| 349 | items = item.split(':') |
|---|
| 350 | if items[1]: |
|---|
| 351 | items[1] = '' |
|---|
| 352 | item = ':'.join(items) |
|---|
| 353 | newVaryList.append(item) |
|---|
| 354 | if newVaryList != firstVaryList and Controls['Copy2Next']: |
|---|
| 355 | # variable lists are expected to match between sequential refinements when Copy2Next is on |
|---|
| 356 | print '**** ERROR - variable list for this histogram does not match previous' |
|---|
| 357 | print ' Copy of variables is not possible' |
|---|
| 358 | print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables' |
|---|
| 359 | combined = list(set(firstVaryList+newVaryList)) |
|---|
| 360 | c = [var for var in combined if var not in newVaryList] |
|---|
| 361 | p = [var for var in combined if var not in firstVaryList] |
|---|
| 362 | line = 'Variables in previous but not in current: ' |
|---|
| 363 | if c: |
|---|
| 364 | for var in c: |
|---|
| 365 | if len(line) > 100: |
|---|
| 366 | print line |
|---|
| 367 | line = ' ' |
|---|
| 368 | line += var + ', ' |
|---|
| 369 | else: |
|---|
| 370 | line += 'none' |
|---|
| 371 | print line |
|---|
| 372 | print '\nPrevious refinement has',len(firstVaryList),'variables' |
|---|
| 373 | line = 'Variables in current but not in previous: ' |
|---|
| 374 | if p: |
|---|
| 375 | for var in p: |
|---|
| 376 | if len(line) > 100: |
|---|
| 377 | print line |
|---|
| 378 | line = ' ' |
|---|
| 379 | line += var + ', ' |
|---|
| 380 | else: |
|---|
| 381 | line += 'none' |
|---|
| 382 | print line |
|---|
| 383 | return False,line |
|---|
| 384 | |
|---|
| 385 | ifPrint = False |
|---|
| 386 | print >>printFile,'\n Refinement results for histogram: v'+histogram |
|---|
| 387 | print >>printFile,135*'-' |
|---|
| 388 | try: |
|---|
| 389 | IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict, |
|---|
| 390 | rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg) |
|---|
| 391 | |
|---|
| 392 | print ' wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%( |
|---|
| 393 | Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']) |
|---|
| 394 | # add the uncertainties into the esd dictionary (sigDict) |
|---|
| 395 | sigDict = dict(zip(varyList,sig)) |
|---|
| 396 | # the uncertainties for dependent constrained parms into the esd dict |
|---|
| 397 | sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict)) |
|---|
| 398 | |
|---|
| 399 | # a dict with values & esds for dependent (constrained) parameters |
|---|
| 400 | depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart |
|---|
| 401 | if i not in varyList} |
|---|
| 402 | newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList)) |
|---|
| 403 | newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList)) |
|---|
| 404 | histRefData = { |
|---|
| 405 | 'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals, |
|---|
| 406 | 'varyListStart':varyListStart, |
|---|
| 407 | 'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict, |
|---|
| 408 | 'newCellDict':newCellDict,'depParmDict':depParmDict, |
|---|
| 409 | 'constraintInfo':constraintInfo, |
|---|
| 410 | 'parmDict':parmDict} |
|---|
| 411 | SeqResult[histogram] = histRefData |
|---|
| 412 | G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True) |
|---|
| 413 | # G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile) |
|---|
| 414 | G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile) |
|---|
| 415 | G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile) |
|---|
| 416 | G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack) |
|---|
| 417 | makeBack = False |
|---|
| 418 | NewparmDict = {} |
|---|
| 419 | # make dict of varied parameters in current histogram, renamed to |
|---|
| 420 | # next histogram, for use in next refinement. |
|---|
| 421 | if Controls['Copy2Next'] and ihst < len(histNames)-1: |
|---|
| 422 | hId = Histo[histogram]['hId'] # current histogram |
|---|
| 423 | nexthId = Histograms[histNames[ihst+1]]['hId'] |
|---|
| 424 | for parm in set(list(varyList)+list(varyListStart)): |
|---|
| 425 | items = parm.split(':') |
|---|
| 426 | if len(items) < 3: continue |
|---|
| 427 | if str(hId) in items[1]: |
|---|
| 428 | items[1] = str(nexthId) |
|---|
| 429 | newparm = ':'.join(items) |
|---|
| 430 | NewparmDict[newparm] = parmDict[parm] |
|---|
| 431 | except G2obj.G2Exception,Msg: |
|---|
| 432 | printFile.close() |
|---|
| 433 | print ' ***** Refinement aborted *****' |
|---|
| 434 | return False,Msg.msg |
|---|
| 435 | G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult) |
|---|
| 436 | printFile.close() |
|---|
| 437 | print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst' |
|---|
| 438 | print ' ***** Sequential refinement successful *****' |
|---|
| 439 | return True,'Success' |
|---|
| 440 | |
|---|
| 441 | def RetDistAngle(DisAglCtls,DisAglData): |
|---|
| 442 | '''Compute and return distances and angles |
|---|
| 443 | |
|---|
| 444 | :param dict DisAglCtls: contains distance/angle radii usually defined using |
|---|
| 445 | :func:`GSASIIgrid.DisAglDialog` |
|---|
| 446 | :param dict DisAglData: contains phase data: |
|---|
| 447 | Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used |
|---|
| 448 | for distance/angle origins and atoms to be used as targets. |
|---|
| 449 | Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`) |
|---|
| 450 | |
|---|
| 451 | :returns: AtomLabels,DistArray,AngArray where: |
|---|
| 452 | |
|---|
| 453 | **AtomLabels** is a dict of atom labels, keys are the atom number |
|---|
| 454 | |
|---|
| 455 | **DistArray** is a dict keyed by the origin atom number where the value is a list |
|---|
| 456 | of distance entries. The value for each distance is a list containing: |
|---|
| 457 | |
|---|
| 458 | 0) the target atom number (int); |
|---|
| 459 | 1) the unit cell offsets added to x,y & z (tuple of int values) |
|---|
| 460 | 2) the symmetry operator number (which may be modified to indicate centering and center of symmetry) |
|---|
| 461 | 3) an interatomic distance in A (float) |
|---|
| 462 | 4) an uncertainty on the distance in A or 0.0 (float) |
|---|
| 463 | |
|---|
| 464 | **AngArray** is a dict keyed by the origin (central) atom number where |
|---|
| 465 | the value is a list of |
|---|
| 466 | angle entries. The value for each angle entry consists of three values: |
|---|
| 467 | |
|---|
| 468 | 0) a distance item reference for one neighbor (int) |
|---|
| 469 | 1) a distance item reference for a second neighbor (int) |
|---|
| 470 | 2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats) |
|---|
| 471 | |
|---|
| 472 | The AngArray distance reference items refer directly to the index of the items in the |
|---|
| 473 | DistArray item for the list of distances for the central atom. |
|---|
| 474 | ''' |
|---|
| 475 | import numpy.ma as ma |
|---|
| 476 | |
|---|
| 477 | SGData = DisAglData['SGData'] |
|---|
| 478 | Cell = DisAglData['Cell'] |
|---|
| 479 | |
|---|
| 480 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
|---|
| 481 | covData = {} |
|---|
| 482 | if 'covData' in DisAglData: |
|---|
| 483 | covData = DisAglData['covData'] |
|---|
| 484 | covMatrix = covData['covMatrix'] |
|---|
| 485 | varyList = covData['varyList'] |
|---|
| 486 | pfx = str(DisAglData['pId'])+'::' |
|---|
| 487 | A = G2lat.cell2A(Cell[:6]) |
|---|
| 488 | cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData) |
|---|
| 489 | names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = '] |
|---|
| 490 | valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)] |
|---|
| 491 | |
|---|
| 492 | Factor = DisAglCtls['Factors'] |
|---|
| 493 | Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii']))) |
|---|
| 494 | indices = (-1,0,1) |
|---|
| 495 | Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) |
|---|
| 496 | origAtoms = DisAglData['OrigAtoms'] |
|---|
| 497 | targAtoms = DisAglData['TargAtoms'] |
|---|
| 498 | AtomLabels = {} |
|---|
| 499 | for Oatom in origAtoms: |
|---|
| 500 | AtomLabels[Oatom[0]] = Oatom[1] |
|---|
| 501 | for Oatom in targAtoms: |
|---|
| 502 | AtomLabels[Oatom[0]] = Oatom[1] |
|---|
| 503 | DistArray = {} |
|---|
| 504 | AngArray = {} |
|---|
| 505 | for Oatom in origAtoms: |
|---|
| 506 | DistArray[Oatom[0]] = [] |
|---|
| 507 | AngArray[Oatom[0]] = [] |
|---|
| 508 | OxyzNames = '' |
|---|
| 509 | IndBlist = [] |
|---|
| 510 | Dist = [] |
|---|
| 511 | Vect = [] |
|---|
| 512 | VectA = [] |
|---|
| 513 | angles = [] |
|---|
| 514 | for Tatom in targAtoms: |
|---|
| 515 | Xvcov = [] |
|---|
| 516 | TxyzNames = '' |
|---|
| 517 | if 'covData' in DisAglData: |
|---|
| 518 | OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])] |
|---|
| 519 | TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])] |
|---|
| 520 | Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix) |
|---|
| 521 | result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False) |
|---|
| 522 | BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0] |
|---|
| 523 | AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1] |
|---|
| 524 | for Txyz,Top,Tunit in result: |
|---|
| 525 | Dx = (Txyz-np.array(Oatom[3:6]))+Units |
|---|
| 526 | dx = np.inner(Amat,Dx) |
|---|
| 527 | dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5) |
|---|
| 528 | IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.)) |
|---|
| 529 | if np.any(IndB): |
|---|
| 530 | for indb in IndB: |
|---|
| 531 | for i in range(len(indb)): |
|---|
| 532 | if str(dx.T[indb][i]) not in IndBlist: |
|---|
| 533 | IndBlist.append(str(dx.T[indb][i])) |
|---|
| 534 | unit = Units[indb][i] |
|---|
| 535 | tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2]) |
|---|
| 536 | pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData) |
|---|
| 537 | sig = 0.0 |
|---|
| 538 | if len(Xvcov): |
|---|
| 539 | sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx))) |
|---|
| 540 | Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig]) |
|---|
| 541 | if (Dist[-1][-2]-AsumR) <= 0.: |
|---|
| 542 | Vect.append(dx.T[indb][i]/Dist[-1][-2]) |
|---|
| 543 | VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top]) |
|---|
| 544 | else: |
|---|
| 545 | Vect.append([0.,0.,0.]) |
|---|
| 546 | VectA.append([]) |
|---|
| 547 | for D in Dist: |
|---|
| 548 | DistArray[Oatom[0]].append(D[1:]) |
|---|
| 549 | Vect = np.array(Vect) |
|---|
| 550 | angles = np.zeros((len(Vect),len(Vect))) |
|---|
| 551 | angsig = np.zeros((len(Vect),len(Vect))) |
|---|
| 552 | for i,veca in enumerate(Vect): |
|---|
| 553 | if np.any(veca): |
|---|
| 554 | for j,vecb in enumerate(Vect): |
|---|
| 555 | if np.any(vecb): |
|---|
| 556 | angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData) |
|---|
| 557 | if i <= j: continue |
|---|
| 558 | AngArray[Oatom[0]].append((i,j, |
|---|
| 559 | G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData))) |
|---|
| 560 | return AtomLabels,DistArray,AngArray |
|---|
| 561 | |
|---|
| 562 | def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout): |
|---|
| 563 | '''Print distances and angles |
|---|
| 564 | |
|---|
| 565 | :param dict DisAglCtls: contains distance/angle radii usually defined using |
|---|
| 566 | :func:`GSASIIgrid.DisAglDialog` |
|---|
| 567 | :param dict DisAglData: contains phase data: |
|---|
| 568 | Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used |
|---|
| 569 | for distance/angle origins and atoms to be used as targets. |
|---|
| 570 | Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`) |
|---|
| 571 | :param file out: file object for output. Defaults to sys.stdout. |
|---|
| 572 | ''' |
|---|
| 573 | import numpy.ma as ma |
|---|
| 574 | def MyPrint(s): |
|---|
| 575 | out.write(s+'\n') |
|---|
| 576 | # print(s,file=out) # use in Python 3 |
|---|
| 577 | |
|---|
| 578 | def ShowBanner(name): |
|---|
| 579 | MyPrint(80*'*') |
|---|
| 580 | MyPrint(' Interatomic Distances and Angles for phase '+name) |
|---|
| 581 | MyPrint((80*'*')+'\n') |
|---|
| 582 | |
|---|
| 583 | ShowBanner(DisAglCtls['Name']) |
|---|
| 584 | SGData = DisAglData['SGData'] |
|---|
| 585 | SGtext,SGtable = G2spc.SGPrint(SGData) |
|---|
| 586 | for line in SGtext: MyPrint(line) |
|---|
| 587 | if len(SGtable): |
|---|
| 588 | for i,item in enumerate(SGtable[::2]): |
|---|
| 589 | line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30)) |
|---|
| 590 | MyPrint(line) |
|---|
| 591 | else: |
|---|
| 592 | MyPrint(' ( 1) %s'%(SGtable[0])) |
|---|
| 593 | Cell = DisAglData['Cell'] |
|---|
| 594 | |
|---|
| 595 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
|---|
| 596 | covData = {} |
|---|
| 597 | if 'covData' in DisAglData: |
|---|
| 598 | covData = DisAglData['covData'] |
|---|
| 599 | covMatrix = covData['covMatrix'] |
|---|
| 600 | varyList = covData['varyList'] |
|---|
| 601 | pfx = str(DisAglData['pId'])+'::' |
|---|
| 602 | A = G2lat.cell2A(Cell[:6]) |
|---|
| 603 | cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData) |
|---|
| 604 | names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = '] |
|---|
| 605 | valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)] |
|---|
| 606 | line = '\n Unit cell:' |
|---|
| 607 | for name,vals in zip(names,valEsd): |
|---|
| 608 | line += name+vals |
|---|
| 609 | MyPrint(line) |
|---|
| 610 | else: |
|---|
| 611 | MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+ |
|---|
| 612 | ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+ |
|---|
| 613 | ('%.3f'%Cell[5])+' volume = '+('%.3f'%Cell[6])) |
|---|
| 614 | |
|---|
| 615 | AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData) |
|---|
| 616 | origAtoms = DisAglData['OrigAtoms'] |
|---|
| 617 | targAtoms = DisAglData['TargAtoms'] |
|---|
| 618 | for Oatom in origAtoms: |
|---|
| 619 | i = Oatom[0] |
|---|
| 620 | Dist = DistArray[i] |
|---|
| 621 | nDist = len(Dist) |
|---|
| 622 | angles = np.zeros((nDist,nDist)) |
|---|
| 623 | angsig = np.zeros((nDist,nDist)) |
|---|
| 624 | for k,j,tup in AngArray[i]: |
|---|
| 625 | angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup |
|---|
| 626 | line = '' |
|---|
| 627 | for i,x in enumerate(Oatom[3:6]): |
|---|
| 628 | line += ('%12.5f'%x).rstrip('0') |
|---|
| 629 | MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip()) |
|---|
| 630 | MyPrint(80*'*') |
|---|
| 631 | line = '' |
|---|
| 632 | for dist in Dist[:-1]: |
|---|
| 633 | line += '%12s'%(AtomLabels[dist[0]].center(12)) |
|---|
| 634 | MyPrint(' To cell +(sym. op.) dist. '+line.rstrip()) |
|---|
| 635 | for i,dist in enumerate(Dist): |
|---|
| 636 | line = '' |
|---|
| 637 | for j,angle in enumerate(angles[i][0:i]): |
|---|
| 638 | sig = angsig[i][j] |
|---|
| 639 | if angle: |
|---|
| 640 | if sig: |
|---|
| 641 | line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12)) |
|---|
| 642 | else: |
|---|
| 643 | val = '%.3f'%(angle) |
|---|
| 644 | line += '%12s'%(val.center(12)) |
|---|
| 645 | else: |
|---|
| 646 | line += 12*' ' |
|---|
| 647 | if dist[4]: #sig exists! |
|---|
| 648 | val = G2mth.ValEsd(dist[3],dist[4]) |
|---|
| 649 | else: |
|---|
| 650 | val = '%8.4f'%(dist[3]) |
|---|
| 651 | tunit = '[%2d%2d%2d]'% dist[1] |
|---|
| 652 | MyPrint((' %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip()) |
|---|
| 653 | |
|---|
| 654 | def DisAglTor(DATData): |
|---|
| 655 | 'Needs a doc string' |
|---|
| 656 | SGData = DATData['SGData'] |
|---|
| 657 | Cell = DATData['Cell'] |
|---|
| 658 | |
|---|
| 659 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
|---|
| 660 | covData = {} |
|---|
| 661 | pfx = '' |
|---|
| 662 | if 'covData' in DATData: |
|---|
| 663 | covData = DATData['covData'] |
|---|
| 664 | covMatrix = covData['covMatrix'] |
|---|
| 665 | varyList = covData['varyList'] |
|---|
| 666 | pfx = str(DATData['pId'])+'::' |
|---|
| 667 | Datoms = [] |
|---|
| 668 | Oatoms = [] |
|---|
| 669 | for i,atom in enumerate(DATData['Datoms']): |
|---|
| 670 | symop = atom[-1].split('+') |
|---|
| 671 | if len(symop) == 1: |
|---|
| 672 | symop.append('0,0,0') |
|---|
| 673 | symop[0] = int(symop[0]) |
|---|
| 674 | symop[1] = eval(symop[1]) |
|---|
| 675 | atom.append(symop) |
|---|
| 676 | Datoms.append(atom) |
|---|
| 677 | oatom = DATData['Oatoms'][i] |
|---|
| 678 | names = ['','',''] |
|---|
| 679 | if pfx: |
|---|
| 680 | names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])] |
|---|
| 681 | oatom += [names,] |
|---|
| 682 | Oatoms.append(oatom) |
|---|
| 683 | atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms] |
|---|
| 684 | if DATData['Natoms'] == 4: #torsion |
|---|
| 685 | Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData) |
|---|
| 686 | print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig) |
|---|
| 687 | print ' NB: Atom sequence determined by selection order' |
|---|
| 688 | return # done with torsion |
|---|
| 689 | elif DATData['Natoms'] == 3: #angle |
|---|
| 690 | Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData) |
|---|
| 691 | print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig) |
|---|
| 692 | print ' NB: Atom sequence determined by selection order' |
|---|
| 693 | else: #2 atoms - distance |
|---|
| 694 | Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData) |
|---|
| 695 | print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig) |
|---|
| 696 | |
|---|
| 697 | def BestPlane(PlaneData): |
|---|
| 698 | 'Needs a doc string' |
|---|
| 699 | |
|---|
| 700 | def ShowBanner(name): |
|---|
| 701 | print 80*'*' |
|---|
| 702 | print ' Best plane result for phase '+name |
|---|
| 703 | print 80*'*','\n' |
|---|
| 704 | |
|---|
| 705 | ShowBanner(PlaneData['Name']) |
|---|
| 706 | |
|---|
| 707 | Cell = PlaneData['Cell'] |
|---|
| 708 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
|---|
| 709 | Atoms = PlaneData['Atoms'] |
|---|
| 710 | sumXYZ = np.zeros(3) |
|---|
| 711 | XYZ = [] |
|---|
| 712 | Natoms = len(Atoms) |
|---|
| 713 | for atom in Atoms: |
|---|
| 714 | xyz = np.array(atom[3:6]) |
|---|
| 715 | XYZ.append(xyz) |
|---|
| 716 | sumXYZ += xyz |
|---|
| 717 | sumXYZ /= Natoms |
|---|
| 718 | XYZ = np.array(XYZ)-sumXYZ |
|---|
| 719 | XYZ = np.inner(Amat,XYZ).T |
|---|
| 720 | Zmat = np.zeros((3,3)) |
|---|
| 721 | for i,xyz in enumerate(XYZ): |
|---|
| 722 | Zmat += np.outer(xyz.T,xyz) |
|---|
| 723 | print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]) |
|---|
| 724 | Evec,Emat = nl.eig(Zmat) |
|---|
| 725 | Evec = np.sqrt(Evec)/(Natoms-3) |
|---|
| 726 | Order = np.argsort(Evec) |
|---|
| 727 | XYZ = np.inner(XYZ,Emat.T).T |
|---|
| 728 | XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T |
|---|
| 729 | print ' Atoms in Cartesian best plane coordinates:' |
|---|
| 730 | print ' Name X Y Z' |
|---|
| 731 | for i,xyz in enumerate(XYZ): |
|---|
| 732 | print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]) |
|---|
| 733 | print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]) |
|---|
| 734 | |
|---|
| 735 | def main(): |
|---|
| 736 | 'Needs a doc string' |
|---|
| 737 | arg = sys.argv |
|---|
| 738 | if len(arg) > 1: |
|---|
| 739 | GPXfile = arg[1] |
|---|
| 740 | if not ospath.exists(GPXfile): |
|---|
| 741 | print 'ERROR - ',GPXfile," doesn't exist!" |
|---|
| 742 | exit() |
|---|
| 743 | GPXpath = ospath.dirname(arg[1]) |
|---|
| 744 | Refine(GPXfile,None) |
|---|
| 745 | else: |
|---|
| 746 | print 'ERROR - missing filename' |
|---|
| 747 | exit() |
|---|
| 748 | print "Done" |
|---|
| 749 | |
|---|
| 750 | if __name__ == '__main__': |
|---|
| 751 | main() |
|---|
