Changeset 960
- Timestamp:
- Jun 19, 2013 5:51:24 PM (10 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIIO.py
r956 r960 1608 1608 for item in consDict: 1609 1609 constList += consDict[item] 1610 constDict,fixedList,ignored = G2stIO.ProcessConstraints(constList)1611 1610 # now process the constraints 1612 1611 G2mv.InitVars() 1613 varyList = covDict.get('varyList',[]) 1612 constDict,fixedList,ignored = G2stIO.ProcessConstraints(constList) 1613 varyList = covDict.get('varyListStart') 1614 if varyList is None and len(constDict) == 0: 1615 # no constraints can use varyList 1616 varyList = covDict.get('varyList') 1617 elif varyList is None: 1618 # old GPX file from before pre-constraint varyList is saved 1619 print ' *** Old refinement: Please use Calculate/Refine to redo ***' 1620 raise Exception(' *** CIF creation aborted ***') 1621 else: 1622 varyList = list(varyList) 1614 1623 try: 1615 1624 groups,parmlist = G2mv.GroupConstraints(constDict) … … 1626 1635 # and add their uncertainties into the esd dictionary (sigDict) 1627 1636 if covDict.get('covMatrix') is not None: 1628 self.sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'], varyList,self.parmDict))1637 self.sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],covDict['varyList'],self.parmDict)) 1629 1638 1630 1639 def loadTree(self): -
trunk/GSASIImath.py
r957 r960 1292 1292 out += ("({:d})").format(intesd) # add the esd 1293 1293 elif nTZ and '.' in out: 1294 out = out.rstrip('0') # strip digits to right of decimal 1294 out = out.rstrip('0') # strip zeros to right of decimal 1295 out = out.rstrip('.') # and decimal place when not needed 1295 1296 if valoff != 0: 1296 1297 out += ("e{:d}").format(valoff) # add an exponent, when needed -
trunk/GSASIIstrIO.py
r953 r960 83 83 84 84 def ProcessConstraints(constList): 85 "interpret constraints" 85 """Interpret the constraints in the constList input into a dictionary, etc. 86 87 :param list constList: a list of lists where each item in the outer list 88 specifies a constraint of some form. The last item in each inner list 89 determines which of the four constraints types has been input: 90 91 * h (hold): a single variable that will not be varied. It 92 will be removed from the varyList later. 93 * c (constraint): specifies a linear relationship that 94 can be varied as a new grouped variable 95 a fixed value. 96 * f (fixed): specifies a linear relationship that is assigned 97 a fixed value. 98 * e (equivalence): specifies a series of variables where the 99 first variable in the last can be used to generate the 100 values for all the remaining variables. 101 102 :returns: a tuple of (constDict,fixedList,ignored) where: 103 104 * constDict (list) contains the constraint relationships 105 * fixedList (list) contains the fixed values for type 106 of constraint. 107 * ignored (int) counts the number of invalid constraint items 108 (should always be zero!) 109 110 """ 86 111 constDict = [] 87 112 fixedList = [] -
trunk/GSASIIstrMain.py
r956 r960 88 88 G2stIO.GetFprime(calcControls,Histograms) 89 89 # do constraint processing 90 varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed 90 91 try: 91 92 groups,parmlist = G2mv.GroupConstraints(constrDict) … … 176 177 newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList) 177 178 covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals, 178 'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict} 179 'varyListStart':varyListStart, 180 'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict, 181 'newCellDict':newCellDict} 179 182 # add the uncertainties into the esd dictionary (sigDict) 180 183 sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict)) … … 289 292 G2stIO.GetFprime(calcControls,Histo) 290 293 # do constraint processing 294 varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed 291 295 try: 292 296 groups,parmlist = G2mv.GroupConstraints(constrDict) … … 366 370 newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList)) 367 371 covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals, 368 'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict} 372 'varyListStart':varyListStart, 373 'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict, 374 'newCellDict':newCellDict} 369 375 # add the uncertainties into the esd dictionary (sigDict) 370 376 G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True) -
trunk/exports/G2cif.py
r956 r960 80 80 WriteCIFitem('_refine_ls_number_parameters',vars) 81 81 try: 82 GOF = str(self.OverallParms['Covariance']['Rvals']['GOF'])82 GOF = G2mth.ValEsd(self.OverallParms['Covariance']['Rvals']['GOF'],-0.009) 83 83 except: 84 84 GOF = '?' 85 85 WriteCIFitem('_refine_ls_goodness_of_fit_all',GOF) 86 #WriteCIFitem('_refine_ls_number_restraints',TEXT(1:7)) 86 87 #WriteCIFitem('_refine_ls_number_restraints',TEXT) 88 87 89 # other things to consider reporting 88 90 # _refine_ls_number_reflns … … 131 133 print getCallerDocString() 132 134 133 def FormatSH(phase dict):135 def FormatSH(phasenam): 134 136 'Format a full spherical harmonics texture description as a string' 135 # SH Texture137 phasedict = self.Phases[phasenam] # pointer to current phase info 136 138 pfx = str(phasedict['pId'])+'::' 137 139 s = "" … … 160 162 return s 161 163 162 def FormatHAPpo(phasedict): 163 'return the March-Dollase/SH correction for every histogram in the current phase' 164 def FormatHAPpo(phasenam): 165 '''return the March-Dollase/SH correction for every 166 histogram in the current phase formatted into a 167 character string 168 ''' 169 phasedict = self.Phases[phasenam] # pointer to current phase info 164 170 s = '' 165 171 for histogram in sorted(phasedict['Histograms']): 172 if histogram.startswith("HKLF"): continue # powder only 166 173 Histogram = self.Histograms.get(histogram) 167 174 if not Histogram: continue … … 197 204 s += s1 198 205 return s 206 207 def FormatInstProfile(instparmdict): 208 '''Format the instrumental profile parameters with a 209 string description. Will only be called on PWDR histograms 210 ''' 211 #print instparmdict[0].keys() 212 return 'TODO: Instrument profile goes here' 213 214 def FormatPhaseProfile(phasenam): 215 '''Format the phase-related profile parameters (size/strain) 216 with a string description. 217 return an empty string or None there are no 218 powder histograms for this phase. 219 ''' 220 s = '' 221 phasedict = self.Phases[phasenam] # pointer to current phase info 222 for histogram in sorted(phasedict['Histograms']): 223 if histogram.startswith("HKLF"): continue # powder only 224 Histogram = self.Histograms.get(histogram) 225 if not Histogram: continue 226 hapData = phasedict['Histograms'][histogram] 227 return 'TODO: Phase profile goes here' 199 228 200 229 def FmtAtomType(sym): 201 230 'Reformat a GSAS-II atom type symbol to match CIF rules' 202 231 sym = sym.replace('_','') # underscores are not allowed: no isotope designation? 203 # in CIF oxidation statesymbols come after, not before232 # in CIF, oxidation state sign symbols come after, not before 204 233 if '+' in sym: 205 234 sym = sym.replace('+','') + '+' … … 241 270 242 271 def WriteAtomsNuclear(phasenam): 243 phasedict = self.Phases[phasenam] # pointer to current phase info 272 'Write atom positions to CIF' 273 phasedict = self.Phases[phasenam] # pointer to current phase info 244 274 General = phasedict['General'] 275 cx,ct,cs,cia = General['AtomPtrs'] 276 Atoms = phasedict['Atoms'] 277 cfrac = cx+3 278 fpfx = str(phasedict['pId'])+'::Afrac:' 279 for i,at in enumerate(Atoms): 280 fval = self.parmDict.get(fpfx+str(i),at[cfrac]) 281 if fval != 0.0: 282 break 283 else: 284 WriteCIFitem('\n# PHASE HAS NO ATOMS!') 285 return 286 245 287 WriteCIFitem('\n# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS') 246 Atoms = phasedict['Atoms']247 288 WriteCIFitem('loop_ '+ 248 289 '\n\t_atom_site_label'+ … … 256 297 '\n\t_atom_site_symmetry_multiplicity') 257 298 258 varnames = { 3:'Ax',4:'Ay',5:'Az',6:'Afrac',259 10:'AUiso',11:'AU11',12:'AU22',13:'AU33',260 14:'AU12',15:'AU13',16:'AU23'}299 varnames = {cx:'Ax',cx+1:'Ay',cx+2:'Az',cx+3:'Afrac', 300 cia+1:'AUiso',cia+2:'AU11',cia+3:'AU22',cia+4:'AU33', 301 cia+5:'AU12',cia+6:'AU13',cia+7:'AU23'} 261 302 labellist = [] 262 303 263 304 pfx = str(phasedict['pId'])+'::' 305 # loop over all atoms 306 naniso = 0 264 307 for i,at in enumerate(Atoms): 265 print at 266 for i,at in enumerate(Atoms): 267 s = " " 268 s += PutInCol(MakeUniqueLabel(at[0],labellist),6) # label 269 #s += PutInCol(MakeUniqueLabel('A',labellist),6) # label 270 s += PutInCol(FmtAtomType(at[1]),6) # type 271 if at[9] == 'I': 308 s = PutInCol(MakeUniqueLabel(at[ct-1],labellist),6) # label 309 fval = self.parmDict.get(fpfx+str(i),at[cfrac]) 310 if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact) 311 s += PutInCol(FmtAtomType(at[ct]),4) # type 312 if at[cia] == 'I': 272 313 adp = 'Uiso ' 273 314 else: 274 315 adp = 'Uani ' 316 naniso += 1 317 # compute Uequiv crudely 318 # correct: Defined as "1/3 trace of diagonalized U matrix". 319 # SEE cell2GS & Uij2Ueqv to GSASIIlattice. Former is needed to make the GS matrix used by the latter. 275 320 t = 0.0 276 for j in ( 11,12,13):277 var = pfx+varnames[ j]+":"+str(i)278 t += self.parmDict.get(var,at[ j])279 for j in ( 3,4,5,6,10):280 if j in ( 3,4,5):321 for j in (2,3,4): 322 var = pfx+varnames[cia+j]+":"+str(i) 323 t += self.parmDict.get(var,at[cia+j]) 324 for j in (cx,cx+1,cx+2,cx+3,cia+1): 325 if j in (cx,cx+1,cx+2): 281 326 dig = 11 282 327 sigdig = -0.00009 283 328 else: 284 dig = 9329 dig = 10 285 330 sigdig = -0.009 286 331 var = pfx+varnames[j]+":"+str(i) … … 288 333 if dvar not in self.sigDict: 289 334 dvar = var 290 if j == 10 and adp == 'Uani ': 291 # compute Uequiv crudely 292 t = 0.0 293 for k in (11,12,13): 294 var = pfx+varnames[j]+":"+str(i) 295 t += self.parmDict.get(var,at[k]) 335 if j == cia+1 and adp == 'Uani ': 296 336 val = t/3. 297 337 sig = sigdig 298 338 else: 339 #print var,(var in self.parmDict),(var in self.sigDict) 299 340 val = self.parmDict.get(var,at[j]) 300 341 sig = self.sigDict.get(dvar,sigdig) 301 342 s += PutInCol(G2mth.ValEsd(val,sig),dig) 302 343 s += adp 303 print s 304 344 s += PutInCol(at[cs+1],3) 345 WriteCIFitem(s) 346 if naniso == 0: return 347 # now loop over aniso atoms 348 WriteCIFitem('\nloop_' + '\n\t_atom_site_aniso_label' + 349 '\n\t_atom_site_aniso_U_11' + '\n\t_atom_site_aniso_U_12' + 350 '\n\t_atom_site_aniso_U_13' + '\n\t_atom_site_aniso_U_22' + 351 '\n\t_atom_site_aniso_U_23' + '\n\t_atom_site_aniso_U_33') 352 for i,at in enumerate(Atoms): 353 fval = self.parmDict.get(fpfx+str(i),at[cfrac]) 354 if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact) 355 if at[cia] == 'I': continue 356 s = PutInCol(labellist[i],6) # label 357 for j in (2,3,4,5,6,7): 358 sigdig = -0.0009 359 var = pfx+varnames[cia+j]+":"+str(i) 360 val = self.parmDict.get(var,at[cia+j]) 361 sig = self.sigDict.get(var,sigdig) 362 s += PutInCol(G2mth.ValEsd(val,sig),11) 363 WriteCIFitem(s) 364 365 def HillSortElements(elmlist): 366 '''Sort elements in "Hill" order: C, H, others, (where others 367 are alphabetical). 368 369 :params list elmlist: a list of element strings 370 371 :returns: a sorted list of element strings 372 ''' 373 newlist = [] 374 oldlist = elmlist[:] 375 for elm in ('C','H'): 376 if elm in elmlist: 377 newlist.append(elm) 378 oldlist.pop(oldlist.index(elm)) 379 return newlist+sorted(oldlist) 380 381 def WriteComposition(phasenam): 382 '''determine the composition for the unit cell, crudely determine Z and 383 then compute the composition in formula units 384 ''' 385 phasedict = self.Phases[phasenam] # pointer to current phase info 386 General = phasedict['General'] 387 Z = General.get('cellZ',0.0) 388 cx,ct,cs,cia = General['AtomPtrs'] 389 Atoms = phasedict['Atoms'] 390 fpfx = str(phasedict['pId'])+'::Afrac:' 391 cfrac = cx+3 392 cmult = cs+1 393 compDict = {} # combines H,D & T 394 sitemultlist = [] 395 massDict = dict(zip(General['AtomTypes'],General['AtomMass'])) 396 cellmass = 0 397 for i,at in enumerate(Atoms): 398 atype = at[ct].strip() 399 if atype.find('-') != -1: atype = atype.split('-')[0] 400 if atype.find('+') != -1: atype = atype.split('+')[0] 401 atype = atype[0].upper()+atype[1:2].lower() # force case conversion 402 if atype == "D" or atype == "D": atype = "H" 403 fvar = fpfx+str(i) 404 fval = self.parmDict.get(fvar,at[cfrac]) 405 mult = at[cmult] 406 if not massDict.get(at[ct]): 407 print 'No mass found for atom type '+at[ct] 408 print 'Will not compute cell contents for phase '+phasenam 409 return 410 cellmass += massDict[at[ct]]*mult*fval 411 compDict[atype] = compDict.get(atype,0.0) + mult*fval 412 if fval == 1: sitemultlist.append(mult) 413 if len(compDict.keys()) == 0: return # no elements! 414 if Z < 1: # Z has not been computed or set by user 415 Z = 1 416 for i in range(2,min(sitemultlist)+1): 417 for m in sitemultlist: 418 if m % i != 0: 419 break 420 else: 421 Z = i 422 General['cellZ'] = Z # save it 423 424 # when scattering factors are included in the CIF, this needs to be 425 # added to the loop here but only in the one-block case. 426 # For multiblock CIFs, scattering factors go in the histogram 427 # blocks (for all atoms in all appropriate phases) 428 429 #if oneblock: # add scattering factors for current phase here 430 WriteCIFitem('\nloop_ _atom_type_symbol _atom_type_number_in_cell') 431 formula = '' 432 reload(G2mth) 433 for elem in HillSortElements(compDict.keys()): 434 WriteCIFitem(' ' + PutInCol(elem,4) + 435 G2mth.ValEsd(compDict[elem],-0.009,True)) 436 if formula: formula += " " 437 formula += elem 438 if compDict[elem] == Z: continue 439 formula += G2mth.ValEsd(compDict[elem]/Z,-0.009,True) 440 WriteCIFitem( '\n# Note that Z affects _cell_formula_sum and _weight') 441 WriteCIFitem( '_cell_formula_units_Z',str(Z)) 442 WriteCIFitem( '_chemical_formula_sum',formula) 443 WriteCIFitem( '_chemical_formula_weight', 444 G2mth.ValEsd(cellmass/Z,-0.09,True)) 445 305 446 def WritePhaseInfo(phasenam): 306 447 # see writepha.for … … 309 450 phasedict = self.Phases[phasenam] # pointer to current phase info 310 451 WriteCIFitem('_pd_phase_name', phasenam) 311 # replace some of these later312 #General = phasedict['General']313 #SGData = phasedict['General']['SGData']314 #cell = General['Cell']315 452 pfx = str(phasedict['pId'])+'::' 316 #covData = self.OverallParms['Covariance']317 453 A,sigA = G2stIO.cellFill(pfx,phasedict['General']['SGData'],self.parmDict,self.sigDict) 318 454 cellSig = G2stIO.getCellEsd(pfx, … … 346 482 WriteCIFitem(' {:3d} {:}'.format(i,op.lower())) 347 483 348 # preferred orientation (always reported by phase)349 SH = FormatSH(phasedict)350 MD = FormatHAPpo(phasedict)351 if SH and MD:352 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD)353 elif SH or MD:354 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD)355 else:356 WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none')357 358 484 # loop over histogram(s) used in this phase 359 if oneblock: 360 # report all profile information here (include histogram & phase) 361 # _pd_proc_ls_profile_function 362 pass 363 else: 364 # pointers to histograms used in this phase 485 if not oneblock and not self.quickmode: 486 # report pointers to the histograms used in this phase 365 487 histlist = [] 366 488 for hist in self.Phases[phasenam]['Histograms']: … … 383 505 WriteCIFitem('loop_ _pd_block_diffractogram_id') 384 506 385 # report sample profile information here (phase only) 386 # _pd_proc_ls_profile_function 387 507 # report atom params 388 508 if phasedict['General']['Type'] == 'nuclear': #this needs macromolecular variant, etc! 389 509 WriteAtomsNuclear(phasenam) 390 510 else: 391 511 raise Exception,"no export for mm coordinates implemented" 392 393 394 raise Exception,'Testing' 395 396 WriteCIFitem('loop_' + '\n\t_atom_site_aniso_label' + 397 '\n\t_atom_site_aniso_U_11' + '\n\t_atom_site_aniso_U_12' + 398 '\n\t_atom_site_aniso_U_13' + '\n\t_atom_site_aniso_U_22' + 399 '\n\t_atom_site_aniso_U_23' + '\n\t_atom_site_aniso_U_33') 400 401 # process the chemical formula: pick a Z value & generate molecular weight 402 # find the maximum possible Z value 403 404 # order the elements in "Hill" order: C, H, D, T & alphabetical or alphabetical 405 if not oneblock: # in single block, this is combined with the scattering factors 406 WriteCIFitem('loop_ _atom_type_symbol _atom_type_number_in_cell') 407 408 WriteCIFitem('# If you change Z, be sure to change all 3 of the following') 409 WriteCIFitem( '_chemical_formula_sum',text) 410 #WRITE(TEXT,'(F15.2)') ATMASS 411 WriteCIFitem( '_chemical_formula_weight',text) 412 #WRITE(TEXT,'(I4)') Z 413 WriteCIFitem( '_cell_formula_units_Z',text) 512 # report cell contents 513 WriteComposition(phasenam) 514 # if not self.quickmode: 515 # report distances and angles 516 # WriteDistances(phasenam) 517 518 #raise Exception,'Testing' 414 519 415 520 #C now loop over interatomic distances for this phase … … 714 819 ): return 715 820 821 if oneblock and not self.quickmode: 822 # select a dataset to use (there should only be one set in one block, 823 # but take whatever comes 1st 824 for hist in self.Histograms: 825 histblk = self.Histograms[hist] 826 if hist.startswith("PWDR"): 827 instnam = histblk["Sample Parameters"]['InstrName'] 828 break # ignore all but 1st data histogram 829 elif hist.startswith("HKLF"): 830 instnam = histblk["Instrument Parameters"][0]['InstrName'] 831 break # ignore all but 1st data histogram 716 832 #====================================================================== 717 833 # Start writing the CIF - single block … … 724 840 phaseblk = self.Phases[phasenam] # pointer to current phase info 725 841 if not self.quickmode: 726 # select data, should only be one set in oneblock, but take whatever comes 1st727 for hist in self.Histograms:728 histblk = self.Histograms[hist]729 if hist.startswith("PWDR"):730 instnam = histblk["Sample Parameters"]['InstrName']731 break # ignore all but 1st data histogram732 elif hist.startswith("HKLF"):733 instnam = histblk["Instrument Parameters"][0]['InstrName']734 break # ignore all but 1st data histogram735 842 instnam = instnam.replace(' ','') 736 843 WriteCIFitem('_pd_block_id', … … 743 850 WritePhaseTemplate() 744 851 WritePhaseInfo(phasenam) 745 if not self.quickmode: 746 if hist.startswith("PWDR"): 747 WritePowderTemplate() 748 WritePowderData(hist) 749 elif hist.startswith("HKLF"): 750 WriteSnglXtalTemplate() 751 WriteSingleXtalData(hist) 852 if hist.startswith("PWDR") and not self.quickmode: 853 # preferred orientation 854 SH = FormatSH(phasenam) 855 MD = FormatHAPpo(phasenam) 856 if SH and MD: 857 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD) 858 elif SH or MD: 859 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD) 860 else: 861 WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none') 862 # report profile, since one-block: include both histogram and phase info 863 WriteCIFitem('_pd_proc_ls_profile_function', 864 FormatInstProfile(histblk["Instrument Parameters"]) 865 + '\n' + 866 FormatPhaseProfile(phasenam)) 867 WritePowderTemplate() 868 WritePowderData(hist) 869 elif hist.startswith("HKLF") and not self.quickmode: 870 WriteSnglXtalTemplate() 871 WriteSingleXtalData(hist) 752 872 else: 753 873 #====================================================================== … … 809 929 WritePhaseTemplate() 810 930 WritePhaseInfo(phasenam) 811 931 # preferred orientation 932 SH = FormatSH(phasenam) 933 MD = FormatHAPpo(phasenam) 934 if SH and MD: 935 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + '\n' + MD) 936 elif SH or MD: 937 WriteCIFitem('_pd_proc_ls_pref_orient_corr', SH + MD) 938 else: 939 WriteCIFitem('_pd_proc_ls_pref_orient_corr', 'none') 940 # report sample profile terms 941 PP = FormatPhaseProfile(phasenam) 942 if PP: 943 WriteCIFitem('_pd_proc_ls_profile_function',PP) 944 812 945 #============================================================ 813 946 # loop over histograms, exporting them … … 818 951 WriteCIFitem('\ndata_'+self.CIFname+"_pwd_"+str(i)) 819 952 #instnam = histblk["Sample Parameters"]['InstrName'] 953 # report instrumental profile terms 954 WriteCIFitem('_pd_proc_ls_profile_function', 955 FormatInstProfile(histblk["Instrument Parameters"])) 820 956 WriteCIFitem('# Information for histogram '+str(i)+': '+ 821 957 hist)
Note: See TracChangeset
for help on using the changeset viewer.