468 | | |
469 | | GPXpath,GPXname = ospath.split(GPXfile) |
470 | | |
471 | | if GPXpath == '': GPXpath = '.' |
472 | | |
473 | | Name = ospath.splitext(GPXname)[0] |
474 | | |
475 | | files = os.listdir(GPXpath) |
476 | | |
477 | | last = 0 |
478 | | |
479 | | for name in files: |
480 | | |
481 | | name = name.split('.') |
482 | | |
483 | | if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]: |
484 | | |
485 | | if makeBack: |
486 | | |
487 | | last = max(last,int(name[-2].strip('bak'))+1) |
488 | | |
| 259 | GPXback = getBackupName(GPXfile,makeBack) |
| 260 | dfu.copy_file(GPXfile,GPXback) |
| 261 | return GPXback |
| 262 | |
| 263 | def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True): |
| 264 | ''' Updates gpxfile from all histograms that are found in any phase |
| 265 | and any phase that used a histogram |
| 266 | input: |
| 267 | GPXfile = .gpx full file name |
| 268 | Histograms = dictionary of histograms as {name:data,...} |
| 269 | Phases = dictionary of phases that use histograms |
| 270 | CovData = dictionary of refined variables, varyList, & covariance matrix |
| 271 | makeBack = True if new backup of .gpx file is to be made; else use the last one made |
| 272 | ''' |
| 273 | |
| 274 | GPXback = GPXBackup(GPXfile,makeBack) |
| 275 | print '\n',135*'-' |
| 276 | print 'Read from file:',GPXback |
| 277 | print 'Save to file :',GPXfile |
| 278 | infile = open(GPXback,'rb') |
| 279 | outfile = open(GPXfile,'wb') |
| 280 | while True: |
| 281 | try: |
| 282 | data = cPickle.load(infile) |
| 283 | except EOFError: |
| 284 | break |
| 285 | datum = data[0] |
| 286 | # print 'read: ',datum[0] |
| 287 | if datum[0] == 'Phases': |
| 288 | for iphase in range(len(data)): |
| 289 | if data[iphase][0] in Phases: |
| 290 | phaseName = data[iphase][0] |
| 291 | data[iphase][1].update(Phases[phaseName]) |
| 292 | elif datum[0] == 'Covariance': |
| 293 | data[0][1] = CovData |
| 294 | try: |
| 295 | histogram = Histograms[datum[0]] |
| 296 | # print 'found ',datum[0] |
| 297 | data[0][1][1] = histogram['Data'] |
| 298 | for datus in data[1:]: |
| 299 | # print ' read: ',datus[0] |
| 300 | if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: |
| 301 | datus[1] = histogram[datus[0]] |
| 302 | except KeyError: |
| 303 | pass |
| 304 | |
| 305 | cPickle.dump(data,outfile,1) |
| 306 | infile.close() |
| 307 | outfile.close() |
| 308 | print 'GPX file save successful' |
| 309 | |
| 310 | def SetSeqResult(GPXfile,Histograms,SeqResult): |
| 311 | GPXback = GPXBackup(GPXfile) |
| 312 | print '\n',135*'-' |
| 313 | print 'Read from file:',GPXback |
| 314 | print 'Save to file :',GPXfile |
| 315 | infile = open(GPXback,'rb') |
| 316 | outfile = open(GPXfile,'wb') |
| 317 | while True: |
| 318 | try: |
| 319 | data = cPickle.load(infile) |
| 320 | except EOFError: |
| 321 | break |
| 322 | datum = data[0] |
| 323 | if datum[0] == 'Sequental results': |
| 324 | data[0][1] = SeqResult |
| 325 | try: |
| 326 | histogram = Histograms[datum[0]] |
| 327 | data[0][1][1] = histogram['Data'] |
| 328 | for datus in data[1:]: |
| 329 | if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: |
| 330 | datus[1] = histogram[datus[0]] |
| 331 | except KeyError: |
| 332 | pass |
| 333 | |
| 334 | cPickle.dump(data,outfile,1) |
| 335 | infile.close() |
| 336 | outfile.close() |
| 337 | print 'GPX file save successful' |
| 338 | |
| 339 | def GetPWDRdata(GPXfile,PWDRname): |
| 340 | ''' Returns powder data from GSASII gpx file |
| 341 | input: |
| 342 | GPXfile = .gpx full file name |
| 343 | PWDRname = powder histogram name as obtained from GetHistogramNames |
| 344 | return: |
| 345 | PWDRdata = powder data dictionary with: |
| 346 | Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters |
| 347 | |
| 348 | ''' |
| 349 | file = open(GPXfile,'rb') |
| 350 | PWDRdata = {} |
| 351 | while True: |
| 352 | try: |
| 353 | data = cPickle.load(file) |
| 354 | except EOFError: |
| 355 | break |
| 356 | datum = data[0] |
| 357 | if datum[0] == PWDRname: |
| 358 | PWDRdata['Data'] = datum[1][1] #powder data arrays |
| 359 | PWDRdata[data[2][0]] = data[2][1] #Limits |
| 360 | PWDRdata[data[3][0]] = data[3][1] #Background |
| 361 | PWDRdata[data[4][0]] = data[4][1] #Instrument parameters |
| 362 | PWDRdata[data[5][0]] = data[5][1] #Sample parameters |
| 363 | try: |
| 364 | PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing |
| 365 | except IndexError: |
| 366 | PWDRdata['Reflection lists'] = {} |
| 367 | file.close() |
| 368 | return PWDRdata |
| 369 | |
| 370 | def GetHKLFdata(GPXfile,HKLFname): |
| 371 | ''' Returns single crystal data from GSASII gpx file |
| 372 | input: |
| 373 | GPXfile = .gpx full file name |
| 374 | HKLFname = single crystal histogram name as obtained from GetHistogramNames |
| 375 | return: |
| 376 | HKLFdata = single crystal data list of reflections: for each reflection: |
| 377 | HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase] |
| 378 | ''' |
| 379 | file = open(GPXfile,'rb') |
| 380 | HKLFdata = [] |
| 381 | while True: |
| 382 | try: |
| 383 | data = cPickle.load(file) |
| 384 | except EOFError: |
| 385 | break |
| 386 | datum = data[0] |
| 387 | if datum[0] == HKLFname: |
| 388 | HKLFdata = datum[1:][0] |
| 389 | file.close() |
| 390 | return HKLFdata |
| 391 | |
| 392 | def ShowBanner(): |
| 393 | print 80*'*' |
| 394 | print ' General Structure Analysis System-II Crystal Structure Refinement' |
| 395 | print ' by Robert B. Von Dreele & Brian H. Toby' |
| 396 | print ' Argonne National Laboratory(C), 2010' |
| 397 | print ' This product includes software developed by the UChicago Argonne, LLC,' |
| 398 | print ' as Operator of Argonne National Laboratory.' |
| 399 | print 80*'*','\n' |
| 400 | |
| 401 | def ShowControls(Controls): |
| 402 | print ' Least squares controls:' |
| 403 | print ' Refinement type: ',Controls['deriv type'] |
| 404 | if 'Hessian' in Controls['deriv type']: |
| 405 | print ' Maximum number of cycles:',Controls['max cyc'] |
| 406 | print ' Maximum histogram processes:',Controls['max Hprocess'] |
| 407 | print ' Maximum reflection processes:',Controls['max Rprocess'] |
| 408 | else: |
| 409 | print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M']) |
| 410 | print ' Initial shift factor: ','%.3f'%(Controls['shift factor']) |
| 411 | |
| 412 | def GetFFtable(General): |
| 413 | ''' returns a dictionary of form factor data for atom types found in General |
| 414 | input: |
| 415 | General = dictionary of phase info.; includes AtomTypes |
| 416 | return: |
| 417 | FFtable = dictionary of form factor data; key is atom type |
| 418 | ''' |
| 419 | atomTypes = General['AtomTypes'] |
| 420 | FFtable = {} |
| 421 | for El in atomTypes: |
| 422 | FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0]) |
| 423 | for item in FFs: |
| 424 | if item['Symbol'] == El.upper(): |
| 425 | FFtable[El] = item |
| 426 | return FFtable |
| 427 | |
| 428 | def GetBLtable(General): |
| 429 | ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General |
| 430 | input: |
| 431 | General = dictionary of phase info.; includes AtomTypes & Isotopes |
| 432 | return: |
| 433 | BLtable = dictionary of scattering length data; key is atom type |
| 434 | ''' |
| 435 | atomTypes = General['AtomTypes'] |
| 436 | BLtable = {} |
| 437 | isotopes = General['Isotopes'] |
| 438 | isotope = General['Isotope'] |
| 439 | for El in atomTypes: |
| 440 | BLtable[El] = [isotope[El],isotopes[El][isotope[El]]] |
| 441 | return BLtable |
| 442 | |
| 443 | def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary): |
| 444 | if SGLaue in ['-1','2/m','mmm']: |
| 445 | return #no Pawley symmetry required constraints |
| 446 | for i,varyI in enumerate(pawleyVary): |
| 447 | refI = int(varyI.split(':')[-1]) |
| 448 | ih,ik,il = PawleyRef[refI][:3] |
| 449 | for varyJ in pawleyVary[0:i]: |
| 450 | refJ = int(varyJ.split(':')[-1]) |
| 451 | jh,jk,jl = PawleyRef[refJ][:3] |
| 452 | if SGLaue in ['4/m','4/mmm']: |
| 453 | isum = ih**2+ik**2 |
| 454 | jsum = jh**2+jk**2 |
| 455 | if abs(il) == abs(jl) and isum == jsum: |
| 456 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
| 457 | elif SGLaue in ['3R','3mR']: |
| 458 | isum = ih**2+ik**2+il**2 |
| 459 | jsum = jh**2+jk**2*jl**2 |
| 460 | isum2 = ih*ik+ih*il+ik*il |
| 461 | jsum2 = jh*jk+jh*jl+jk*jl |
| 462 | if isum == jsum and isum2 == jsum2: |
| 463 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
| 464 | elif SGLaue in ['3','3m1','31m','6/m','6/mmm']: |
| 465 | isum = ih**2+ik**2+ih*ik |
| 466 | jsum = jh**2+jk**2+jh*jk |
| 467 | if abs(il) == abs(jl) and isum == jsum: |
| 468 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
| 469 | elif SGLaue in ['m3','m3m']: |
| 470 | isum = ih**2+ik**2+il**2 |
| 471 | jsum = jh**2+jk**2+jl**2 |
| 472 | if isum == jsum: |
| 473 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
| 474 | |
| 475 | def cellVary(pfx,SGData): |
| 476 | if SGData['SGLaue'] in ['-1',]: |
| 477 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
| 478 | elif SGData['SGLaue'] in ['2/m',]: |
| 479 | if SGData['SGUniq'] == 'a': |
| 480 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3'] |
| 481 | elif SGData['SGUniq'] == 'b': |
| 482 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4'] |
| 483 | else: |
| 484 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5'] |
| 485 | elif SGData['SGLaue'] in ['mmm',]: |
| 486 | return [pfx+'A0',pfx+'A1',pfx+'A2'] |
| 487 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
| 488 | return [pfx+'A0',pfx+'A2'] |
| 489 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
| 490 | return [pfx+'A0',pfx+'A2'] |
| 491 | elif SGData['SGLaue'] in ['3R', '3mR']: |
| 492 | return [pfx+'A0',pfx+'A3'] |
| 493 | elif SGData['SGLaue'] in ['m3m','m3']: |
| 494 | return [pfx+'A0',] |
| 495 | |
| 496 | def GetPhaseData(PhaseData,Print=True): |
| 497 | |
| 498 | def PrintFFtable(FFtable): |
| 499 | print '\n X-ray scattering factors:' |
| 500 | print ' Symbol fa fb fc' |
| 501 | print 99*'-' |
| 502 | for Ename in FFtable: |
| 503 | ffdata = FFtable[Ename] |
| 504 | fa = ffdata['fa'] |
| 505 | fb = ffdata['fb'] |
| 506 | print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \ |
| 507 | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc']) |
| 508 | |
| 509 | def PrintBLtable(BLtable): |
| 510 | print '\n Neutron scattering factors:' |
| 511 | print ' Symbol isotope mass b resonant terms' |
| 512 | print 99*'-' |
| 513 | for Ename in BLtable: |
| 514 | bldata = BLtable[Ename] |
| 515 | isotope = bldata[0] |
| 516 | mass = bldata[1][0] |
| 517 | blen = bldata[1][1] |
| 518 | bres = [] |
| 519 | if len(bldata[1]) > 2: |
| 520 | bres = bldata[1][2:] |
| 521 | line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen) |
| 522 | for item in bres: |
| 523 | line += '%10.5g'%(item) |
| 524 | print line |
| 525 | |
| 526 | def PrintAtoms(General,Atoms): |
| 527 | print '\n Atoms:' |
| 528 | line = ' name type refine? x y z '+ \ |
| 529 | ' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23' |
| 530 | if General['Type'] == 'magnetic': |
| 531 | line += ' Mx My Mz' |
| 532 | elif General['Type'] == 'macromolecular': |
| 533 | line = ' res no residue chain '+line |
| 534 | print line |
| 535 | if General['Type'] == 'nuclear': |
| 536 | print 135*'-' |
| 537 | for i,at in enumerate(Atoms): |
| 538 | line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \ |
| 539 | '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9]) |
| 540 | if at[9] == 'I': |
| 541 | line += '%8.4f'%(at[10])+48*' ' |
| 542 | else: |
| 543 | line += 8*' ' |
| 544 | for j in range(6): |
| 545 | line += '%8.4f'%(at[11+j]) |
| 546 | print line |
| 547 | |
| 548 | def PrintTexture(textureData): |
| 549 | topstr = '\n Spherical harmonics texture: Order:' + \ |
| 550 | str(textureData['Order']) |
| 551 | if textureData['Order']: |
| 552 | print topstr+' Refine? '+str(textureData['SH Coeff'][0]) |
| 553 | else: |
| 554 | print topstr |
| 555 | return |
| 556 | names = ['omega','chi','phi'] |
| 557 | line = '\n' |
| 558 | for name in names: |
| 559 | line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0]) |
| 560 | print line |
| 561 | print '\n Texture coefficients:' |
| 562 | ptlbls = ' names :' |
| 563 | ptstr = ' values:' |
| 564 | SHcoeff = textureData['SH Coeff'][1] |
| 565 | for item in SHcoeff: |
| 566 | ptlbls += '%12s'%(item) |
| 567 | ptstr += '%12.4f'%(SHcoeff[item]) |
| 568 | print ptlbls |
| 569 | print ptstr |
| 570 | |
| 571 | if Print: print ' Phases:' |
| 572 | phaseVary = [] |
| 573 | phaseDict = {} |
| 574 | phaseConstr = {} |
| 575 | pawleyLookup = {} |
| 576 | FFtables = {} #scattering factors - xrays |
| 577 | BLtables = {} # neutrons |
| 578 | Natoms = {} |
| 579 | AtMults = {} |
| 580 | AtIA = {} |
| 581 | shModels = ['cylindrical','none','shear - 2/m','rolling - mmm'] |
| 582 | SamSym = dict(zip(shModels,['0','-1','2/m','mmm'])) |
| 583 | for name in PhaseData: |
| 584 | General = PhaseData[name]['General'] |
| 585 | pId = PhaseData[name]['pId'] |
| 586 | pfx = str(pId)+'::' |
| 587 | FFtable = GetFFtable(General) |
| 588 | BLtable = GetBLtable(General) |
| 589 | FFtables.update(FFtable) |
| 590 | BLtables.update(BLtable) |
| 591 | Atoms = PhaseData[name]['Atoms'] |
| 592 | try: |
| 593 | PawleyRef = PhaseData[name]['Pawley ref'] |
| 594 | except KeyError: |
| 595 | PawleyRef = [] |
| 596 | SGData = General['SGData'] |
| 597 | SGtext = G2spc.SGPrint(SGData) |
| 598 | cell = General['Cell'] |
| 599 | A = G2lat.cell2A(cell[1:7]) |
| 600 | phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]}) |
| 601 | if cell[0]: |
| 602 | phaseVary += cellVary(pfx,SGData) |
| 603 | Natoms[pfx] = 0 |
| 604 | if Atoms: |
| 605 | if General['Type'] == 'nuclear': |
| 606 | Natoms[pfx] = len(Atoms) |
| 607 | for i,at in enumerate(Atoms): |
| 608 | phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8], |
| 609 | pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5], |
| 610 | pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z |
| 611 | pfx+'AI/A:'+str(i):at[9],}) |
| 612 | if at[9] == 'I': |
| 613 | phaseDict[pfx+'AUiso:'+str(i)] = at[10] |
| 614 | else: |
| 615 | phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13], |
| 616 | pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]}) |
| 617 | if 'F' in at[2]: |
| 618 | phaseVary.append(pfx+'Afrac:'+str(i)) |
| 619 | if 'X' in at[2]: |
| 620 | xId,xCoef = G2spc.GetCSxinel(at[7]) |
| 621 | delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)] |
| 622 | for j in range(3): |
| 623 | if xId[j] > 0: |
| 624 | phaseVary.append(delnames[j]) |
| 625 | for k in range(j): |
| 626 | if xId[j] == xId[k]: |
| 627 | G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) |
| 628 | if 'U' in at[2]: |
| 629 | if at[9] == 'I': |
| 630 | phaseVary.append(pfx+'AUiso:'+str(i)) |
| 631 | else: |
| 632 | uId,uCoef = G2spc.GetCSuinel(at[7])[:2] |
| 633 | names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i), |
| 634 | pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)] |
| 635 | for j in range(6): |
| 636 | if uId[j] > 0: |
| 637 | phaseVary.append(names[j]) |
| 638 | for k in range(j): |
| 639 | if uId[j] == uId[k]: |
| 640 | G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),)) |
| 641 | # elif General['Type'] == 'magnetic': |
| 642 | # elif General['Type'] == 'macromolecular': |
| 643 | |
| 644 | |
| 645 | textureData = General['SH Texture'] |
| 646 | if textureData['Order']: |
| 647 | phaseDict[pfx+'SHorder'] = textureData['Order'] |
| 648 | phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']] |
| 649 | for name in ['omega','chi','phi']: |
| 650 | phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1] |
| 651 | if textureData['Sample '+name][0]: |
| 652 | phaseVary.append(pfx+'SH '+name) |
| 653 | for name in textureData['SH Coeff'][1]: |
| 654 | phaseDict[pfx+name] = textureData['SH Coeff'][1][name] |
| 655 | if textureData['SH Coeff'][0]: |
| 656 | phaseVary.append(pfx+name) |
| 657 | |
| 658 | if Print: |
| 659 | print '\n Phase name: ',General['Name'] |
| 660 | print 135*'-' |
| 661 | PrintFFtable(FFtable) |
| 662 | PrintBLtable(BLtable) |
| 663 | print '' |
| 664 | for line in SGtext: print line |
| 665 | PrintAtoms(General,Atoms) |
| 666 | print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \ |
| 667 | ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \ |
| 668 | '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0] |
| 669 | PrintTexture(textureData) |
| 670 | |
| 671 | elif PawleyRef: |
| 672 | pawleyVary = [] |
| 673 | for i,refl in enumerate(PawleyRef): |
| 674 | phaseDict[pfx+'PWLref:'+str(i)] = refl[6] |
| 675 | pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i |
| 676 | if refl[5]: |
| 677 | pawleyVary.append(pfx+'PWLref:'+str(i)) |
| 678 | GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence |
| 679 | phaseVary += pawleyVary |
| 680 | |
| 681 | return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables |
| 682 | |
| 683 | def cellFill(pfx,SGData,parmDict,sigDict): |
| 684 | if SGData['SGLaue'] in ['-1',]: |
| 685 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
| 686 | parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']] |
| 687 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
| 688 | sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']] |
| 689 | elif SGData['SGLaue'] in ['2/m',]: |
| 690 | if SGData['SGUniq'] == 'a': |
| 691 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
| 692 | parmDict[pfx+'A3'],0,0] |
| 693 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
| 694 | sigDict[pfx+'A3'],0,0] |
| 695 | elif SGData['SGUniq'] == 'b': |
| 696 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
| 697 | 0,parmDict[pfx+'A4'],0] |
| 698 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
| 699 | 0,sigDict[pfx+'A4'],0] |
| 700 | else: |
| 701 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
| 702 | 0,0,parmDict[pfx+'A5']] |
| 703 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
| 704 | 0,0,sigDict[pfx+'A5']] |
| 705 | elif SGData['SGLaue'] in ['mmm',]: |
| 706 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0] |
| 707 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0] |
| 708 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
| 709 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0] |
| 710 | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
| 711 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
| 712 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'], |
| 713 | parmDict[pfx+'A0'],0,0] |
| 714 | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
| 715 | elif SGData['SGLaue'] in ['3R', '3mR']: |
| 716 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'], |
| 717 | parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']] |
| 718 | sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0] |
| 719 | elif SGData['SGLaue'] in ['m3m','m3']: |
| 720 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0] |
| 721 | sigA = [sigDict[pfx+'A0'],0,0,0,0,0] |
| 722 | return A,sigA |
| 723 | |
| 724 | def getCellEsd(pfx,SGData,A,covData): |
| 725 | dpr = 180./np.pi |
| 726 | rVsq = G2lat.calc_rVsq(A) |
| 727 | G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors |
| 728 | cell = np.array(G2lat.Gmat2cell(g)) #real cell |
| 729 | cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell |
| 730 | scos = cosd(cellst[3:6]) |
| 731 | ssin = sind(cellst[3:6]) |
| 732 | scot = scos/ssin |
| 733 | rcos = cosd(cell[3:6]) |
| 734 | rsin = sind(cell[3:6]) |
| 735 | rcot = rcos/rsin |
| 736 | RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
| 737 | varyList = covData['varyList'] |
| 738 | covMatrix = covData['covMatrix'] |
| 739 | vcov = G2mth.getVCov(RMnames,varyList,covMatrix) |
| 740 | Ax = np.array(A) |
| 741 | Ax[3:] /= 2. |
| 742 | drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2, |
| 743 | Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]]) |
| 744 | srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T)) |
| 745 | Vol = 1/np.sqrt(rVsq) |
| 746 | sigVol = Vol**3*np.sqrt(srcvlsq)/2. |
| 747 | R123 = Ax[0]*Ax[1]*Ax[2] |
| 748 | dsasdg = np.zeros((3,6)) |
| 749 | dadg = np.zeros((6,6)) |
| 750 | for i0 in range(3): #0 1 2 |
| 751 | i1 = (i0+1)%3 #1 2 0 |
| 752 | i2 = (i1+1)%3 #2 0 1 |
| 753 | i3 = 5-i2 #3 5 4 |
| 754 | i4 = 5-i1 #4 3 5 |
| 755 | i5 = 5-i0 #5 4 3 |
| 756 | dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1] |
| 757 | dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2] |
| 758 | dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2]) |
| 759 | denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2) |
| 760 | denom = np.sqrt(denmsq) |
| 761 | dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2) |
| 762 | dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2) |
| 763 | dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2) |
| 764 | dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2) |
| 765 | dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4]) |
| 766 | dadg[i5][i5] = -Ax[i0]/denom |
| 767 | for i0 in range(3): |
| 768 | i1 = (i0+1)%3 |
| 769 | i2 = (i1+1)%3 |
| 770 | i3 = 5-i2 |
| 771 | for ij in range(6): |
| 772 | dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1]) |
| 773 | if ij == i0: |
| 774 | dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0] |
| 775 | dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr |
| 776 | sigMat = np.inner(dadg,np.inner(vcov,dadg.T)) |
| 777 | var = np.diag(sigMat) |
| 778 | CS = np.where(var>0.,np.sqrt(var),0.) |
| 779 | cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order |
| 780 | return cellSig |
| 781 | |
| 782 | def SetPhaseData(parmDict,sigDict,Phases,covData): |
| 783 | |
| 784 | def PrintAtomsAndSig(General,Atoms,atomsSig): |
| 785 | print '\n Atoms:' |
| 786 | line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23' |
| 787 | if General['Type'] == 'magnetic': |
| 788 | line += ' Mx My Mz' |
| 789 | elif General['Type'] == 'macromolecular': |
| 790 | line = ' res no residue chain '+line |
| 791 | print line |
| 792 | if General['Type'] == 'nuclear': |
| 793 | print 135*'-' |
| 794 | fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f', |
| 795 | 11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'} |
| 796 | noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']} |
| 797 | for i,at in enumerate(Atoms): |
| 798 | name = fmt[0]%(at[0])+fmt[1]%(at[1])+':' |
| 799 | valstr = ' values:' |
| 800 | sigstr = ' sig :' |
| 801 | for ind in [3,4,5,6]: |
| 802 | sigind = str(i)+':'+str(ind) |
| 803 | valstr += fmt[ind]%(at[ind]) |
| 804 | if sigind in atomsSig: |
| 805 | sigstr += fmt[ind]%(atomsSig[sigind]) |
| 806 | else: |
| 807 | sigstr += noFXsig[ind][1]%(noFXsig[ind][0]) |
| 808 | if at[9] == 'I': |
| 809 | valstr += fmt[10]%(at[10]) |
| 810 | if str(i)+':10' in atomsSig: |
| 811 | sigstr += fmt[10]%(atomsSig[str(i)+':10']) |
| 812 | else: |
| 813 | sigstr += 8*' ' |
| 814 | else: |
| 815 | valstr += 8*' ' |
| 816 | sigstr += 8*' ' |
| 817 | for ind in [11,12,13,14,15,16]: |
| 818 | sigind = str(i)+':'+str(ind) |
| 819 | valstr += fmt[ind]%(at[ind]) |
| 820 | if sigind in atomsSig: |
| 821 | sigstr += fmt[ind]%(atomsSig[sigind]) |
| 822 | else: |
| 823 | sigstr += 8*' ' |
| 824 | print name |
| 825 | print valstr |
| 826 | print sigstr |
| 827 | |
| 828 | def PrintSHtextureAndSig(textureData,SHtextureSig): |
| 829 | print '\n Spherical harmonics texture: Order:' + str(textureData['Order']) |
| 830 | names = ['omega','chi','phi'] |
| 831 | namstr = ' names :' |
| 832 | ptstr = ' values:' |
| 833 | sigstr = ' esds :' |
| 834 | for name in names: |
| 835 | namstr += '%12s'%(name) |
| 836 | ptstr += '%12.3f'%(textureData['Sample '+name][1]) |
| 837 | if 'Sample '+name in SHtextureSig: |
| 838 | sigstr += '%12.3f'%(SHtextureSig['Sample '+name]) |
490 | | |
491 | | last = max(last,int(name[-2].strip('bak'))) |
492 | | |
493 | | GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx') |
494 | | |
495 | | dfu.copy_file(GPXfile,GPXback) |
496 | | |
497 | | return GPXback |
498 | | |
499 | | |
500 | | |
501 | | def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True): |
502 | | |
503 | | ''' Updates gpxfile from all histograms that are found in any phase |
504 | | |
505 | | and any phase that used a histogram |
506 | | |
| 840 | sigstr += 12*' ' |
| 841 | print namstr |
| 842 | print ptstr |
| 843 | print sigstr |
| 844 | print '\n Texture coefficients:' |
| 845 | namstr = ' names :' |
| 846 | ptstr = ' values:' |
| 847 | sigstr = ' esds :' |
| 848 | SHcoeff = textureData['SH Coeff'][1] |
| 849 | for name in SHcoeff: |
| 850 | namstr += '%12s'%(name) |
| 851 | ptstr += '%12.3f'%(SHcoeff[name]) |
| 852 | if name in SHtextureSig: |
| 853 | sigstr += '%12.3f'%(SHtextureSig[name]) |
| 854 | else: |
| 855 | sigstr += 12*' ' |
| 856 | print namstr |
| 857 | print ptstr |
| 858 | print sigstr |
| 859 | |
| 860 | |
| 861 | print '\n Phases:' |
| 862 | for phase in Phases: |
| 863 | print ' Result for phase: ',phase |
| 864 | Phase = Phases[phase] |
| 865 | General = Phase['General'] |
| 866 | SGData = General['SGData'] |
| 867 | Atoms = Phase['Atoms'] |
| 868 | cell = General['Cell'] |
| 869 | pId = Phase['pId'] |
| 870 | pfx = str(pId)+'::' |
| 871 | if cell[0]: |
| 872 | A,sigA = cellFill(pfx,SGData,parmDict,sigDict) |
| 873 | cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol |
| 874 | print ' Reciprocal metric tensor: ' |
| 875 | ptfmt = "%15.9f" |
| 876 | names = ['A11','A22','A33','A12','A13','A23'] |
| 877 | namstr = ' names :' |
| 878 | ptstr = ' values:' |
| 879 | sigstr = ' esds :' |
| 880 | for name,a,siga in zip(names,A,sigA): |
| 881 | namstr += '%15s'%(name) |
| 882 | ptstr += ptfmt%(a) |
| 883 | if siga: |
| 884 | sigstr += ptfmt%(siga) |
| 885 | else: |
| 886 | sigstr += 15*' ' |
| 887 | print namstr |
| 888 | print ptstr |
| 889 | print sigstr |
| 890 | cell[1:7] = G2lat.A2cell(A) |
| 891 | cell[7] = G2lat.calc_V(A) |
| 892 | print ' New unit cell:' |
| 893 | ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"] |
| 894 | names = ['a','b','c','alpha','beta','gamma','Volume'] |
| 895 | namstr = ' names :' |
| 896 | ptstr = ' values:' |
| 897 | sigstr = ' esds :' |
| 898 | for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig): |
| 899 | namstr += '%12s'%(name) |
| 900 | ptstr += fmt%(a) |
| 901 | if siga: |
| 902 | sigstr += fmt%(siga) |
| 903 | else: |
| 904 | sigstr += 12*' ' |
| 905 | print namstr |
| 906 | print ptstr |
| 907 | print sigstr |
| 908 | |
| 909 | if 'Pawley' in Phase['General']['Type']: |
| 910 | pawleyRef = Phase['Pawley ref'] |
| 911 | for i,refl in enumerate(pawleyRef): |
| 912 | key = pfx+'PWLref:'+str(i) |
| 913 | refl[6] = abs(parmDict[key]) #suppress negative Fsq |
| 914 | if key in sigDict: |
| 915 | refl[7] = sigDict[key] |
| 916 | else: |
| 917 | refl[7] = 0 |
| 918 | else: |
| 919 | atomsSig = {} |
| 920 | if General['Type'] == 'nuclear': |
| 921 | for i,at in enumerate(Atoms): |
| 922 | names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i), |
| 923 | 10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i), |
| 924 | 14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)} |
| 925 | for ind in [3,4,5,6]: |
| 926 | at[ind] = parmDict[names[ind]] |
| 927 | if ind in [3,4,5]: |
| 928 | name = names[ind].replace('A','dA') |
| 929 | else: |
| 930 | name = names[ind] |
| 931 | if name in sigDict: |
| 932 | atomsSig[str(i)+':'+str(ind)] = sigDict[name] |
| 933 | if at[9] == 'I': |
| 934 | at[10] = parmDict[names[10]] |
| 935 | if names[10] in sigDict: |
| 936 | atomsSig[str(i)+':10'] = sigDict[names[10]] |
| 937 | else: |
| 938 | for ind in [11,12,13,14,15,16]: |
| 939 | at[ind] = parmDict[names[ind]] |
| 940 | if names[ind] in sigDict: |
| 941 | atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]] |
| 942 | PrintAtomsAndSig(General,Atoms,atomsSig) |
| 943 | |
| 944 | textureData = General['SH Texture'] |
| 945 | if textureData['Order']: |
| 946 | SHtextureSig = {} |
| 947 | for name in ['omega','chi','phi']: |
| 948 | aname = pfx+'SH '+name |
| 949 | textureData['Sample '+name][1] = parmDict[aname] |
| 950 | if aname in sigDict: |
| 951 | SHtextureSig['Sample '+name] = sigDict[aname] |
| 952 | for name in textureData['SH Coeff'][1]: |
| 953 | aname = pfx+name |
| 954 | textureData['SH Coeff'][1][name] = parmDict[aname] |
| 955 | if aname in sigDict: |
| 956 | SHtextureSig[name] = sigDict[aname] |
| 957 | PrintSHtextureAndSig(textureData,SHtextureSig) |
| 958 | |
| 959 | def GetHistogramPhaseData(Phases,Histograms,Print=True): |
| 960 | |
| 961 | def PrintSize(hapData): |
| 962 | if hapData[0] in ['isotropic','uniaxial']: |
| 963 | line = '\n Size model : %9s'%(hapData[0]) |
| 964 | line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
| 965 | if hapData[0] == 'uniaxial': |
| 966 | line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
| 967 | print line |
| 968 | else: |
| 969 | print '\n Size model : %s'%(hapData[0]) |
| 970 | Snames = ['S11','S22','S33','S12','S13','S23'] |
| 971 | ptlbls = ' names :' |
| 972 | ptstr = ' values:' |
| 973 | varstr = ' refine:' |
| 974 | for i,name in enumerate(Snames): |
| 975 | ptlbls += '%12s' % (name) |
| 976 | ptstr += '%12.6f' % (hapData[4][i]) |
| 977 | varstr += '%12s' % (str(hapData[5][i])) |
| 978 | print ptlbls |
| 979 | print ptstr |
| 980 | print varstr |
| 981 | |
| 982 | def PrintMuStrain(hapData,SGData): |
| 983 | if hapData[0] in ['isotropic','uniaxial']: |
| 984 | line = '\n Mustrain model: %9s'%(hapData[0]) |
| 985 | line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
| 986 | if hapData[0] == 'uniaxial': |
| 987 | line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
| 988 | print line |
| 989 | else: |
| 990 | print '\n Mustrain model: %s'%(hapData[0]) |
| 991 | Snames = G2spc.MustrainNames(SGData) |
| 992 | ptlbls = ' names :' |
| 993 | ptstr = ' values:' |
| 994 | varstr = ' refine:' |
| 995 | for i,name in enumerate(Snames): |
| 996 | ptlbls += '%12s' % (name) |
| 997 | ptstr += '%12.6f' % (hapData[4][i]) |
| 998 | varstr += '%12s' % (str(hapData[5][i])) |
| 999 | print ptlbls |
| 1000 | print ptstr |
| 1001 | print varstr |
| 1002 | |
| 1003 | def PrintHStrain(hapData,SGData): |
| 1004 | print '\n Hydrostatic/elastic strain: ' |
| 1005 | Hsnames = G2spc.HStrainNames(SGData) |
| 1006 | ptlbls = ' names :' |
| 1007 | ptstr = ' values:' |
| 1008 | varstr = ' refine:' |
| 1009 | for i,name in enumerate(Hsnames): |
| 1010 | ptlbls += '%12s' % (name) |
| 1011 | ptstr += '%12.6f' % (hapData[0][i]) |
| 1012 | varstr += '%12s' % (str(hapData[1][i])) |
| 1013 | print ptlbls |
| 1014 | print ptstr |
| 1015 | print varstr |
| 1016 | |
| 1017 | def PrintSHPO(hapData): |
| 1018 | print '\n Spherical harmonics preferred orientation: Order:' + \ |
| 1019 | str(hapData[4])+' Refine? '+str(hapData[2]) |
| 1020 | ptlbls = ' names :' |
| 1021 | ptstr = ' values:' |
| 1022 | for item in hapData[5]: |
| 1023 | ptlbls += '%12s'%(item) |
| 1024 | ptstr += '%12.3f'%(hapData[5][item]) |
| 1025 | print ptlbls |
| 1026 | print ptstr |
| 1027 | |
| 1028 | hapDict = {} |
| 1029 | hapVary = [] |
| 1030 | controlDict = {} |
| 1031 | poType = {} |
| 1032 | poAxes = {} |
| 1033 | spAxes = {} |
| 1034 | spType = {} |
| 1035 | |
| 1036 | for phase in Phases: |
| 1037 | HistoPhase = Phases[phase]['Histograms'] |
| 1038 | SGData = Phases[phase]['General']['SGData'] |
| 1039 | cell = Phases[phase]['General']['Cell'][1:7] |
| 1040 | A = G2lat.cell2A(cell) |
| 1041 | pId = Phases[phase]['pId'] |
| 1042 | histoList = HistoPhase.keys() |
| 1043 | histoList.sort() |
| 1044 | for histogram in histoList: |
| 1045 | try: |
| 1046 | Histogram = Histograms[histogram] |
| 1047 | except KeyError: |
| 1048 | #skip if histogram not included e.g. in a sequential refinement |
| 1049 | continue |
| 1050 | hapData = HistoPhase[histogram] |
| 1051 | hId = Histogram['hId'] |
| 1052 | limits = Histogram['Limits'][1] |
| 1053 | inst = Histogram['Instrument Parameters'] |
| 1054 | inst = dict(zip(inst[3],inst[1])) |
| 1055 | Zero = inst['Zero'] |
| 1056 | if 'C' in inst['Type']: |
| 1057 | try: |
| 1058 | wave = inst['Lam'] |
| 1059 | except KeyError: |
| 1060 | wave = inst['Lam1'] |
| 1061 | dmin = wave/(2.0*sind(limits[1]/2.0)) |
| 1062 | pfx = str(pId)+':'+str(hId)+':' |
| 1063 | for item in ['Scale','Extinction']: |
| 1064 | hapDict[pfx+item] = hapData[item][0] |
| 1065 | if hapData[item][1]: |
| 1066 | hapVary.append(pfx+item) |
| 1067 | names = G2spc.HStrainNames(SGData) |
| 1068 | for i,name in enumerate(names): |
| 1069 | hapDict[pfx+name] = hapData['HStrain'][0][i] |
| 1070 | if hapData['HStrain'][1][i]: |
| 1071 | hapVary.append(pfx+name) |
| 1072 | controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0] |
| 1073 | if hapData['Pref.Ori.'][0] == 'MD': |
| 1074 | hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1] |
| 1075 | controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3] |
| 1076 | if hapData['Pref.Ori.'][2]: |
| 1077 | hapVary.append(pfx+'MD') |
| 1078 | else: #'SH' spherical harmonics |
| 1079 | controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4] |
| 1080 | controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5]) |
| 1081 | for item in hapData['Pref.Ori.'][5]: |
| 1082 | hapDict[pfx+item] = hapData['Pref.Ori.'][5][item] |
| 1083 | if hapData['Pref.Ori.'][2]: |
| 1084 | hapVary.append(pfx+item) |
| 1085 | for item in ['Mustrain','Size']: |
| 1086 | controlDict[pfx+item+'Type'] = hapData[item][0] |
| 1087 | if hapData[item][0] in ['isotropic','uniaxial']: |
| 1088 | hapDict[pfx+item+':i'] = hapData[item][1][0] |
| 1089 | if hapData[item][2][0]: |
| 1090 | hapVary.append(pfx+item+':i') |
| 1091 | if hapData[item][0] == 'uniaxial': |
| 1092 | controlDict[pfx+item+'Axis'] = hapData[item][3] |
| 1093 | hapDict[pfx+item+':a'] = hapData[item][1][1] |
| 1094 | if hapData[item][2][1]: |
| 1095 | hapVary.append(pfx+item+':a') |
| 1096 | else: #generalized for mustrain or ellipsoidal for size |
| 1097 | if item == 'Mustrain': |
| 1098 | names = G2spc.MustrainNames(SGData) |
| 1099 | pwrs = [] |
| 1100 | for name in names: |
| 1101 | h,k,l = name[1:] |
| 1102 | pwrs.append([int(h),int(k),int(l)]) |
| 1103 | controlDict[pfx+'MuPwrs'] = pwrs |
| 1104 | for i in range(len(hapData[item][4])): |
| 1105 | sfx = ':'+str(i) |
| 1106 | hapDict[pfx+item+sfx] = hapData[item][4][i] |
| 1107 | if hapData[item][5][i]: |
| 1108 | hapVary.append(pfx+item+sfx) |
| 1109 | |
| 1110 | if Print: |
| 1111 | print '\n Phase: ',phase,' in histogram: ',histogram |
| 1112 | print 135*'-' |
| 1113 | print ' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1] |
| 1114 | print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1] |
| 1115 | if hapData['Pref.Ori.'][0] == 'MD': |
| 1116 | Ax = hapData['Pref.Ori.'][3] |
| 1117 | print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \ |
| 1118 | ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2]) |
| 1119 | else: #'SH' for spherical harmonics |
| 1120 | PrintSHPO(hapData['Pref.Ori.']) |
| 1121 | PrintSize(hapData['Size']) |
| 1122 | PrintMuStrain(hapData['Mustrain'],SGData) |
| 1123 | PrintHStrain(hapData['HStrain'],SGData) |
| 1124 | HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A)) |
| 1125 | refList = [] |
| 1126 | for h,k,l,d in HKLd: |
| 1127 | ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData) |
| 1128 | if ext: |
| 1129 | continue |
| 1130 | if 'C' in inst['Type']: |
| 1131 | pos = 2.0*asind(wave/(2.0*d))+Zero |
| 1132 | if limits[0] < pos < limits[1]: |
| 1133 | refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}]) |
| 1134 | #last item should contain form factor values by atom type |
| 1135 | else: |
| 1136 | raise ValueError |
| 1137 | Histogram['Reflection Lists'][phase] = refList |
| 1138 | return hapVary,hapDict,controlDict |
| 1139 | |
| 1140 | def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True): |
| 1141 | |
| 1142 | def PrintSizeAndSig(hapData,sizeSig): |
| 1143 | line = '\n Size model: %9s'%(hapData[0]) |
| 1144 | refine = False |
| 1145 | if hapData[0] in ['isotropic','uniaxial']: |
| 1146 | line += ' equatorial:%12.3f'%(hapData[1][0]) |
| 1147 | if sizeSig[0][0]: |
| 1148 | line += ', sig: %8.3f'%(sizeSig[0][0]) |
| 1149 | refine = True |
| 1150 | if hapData[0] == 'uniaxial': |
| 1151 | line += ' axial:%12.3f'%(hapData[1][1]) |
| 1152 | if sizeSig[0][1]: |
| 1153 | refine = True |
| 1154 | line += ', sig: %8.3f'%(sizeSig[0][1]) |
| 1155 | if refine: |
| 1156 | print line |
| 1157 | else: |
| 1158 | Snames = ['S11','S22','S33','S12','S13','S23'] |
| 1159 | ptlbls = ' name :' |
| 1160 | ptstr = ' value :' |
| 1161 | sigstr = ' sig :' |
| 1162 | for i,name in enumerate(Snames): |
| 1163 | ptlbls += '%12s' % (name) |
| 1164 | ptstr += '%12.6f' % (hapData[4][i]) |
| 1165 | if sizeSig[1][i]: |
| 1166 | refine = True |
| 1167 | sigstr += '%12.6f' % (sizeSig[1][i]) |
| 1168 | else: |
| 1169 | sigstr += 12*' ' |
| 1170 | if refine: |
| 1171 | print line |
| 1172 | print ptlbls |
| 1173 | print ptstr |
| 1174 | print sigstr |
| 1175 | |
| 1176 | def PrintMuStrainAndSig(hapData,mustrainSig,SGData): |
| 1177 | line = '\n Mustrain model: %9s'%(hapData[0]) |
| 1178 | refine = False |
| 1179 | if hapData[0] in ['isotropic','uniaxial']: |
| 1180 | line += ' equatorial:%12.1f'%(hapData[1][0]) |
| 1181 | if mustrainSig[0][0]: |
| 1182 | line += ', sig: %8.1f'%(mustrainSig[0][0]) |
| 1183 | refine = True |
| 1184 | if hapData[0] == 'uniaxial': |
| 1185 | line += ' axial:%12.1f'%(hapData[1][1]) |
| 1186 | if mustrainSig[0][1]: |
| 1187 | line += ', sig: %8.1f'%(mustrainSig[0][1]) |
| 1188 | if refine: |
| 1189 | print line |
| 1190 | else: |
| 1191 | Snames = G2spc.MustrainNames(SGData) |
| 1192 | ptlbls = ' name :' |
| 1193 | ptstr = ' value :' |
| 1194 | sigstr = ' sig :' |
| 1195 | for i,name in enumerate(Snames): |
| 1196 | ptlbls += '%12s' % (name) |
| 1197 | ptstr += '%12.6f' % (hapData[4][i]) |
| 1198 | if mustrainSig[1][i]: |
| 1199 | refine = True |
| 1200 | sigstr += '%12.6f' % (mustrainSig[1][i]) |
| 1201 | else: |
| 1202 | sigstr += 12*' ' |
| 1203 | if refine: |
| 1204 | print line |
| 1205 | print ptlbls |
| 1206 | print ptstr |
| 1207 | print sigstr |
| 1208 | |
| 1209 | def PrintHStrainAndSig(hapData,strainSig,SGData): |
| 1210 | Hsnames = G2spc.HStrainNames(SGData) |
| 1211 | ptlbls = ' name :' |
| 1212 | ptstr = ' value :' |
| 1213 | sigstr = ' sig :' |
| 1214 | refine = False |
| 1215 | for i,name in enumerate(Hsnames): |
| 1216 | ptlbls += '%12s' % (name) |
| 1217 | ptstr += '%12.6g' % (hapData[0][i]) |
| 1218 | if name in strainSig: |
| 1219 | refine = True |
| 1220 | sigstr += '%12.6g' % (strainSig[name]) |
| 1221 | else: |
| 1222 | sigstr += 12*' ' |
| 1223 | if refine: |
| 1224 | print '\n Hydrostatic/elastic strain: ' |
| 1225 | print ptlbls |
| 1226 | print ptstr |
| 1227 | print sigstr |
| 1228 | |
| 1229 | def PrintSHPOAndSig(hapData,POsig): |
| 1230 | print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4]) |
| 1231 | ptlbls = ' names :' |
| 1232 | ptstr = ' values:' |
| 1233 | sigstr = ' sig :' |
| 1234 | for item in hapData[5]: |
| 1235 | ptlbls += '%12s'%(item) |
| 1236 | ptstr += '%12.3f'%(hapData[5][item]) |
| 1237 | if item in POsig: |
| 1238 | sigstr += '%12.3f'%(POsig[item]) |
| 1239 | else: |
| 1240 | sigstr += 12*' ' |
| 1241 | print ptlbls |
| 1242 | print ptstr |
| 1243 | print sigstr |
| 1244 | |
| 1245 | for phase in Phases: |
| 1246 | HistoPhase = Phases[phase]['Histograms'] |
| 1247 | SGData = Phases[phase]['General']['SGData'] |
| 1248 | pId = Phases[phase]['pId'] |
| 1249 | histoList = HistoPhase.keys() |
| 1250 | histoList.sort() |
| 1251 | for histogram in histoList: |
| 1252 | try: |
| 1253 | Histogram = Histograms[histogram] |
| 1254 | except KeyError: |
| 1255 | #skip if histogram not included e.g. in a sequential refinement |
| 1256 | continue |
| 1257 | print '\n Phase: ',phase,' in histogram: ',histogram |
| 1258 | print 130*'-' |
| 1259 | hapData = HistoPhase[histogram] |
| 1260 | hId = Histogram['hId'] |
| 1261 | pfx = str(pId)+':'+str(hId)+':' |
| 1262 | print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \ |
| 1263 | %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref']) |
| 1264 | |
| 1265 | PhFrExtPOSig = {} |
| 1266 | for item in ['Scale','Extinction']: |
| 1267 | hapData[item][0] = parmDict[pfx+item] |
| 1268 | if pfx+item in sigDict: |
| 1269 | PhFrExtPOSig[item] = sigDict[pfx+item] |
| 1270 | if hapData['Pref.Ori.'][0] == 'MD': |
| 1271 | hapData['Pref.Ori.'][1] = parmDict[pfx+'MD'] |
| 1272 | if pfx+'MD' in sigDict: |
| 1273 | PhFrExtPOSig['MD'] = sigDict[pfx+'MD'] |
| 1274 | else: #'SH' spherical harmonics |
| 1275 | for item in hapData['Pref.Ori.'][5]: |
| 1276 | hapData['Pref.Ori.'][5][item] = parmDict[pfx+item] |
| 1277 | if pfx+item in sigDict: |
| 1278 | PhFrExtPOSig[item] = sigDict[pfx+item] |
| 1279 | if Print: |
| 1280 | if 'Scale' in PhFrExtPOSig: |
| 1281 | print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale']) |
| 1282 | if 'Extinction' in PhFrExtPOSig: |
| 1283 | print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction']) |
| 1284 | if hapData['Pref.Ori.'][0] == 'MD': |
| 1285 | if 'MD' in PhFrExtPOSig: |
| 1286 | print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD']) |
| 1287 | else: |
| 1288 | PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig) |
| 1289 | SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]], |
| 1290 | 'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]], |
| 1291 | 'HStrain':{}} |
| 1292 | for item in ['Mustrain','Size']: |
| 1293 | if hapData[item][0] in ['isotropic','uniaxial']: |
| 1294 | hapData[item][1][0] = parmDict[pfx+item+':i'] |
| 1295 | if item == 'Size': |
| 1296 | hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0])) |
| 1297 | if pfx+item+':i' in sigDict: |
| 1298 | SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i'] |
| 1299 | if hapData[item][0] == 'uniaxial': |
| 1300 | hapData[item][1][1] = parmDict[pfx+item+':a'] |
| 1301 | if item == 'Size': |
| 1302 | hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1])) |
| 1303 | if pfx+item+':a' in sigDict: |
| 1304 | SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a'] |
| 1305 | else: #generalized for mustrain or ellipsoidal for size |
| 1306 | for i in range(len(hapData[item][4])): |
| 1307 | sfx = ':'+str(i) |
| 1308 | hapData[item][4][i] = parmDict[pfx+item+sfx] |
| 1309 | if pfx+item+sfx in sigDict: |
| 1310 | SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx] |
| 1311 | names = G2spc.HStrainNames(SGData) |
| 1312 | for i,name in enumerate(names): |
| 1313 | hapData['HStrain'][0][i] = parmDict[pfx+name] |
| 1314 | if pfx+name in sigDict: |
| 1315 | SizeMuStrSig['HStrain'][name] = sigDict[pfx+name] |
| 1316 | if Print: |
| 1317 | PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size']) |
| 1318 | PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData) |
| 1319 | PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData) |
| 1320 | |
| 1321 | def GetHistogramData(Histograms,Print=True): |
| 1322 | |
| 1323 | def GetBackgroundParms(hId,Background): |
| 1324 | Back = Background[0] |
| 1325 | Debye = Background[1] |
| 1326 | bakType,bakFlag = Back[:2] |
| 1327 | backVals = Back[3:] |
| 1328 | backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))] |
| 1329 | backDict = dict(zip(backNames,backVals)) |
| 1330 | backVary = [] |
| 1331 | if bakFlag: |
| 1332 | backVary = backNames |
| 1333 | backDict[':'+str(hId)+':nDebye'] = Debye['nDebye'] |
| 1334 | debyeDict = {} |
| 1335 | debyeList = [] |
| 1336 | for i in range(Debye['nDebye']): |
| 1337 | debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)] |
| 1338 | debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2]))) |
| 1339 | debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2]) |
| 1340 | debyeVary = [] |
| 1341 | for item in debyeList: |
| 1342 | if item[1]: |
| 1343 | debyeVary.append(item[0]) |
| 1344 | backDict.update(debyeDict) |
| 1345 | backVary += debyeVary |
| 1346 | return bakType,backDict,backVary |
| 1347 | |
| 1348 | def GetInstParms(hId,Inst): |
| 1349 | insVals,insFlags,insNames = Inst[1:4] |
| 1350 | dataType = insVals[0] |
| 1351 | instDict = {} |
| 1352 | insVary = [] |
| 1353 | pfx = ':'+str(hId)+':' |
| 1354 | for i,flag in enumerate(insFlags): |
| 1355 | insName = pfx+insNames[i] |
| 1356 | instDict[insName] = insVals[i] |
| 1357 | if flag: |
| 1358 | insVary.append(insName) |
| 1359 | instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001) |
| 1360 | instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001) |
| 1361 | instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005) |
| 1362 | return dataType,instDict,insVary |
| 1363 | |
| 1364 | def GetSampleParms(hId,Sample): |
| 1365 | sampVary = [] |
| 1366 | hfx = ':'+str(hId)+':' |
| 1367 | sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'], |
| 1368 | hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']} |
| 1369 | Type = Sample['Type'] |
| 1370 | if 'Bragg' in Type: #Bragg-Brentano |
| 1371 | for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering? |
| 1372 | sampDict[hfx+item] = Sample[item][0] |
| 1373 | if Sample[item][1]: |
| 1374 | sampVary.append(hfx+item) |
| 1375 | elif 'Debye' in Type: #Debye-Scherrer |
| 1376 | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
| 1377 | sampDict[hfx+item] = Sample[item][0] |
| 1378 | if Sample[item][1]: |
| 1379 | sampVary.append(hfx+item) |
| 1380 | return Type,sampDict,sampVary |
| 1381 | |
| 1382 | def PrintBackground(Background): |
| 1383 | Back = Background[0] |
| 1384 | Debye = Background[1] |
| 1385 | print '\n Background function: ',Back[0],' Refine?',bool(Back[1]) |
| 1386 | line = ' Coefficients: ' |
| 1387 | for i,back in enumerate(Back[3:]): |
| 1388 | line += '%10.3f'%(back) |
| 1389 | if i and not i%10: |
| 1390 | line += '\n'+15*' ' |
| 1391 | print line |
| 1392 | if Debye['nDebye']: |
| 1393 | print '\n Debye diffuse scattering coefficients' |
| 1394 | parms = ['DebyeA','DebyeR','DebyeU'] |
| 1395 | line = ' names :' |
| 1396 | for parm in parms: |
| 1397 | line += '%16s'%(parm) |
| 1398 | print line |
| 1399 | for j,term in enumerate(Debye['debyeTerms']): |
| 1400 | line = ' term'+'%2d'%(j)+':' |
| 1401 | for i in range(3): |
| 1402 | line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1])) |
| 1403 | print line |
| 1404 | |
| 1405 | def PrintInstParms(Inst): |
| 1406 | print '\n Instrument Parameters:' |
| 1407 | ptlbls = ' name :' |
| 1408 | ptstr = ' value :' |
| 1409 | varstr = ' refine:' |
| 1410 | instNames = Inst[3][1:] |
| 1411 | for i,name in enumerate(instNames): |
| 1412 | ptlbls += '%12s' % (name) |
| 1413 | ptstr += '%12.6f' % (Inst[1][i+1]) |
| 1414 | if name in ['Lam1','Lam2','Azimuth']: |
| 1415 | varstr += 12*' ' |
| 1416 | else: |
| 1417 | varstr += '%12s' % (str(bool(Inst[2][i+1]))) |
| 1418 | print ptlbls |
| 1419 | print ptstr |
| 1420 | print varstr |
| 1421 | |
| 1422 | def PrintSampleParms(Sample): |
| 1423 | print '\n Sample Parameters:' |
| 1424 | print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \ |
| 1425 | (Sample['Omega'],Sample['Chi'],Sample['Phi']) |
| 1426 | ptlbls = ' name :' |
| 1427 | ptstr = ' value :' |
| 1428 | varstr = ' refine:' |
| 1429 | if 'Bragg' in Sample['Type']: |
| 1430 | for item in ['Scale','Shift','Transparency']: |
| 1431 | ptlbls += '%14s'%(item) |
| 1432 | ptstr += '%14.4f'%(Sample[item][0]) |
| 1433 | varstr += '%14s'%(str(bool(Sample[item][1]))) |
| 1434 | |
| 1435 | elif 'Debye' in Type: #Debye-Scherrer |
| 1436 | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
| 1437 | ptlbls += '%14s'%(item) |
| 1438 | ptstr += '%14.4f'%(Sample[item][0]) |
| 1439 | varstr += '%14s'%(str(bool(Sample[item][1]))) |
| 1440 | |
| 1441 | print ptlbls |
| 1442 | print ptstr |
| 1443 | print varstr |
| 1444 | |
| 1445 | |
| 1446 | histDict = {} |
| 1447 | histVary = [] |
| 1448 | controlDict = {} |
| 1449 | histoList = Histograms.keys() |
| 1450 | histoList.sort() |
| 1451 | for histogram in histoList: |
| 1452 | Histogram = Histograms[histogram] |
| 1453 | hId = Histogram['hId'] |
| 1454 | pfx = ':'+str(hId)+':' |
| 1455 | controlDict[pfx+'Limits'] = Histogram['Limits'][1] |
| 1456 | |
| 1457 | Background = Histogram['Background'] |
| 1458 | Type,bakDict,bakVary = GetBackgroundParms(hId,Background) |
| 1459 | controlDict[pfx+'bakType'] = Type |
| 1460 | histDict.update(bakDict) |
| 1461 | histVary += bakVary |
| 1462 | |
| 1463 | Inst = Histogram['Instrument Parameters'] |
| 1464 | Type,instDict,insVary = GetInstParms(hId,Inst) |
| 1465 | controlDict[pfx+'histType'] = Type |
| 1466 | if pfx+'Lam1' in instDict: |
| 1467 | controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1'] |
| 1468 | else: |
| 1469 | controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam'] |
| 1470 | histDict.update(instDict) |
| 1471 | histVary += insVary |
| 1472 | |
| 1473 | Sample = Histogram['Sample Parameters'] |
| 1474 | Type,sampDict,sampVary = GetSampleParms(hId,Sample) |
| 1475 | controlDict[pfx+'instType'] = Type |
| 1476 | histDict.update(sampDict) |
| 1477 | histVary += sampVary |
| 1478 | |
| 1479 | if Print: |
| 1480 | print '\n Histogram: ',histogram,' histogram Id: ',hId |
| 1481 | print 135*'-' |
| 1482 | Units = {'C':' deg','T':' msec'} |
| 1483 | units = Units[controlDict[pfx+'histType'][2]] |
| 1484 | Limits = controlDict[pfx+'Limits'] |
| 1485 | print ' Instrument type: ',Sample['Type'] |
| 1486 | print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units) |
| 1487 | PrintSampleParms(Sample) |
| 1488 | PrintInstParms(Inst) |
| 1489 | PrintBackground(Background) |
| 1490 | |
| 1491 | return histVary,histDict,controlDict |
| 1492 | |
| 1493 | def SetHistogramData(parmDict,sigDict,Histograms,Print=True): |
| 1494 | |
| 1495 | def SetBackgroundParms(pfx,Background,parmDict,sigDict): |
| 1496 | Back = Background[0] |
| 1497 | Debye = Background[1] |
| 1498 | lenBack = len(Back[3:]) |
| 1499 | backSig = [0 for i in range(lenBack+3*Debye['nDebye'])] |
| 1500 | for i in range(lenBack): |
| 1501 | Back[3+i] = parmDict[pfx+'Back:'+str(i)] |
| 1502 | if pfx+'Back:'+str(i) in sigDict: |
| 1503 | backSig[i] = sigDict[pfx+'Back:'+str(i)] |
| 1504 | if Debye['nDebye']: |
| 1505 | for i in range(Debye['nDebye']): |
| 1506 | names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)] |
| 1507 | for j,name in enumerate(names): |
| 1508 | Debye['debyeTerms'][i][2*j] = parmDict[name] |
| 1509 | if name in sigDict: |
| 1510 | backSig[lenBack+3*i+j] = sigDict[name] |
| 1511 | return backSig |
| 1512 | |
| 1513 | def SetInstParms(pfx,Inst,parmDict,sigDict): |
| 1514 | insVals,insFlags,insNames = Inst[1:4] |
| 1515 | instSig = [0 for i in range(len(insVals))] |
| 1516 | for i,flag in enumerate(insFlags): |
| 1517 | insName = pfx+insNames[i] |
| 1518 | insVals[i] = parmDict[insName] |
| 1519 | if insName in sigDict: |
| 1520 | instSig[i] = sigDict[insName] |
| 1521 | return instSig |
| 1522 | |
| 1523 | def SetSampleParms(pfx,Sample,parmDict,sigDict): |
| 1524 | if 'Bragg' in Sample['Type']: #Bragg-Brentano |
| 1525 | sampSig = [0 for i in range(3)] |
| 1526 | for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering? |
| 1527 | Sample[item][0] = parmDict[pfx+item] |
| 1528 | if pfx+item in sigDict: |
| 1529 | sampSig[i] = sigDict[pfx+item] |
| 1530 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
| 1531 | sampSig = [0 for i in range(4)] |
| 1532 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
| 1533 | Sample[item][0] = parmDict[pfx+item] |
| 1534 | if pfx+item in sigDict: |
| 1535 | sampSig[i] = sigDict[pfx+item] |
| 1536 | return sampSig |
| 1537 | |
| 1538 | def PrintBackgroundSig(Background,backSig): |
| 1539 | Back = Background[0] |
| 1540 | Debye = Background[1] |
| 1541 | lenBack = len(Back[3:]) |
| 1542 | valstr = ' value : ' |
| 1543 | sigstr = ' sig : ' |
| 1544 | refine = False |
| 1545 | for i,back in enumerate(Back[3:]): |
| 1546 | valstr += '%10.4g'%(back) |
| 1547 | if Back[1]: |
| 1548 | refine = True |
| 1549 | sigstr += '%10.4g'%(backSig[i]) |
| 1550 | else: |
| 1551 | sigstr += 10*' ' |
| 1552 | if refine: |
| 1553 | print '\n Background function: ',Back[0] |
| 1554 | print valstr |
| 1555 | print sigstr |
| 1556 | if Debye['nDebye']: |
| 1557 | ifAny = False |
| 1558 | ptfmt = "%12.5f" |
| 1559 | names = ' names :' |
| 1560 | ptstr = ' values:' |
| 1561 | sigstr = ' esds :' |
| 1562 | for item in sigDict: |
| 1563 | if 'Debye' in item: |
| 1564 | ifAny = True |
| 1565 | names += '%12s'%(item) |
| 1566 | ptstr += ptfmt%(parmDict[item]) |
| 1567 | sigstr += ptfmt%(sigDict[item]) |
| 1568 | if ifAny: |
| 1569 | print '\n Debye diffuse scattering coefficients' |
| 1570 | print names |
| 1571 | print ptstr |
| 1572 | print sigstr |
| 1573 | |
| 1574 | def PrintInstParmsSig(Inst,instSig): |
| 1575 | ptlbls = ' names :' |
| 1576 | ptstr = ' value :' |
| 1577 | sigstr = ' sig :' |
| 1578 | instNames = Inst[3][1:] |
| 1579 | refine = False |
| 1580 | for i,name in enumerate(instNames): |
| 1581 | ptlbls += '%12s' % (name) |
| 1582 | ptstr += '%12.6f' % (Inst[1][i+1]) |
| 1583 | if instSig[i+1]: |
| 1584 | refine = True |
| 1585 | sigstr += '%12.6f' % (instSig[i+1]) |
| 1586 | else: |
| 1587 | sigstr += 12*' ' |
| 1588 | if refine: |
| 1589 | print '\n Instrument Parameters:' |
| 1590 | print ptlbls |
| 1591 | print ptstr |
| 1592 | print sigstr |
| 1593 | |
| 1594 | def PrintSampleParmsSig(Sample,sampleSig): |
| 1595 | ptlbls = ' names :' |
| 1596 | ptstr = ' values:' |
| 1597 | sigstr = ' sig :' |
| 1598 | refine = False |
| 1599 | if 'Bragg' in Sample['Type']: |
| 1600 | for i,item in enumerate(['Scale','Shift','Transparency']): |
| 1601 | ptlbls += '%14s'%(item) |
| 1602 | ptstr += '%14.4f'%(Sample[item][0]) |
| 1603 | if sampleSig[i]: |
| 1604 | refine = True |
| 1605 | sigstr += '%14.4f'%(sampleSig[i]) |
| 1606 | else: |
| 1607 | sigstr += 14*' ' |
| 1608 | |
| 1609 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
| 1610 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
| 1611 | ptlbls += '%14s'%(item) |
| 1612 | ptstr += '%14.4f'%(Sample[item][0]) |
| 1613 | if sampleSig[i]: |
| 1614 | refine = True |
| 1615 | sigstr += '%14.4f'%(sampleSig[i]) |
| 1616 | else: |
| 1617 | sigstr += 14*' ' |
| 1618 | |
| 1619 | if refine: |
| 1620 | print '\n Sample Parameters:' |
| 1621 | print ptlbls |
| 1622 | print ptstr |
| 1623 | print sigstr |
| 1624 | |
| 1625 | histoList = Histograms.keys() |
| 1626 | histoList.sort() |
| 1627 | for histogram in histoList: |
| 1628 | if 'PWDR' in histogram: |
| 1629 | Histogram = Histograms[histogram] |
| 1630 | hId = Histogram['hId'] |
| 1631 | pfx = ':'+str(hId)+':' |
| 1632 | Background = Histogram['Background'] |
| 1633 | backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict) |
| 1634 | |
| 1635 | Inst = Histogram['Instrument Parameters'] |
| 1636 | instSig = SetInstParms(pfx,Inst,parmDict,sigDict) |
| 1637 | |
| 1638 | Sample = Histogram['Sample Parameters'] |
| 1639 | sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict) |
| 1640 | |
| 1641 | print '\n Histogram: ',histogram,' histogram Id: ',hId |
| 1642 | print 135*'-' |
| 1643 | print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs']) |
| 1644 | if Print: |
| 1645 | print ' Instrument type: ',Sample['Type'] |
| 1646 | PrintSampleParmsSig(Sample,sampSig) |
| 1647 | PrintInstParmsSig(Inst,instSig) |
| 1648 | PrintBackgroundSig(Background,backSig) |
| 1649 | |
| 1650 | def GetAtomFXU(pfx,calcControls,parmDict): |
| 1651 | Natoms = calcControls['Natoms'][pfx] |
| 1652 | Tdata = Natoms*[' ',] |
| 1653 | Mdata = np.zeros(Natoms) |
| 1654 | IAdata = Natoms*[' ',] |
| 1655 | Fdata = np.zeros(Natoms) |
| 1656 | FFdata = [] |
| 1657 | BLdata = [] |
| 1658 | Xdata = np.zeros((3,Natoms)) |
| 1659 | dXdata = np.zeros((3,Natoms)) |
| 1660 | Uisodata = np.zeros(Natoms) |
| 1661 | Uijdata = np.zeros((6,Natoms)) |
| 1662 | keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata, |
| 1663 | 'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2], |
| 1664 | 'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata, |
| 1665 | 'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2], |
| 1666 | 'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]} |
| 1667 | for iatm in range(Natoms): |
| 1668 | for key in keys: |
| 1669 | parm = pfx+key+str(iatm) |
| 1670 | if parm in parmDict: |
| 1671 | keys[key][iatm] = parmDict[parm] |
| 1672 | return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata |
| 1673 | |
| 1674 | def getFFvalues(FFtables,SQ): |
| 1675 | FFvals = {} |
| 1676 | for El in FFtables: |
| 1677 | FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0] |
| 1678 | return FFvals |
| 1679 | |
| 1680 | def getBLvalues(BLtables): |
| 1681 | BLvals = {} |
| 1682 | for El in BLtables: |
| 1683 | BLvals[El] = BLtables[El][1] |
| 1684 | return BLvals |
| 1685 | |
| 1686 | def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict): |
| 1687 | ''' Compute structure factors for all h,k,l for phase |
520 | | |
521 | | |
522 | | |
523 | | GPXback = GPXBackup(GPXfile,makeBack) |
524 | | |
525 | | print '\n',135*'-' |
526 | | |
527 | | print 'Read from file:',GPXback |
528 | | |
529 | | print 'Save to file :',GPXfile |
530 | | |
531 | | infile = open(GPXback,'rb') |
532 | | |
533 | | outfile = open(GPXfile,'wb') |
| 1863 | newAtomDict = {} |
| 1864 | for item in parmDict: |
| 1865 | if 'dA' in item: |
| 1866 | parm = ''.join(item.split('d')) |
| 1867 | parmDict[parm] += parmDict[item] |
| 1868 | newAtomDict[item] = [parm,parmDict[parm]] |
| 1869 | return newAtomDict |
| 1870 | |
| 1871 | def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict): |
| 1872 | IFCoup = 'Bragg' in calcControls[hfx+'instType'] |
| 1873 | odfCor = 1.0 |
| 1874 | H = refl[:3] |
| 1875 | cell = G2lat.Gmat2cell(g) |
| 1876 | Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']] |
| 1877 | Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] |
| 1878 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 1879 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs. |
| 1880 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) |
| 1881 | for item in SHnames: |
| 1882 | L,M,N = eval(item.strip('C')) |
| 1883 | Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta) |
| 1884 | Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam) |
| 1885 | Lnorm = G2lat.Lnorm(L) |
| 1886 | odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl |
| 1887 | return odfCor |
| 1888 | |
| 1889 | def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict): |
| 1890 | FORPI = 12.5663706143592 |
| 1891 | IFCoup = 'Bragg' in calcControls[hfx+'instType'] |
| 1892 | odfCor = 1.0 |
| 1893 | dFdODF = {} |
| 1894 | dFdSA = [0,0,0] |
| 1895 | H = refl[:3] |
| 1896 | cell = G2lat.Gmat2cell(g) |
| 1897 | Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']] |
| 1898 | Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] |
| 1899 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 1900 | psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) |
| 1901 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) |
| 1902 | for item in SHnames: |
| 1903 | L,M,N = eval(item.strip('C')) |
| 1904 | Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta) |
| 1905 | Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam) |
| 1906 | Lnorm = G2lat.Lnorm(L) |
| 1907 | odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl |
| 1908 | dFdODF[pfx+item] = Lnorm*Kcl*Ksl |
| 1909 | for i in range(3): |
| 1910 | dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i]) |
| 1911 | return odfCor,dFdODF,dFdSA |
| 1912 | |
| 1913 | def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict): |
| 1914 | odfCor = 1.0 |
| 1915 | H = refl[:3] |
| 1916 | cell = G2lat.Gmat2cell(g) |
| 1917 | Sangl = [0.,0.,0.] |
| 1918 | if 'Bragg' in calcControls[hfx+'instType']: |
| 1919 | Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']] |
| 1920 | IFCoup = True |
| 1921 | else: |
| 1922 | Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']] |
| 1923 | IFCoup = False |
| 1924 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 1925 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. |
| 1926 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) |
| 1927 | for item in SHnames: |
| 1928 | L,N = eval(item.strip('C')) |
| 1929 | Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) |
| 1930 | odfCor += parmDict[phfx+item]*Lnorm*Kcsl |
| 1931 | return odfCor |
| 1932 | |
| 1933 | def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict): |
| 1934 | FORPI = 12.5663706143592 |
| 1935 | odfCor = 1.0 |
| 1936 | dFdODF = {} |
| 1937 | H = refl[:3] |
| 1938 | cell = G2lat.Gmat2cell(g) |
| 1939 | Sangl = [0.,0.,0.] |
| 1940 | if 'Bragg' in calcControls[hfx+'instType']: |
| 1941 | Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']] |
| 1942 | IFCoup = True |
| 1943 | else: |
| 1944 | Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']] |
| 1945 | IFCoup = False |
| 1946 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 1947 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. |
| 1948 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) |
| 1949 | for item in SHnames: |
| 1950 | L,N = eval(item.strip('C')) |
| 1951 | Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) |
| 1952 | odfCor += parmDict[phfx+item]*Lnorm*Kcsl |
| 1953 | dFdODF[phfx+item] = Kcsl*Lnorm |
| 1954 | return odfCor,dFdODF |
| 1955 | |
| 1956 | def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): |
| 1957 | POcorr = 1.0 |
| 1958 | if calcControls[phfx+'poType'] == 'MD': |
| 1959 | MD = parmDict[phfx+'MD'] |
| 1960 | if MD != 1.0: |
| 1961 | MDAxis = calcControls[phfx+'MDAxis'] |
| 1962 | sumMD = 0 |
| 1963 | for H in refl[11]: |
| 1964 | cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) |
| 1965 | A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) |
| 1966 | sumMD += A**3 |
| 1967 | POcorr = sumMD/len(refl[11]) |
| 1968 | else: #spherical harmonics |
| 1969 | if calcControls[pfx+'SHord']: |
| 1970 | POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict) |
| 1971 | return POcorr |
| 1972 | |
| 1973 | def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): |
| 1974 | POcorr = 1.0 |
| 1975 | POderv = {} |
| 1976 | if calcControls[phfx+'poType'] == 'MD': |
| 1977 | MD = parmDict[phfx+'MD'] |
| 1978 | MDAxis = calcControls[phfx+'MDAxis'] |
| 1979 | sumMD = 0 |
| 1980 | sumdMD = 0 |
| 1981 | for H in refl[11]: |
| 1982 | cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) |
| 1983 | A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) |
| 1984 | sumMD += A**3 |
| 1985 | sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2) |
| 1986 | POcorr = sumMD/len(refl[11]) |
| 1987 | POderv[phfx+'MD'] = sumdMD/len(refl[11]) |
| 1988 | else: #spherical harmonics |
| 1989 | if calcControls[pfx+'SHord']: |
| 1990 | POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict) |
| 1991 | return POcorr,POderv |
| 1992 | |
| 1993 | def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): |
| 1994 | Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity |
| 1995 | if 'X' in parmDict[hfx+'Type']: |
| 1996 | Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0] |
| 1997 | Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) |
| 1998 | if pfx+'SHorder' in parmDict: |
| 1999 | Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict) |
| 2000 | refl[13] = Icorr |
| 2001 | |
| 2002 | def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): |
| 2003 | dIdsh = 1./parmDict[hfx+'Scale'] |
| 2004 | dIdsp = 1./parmDict[phfx+'Scale'] |
| 2005 | if 'X' in parmDict[hfx+'Type']: |
| 2006 | pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth']) |
| 2007 | dIdPola /= pola |
| 2008 | else: #'N' |
| 2009 | dIdPola = 0.0 |
| 2010 | POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) |
| 2011 | for iPO in dIdPO: |
| 2012 | dIdPO[iPO] /= POcorr |
| 2013 | dFdODF = {} |
| 2014 | dFdSA = [0,0,0] |
| 2015 | if pfx+'SHorder' in parmDict: |
| 2016 | odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict) |
| 2017 | for iSH in dFdODF: |
| 2018 | dFdODF[iSH] /= odfCor |
| 2019 | for i in range(3): |
| 2020 | dFdSA[i] /= odfCor |
| 2021 | return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA |
| 2022 | |
| 2023 | def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict): |
| 2024 | costh = cosd(refl[5]/2.) |
| 2025 | #crystallite size |
| 2026 | if calcControls[phfx+'SizeType'] == 'isotropic': |
| 2027 | gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh) |
| 2028 | elif calcControls[phfx+'SizeType'] == 'uniaxial': |
| 2029 | H = np.array(refl[:3]) |
| 2030 | P = np.array(calcControls[phfx+'SizeAxis']) |
| 2031 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 2032 | gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh) |
| 2033 | gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2) |
| 2034 | else: #ellipsoidal crystallites |
| 2035 | Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] |
| 2036 | H = np.array(refl[:3]) |
| 2037 | lenR = G2pwd.ellipseSize(H,Sij,GB) |
| 2038 | gam = 1.8*wave/(np.pi*costh*lenR) |
| 2039 | #microstrain |
| 2040 | if calcControls[phfx+'MustrainType'] == 'isotropic': |
| 2041 | gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi |
| 2042 | elif calcControls[phfx+'MustrainType'] == 'uniaxial': |
| 2043 | H = np.array(refl[:3]) |
| 2044 | P = np.array(calcControls[phfx+'MustrainAxis']) |
| 2045 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 2046 | Si = parmDict[phfx+'Mustrain:i'] |
| 2047 | Sa = parmDict[phfx+'Mustrain:a'] |
| 2048 | gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) |
| 2049 | else: #generalized - P.W. Stephens model |
| 2050 | pwrs = calcControls[phfx+'MuPwrs'] |
| 2051 | sum = 0 |
| 2052 | for i,pwr in enumerate(pwrs): |
| 2053 | sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] |
| 2054 | gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum |
| 2055 | return gam |
| 2056 | |
| 2057 | def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict): |
| 2058 | gamDict = {} |
| 2059 | costh = cosd(refl[5]/2.) |
| 2060 | tanth = tand(refl[5]/2.) |
| 2061 | #crystallite size derivatives |
| 2062 | if calcControls[phfx+'SizeType'] == 'isotropic': |
| 2063 | gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh) |
| 2064 | elif calcControls[phfx+'SizeType'] == 'uniaxial': |
| 2065 | H = np.array(refl[:3]) |
| 2066 | P = np.array(calcControls[phfx+'SizeAxis']) |
| 2067 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 2068 | Si = parmDict[phfx+'Size:i'] |
| 2069 | Sa = parmDict[phfx+'Size:a'] |
| 2070 | gami = (1.8*wave/np.pi)/(Si*Sa) |
| 2071 | sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2) |
| 2072 | gam = gami*sqtrm/costh |
| 2073 | gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si |
| 2074 | gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa |
| 2075 | else: #ellipsoidal crystallites |
| 2076 | const = 1.8*wave/(np.pi*costh) |
| 2077 | Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] |
| 2078 | H = np.array(refl[:3]) |
| 2079 | R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) |
| 2080 | for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): |
| 2081 | gamDict[item] = -(const/R**2)*dRdS[i] |
| 2082 | #microstrain derivatives |
| 2083 | if calcControls[phfx+'MustrainType'] == 'isotropic': |
| 2084 | gamDict[phfx+'Mustrain:i'] = 0.018*tanth/np.pi |
| 2085 | elif calcControls[phfx+'MustrainType'] == 'uniaxial': |
| 2086 | H = np.array(refl[:3]) |
| 2087 | P = np.array(calcControls[phfx+'MustrainAxis']) |
| 2088 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 2089 | Si = parmDict[phfx+'Mustrain:i'] |
| 2090 | Sa = parmDict[phfx+'Mustrain:a'] |
| 2091 | gami = 0.018*Si*Sa*tanth/np.pi |
| 2092 | sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) |
| 2093 | gam = gami/sqtrm |
| 2094 | gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3 |
| 2095 | gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3 |
| 2096 | else: #generalized - P.W. Stephens model |
| 2097 | pwrs = calcControls[phfx+'MuPwrs'] |
| 2098 | const = 0.018*refl[4]**2*tanth |
| 2099 | for i,pwr in enumerate(pwrs): |
| 2100 | gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] |
| 2101 | return gamDict |
| 2102 | |
| 2103 | def GetReflPos(refl,wave,G,hfx,calcControls,parmDict): |
| 2104 | h,k,l = refl[:3] |
| 2105 | dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G) |
| 2106 | d = np.sqrt(dsq) |
| 2107 | |
| 2108 | refl[4] = d |
| 2109 | pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero'] |
| 2110 | const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns |
| 2111 | if 'Bragg' in calcControls[hfx+'instType']: |
| 2112 | pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \ |
| 2113 | parmDict[hfx+'Transparency']*sind(pos)*100.0) #trans(=1/mueff) in cm |
| 2114 | else: #Debye-Scherrer - simple but maybe not right |
| 2115 | pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos)) |
| 2116 | return pos |
| 2117 | |
| 2118 | def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict): |
| 2119 | dpr = 180./np.pi |
| 2120 | h,k,l = refl[:3] |
| 2121 | dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A) |
| 2122 | dst = np.sqrt(dstsq) |
| 2123 | pos = refl[5]-parmDict[hfx+'Zero'] |
| 2124 | const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0) |
| 2125 | dpdw = const*dst |
| 2126 | dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) |
| 2127 | dpdA *= const*wave/(2.0*dst) |
| 2128 | dpdZ = 1.0 |
| 2129 | const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns |
| 2130 | if 'Bragg' in calcControls[hfx+'instType']: |
| 2131 | dpdSh = -4.*const*cosd(pos/2.0) |
| 2132 | dpdTr = -const*sind(pos)*100.0 |
| 2133 | return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0. |
| 2134 | else: #Debye-Scherrer - simple but maybe not right |
| 2135 | dpdXd = -const*cosd(pos) |
| 2136 | dpdYd = -const*sind(pos) |
| 2137 | return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd |
| 2138 | |
| 2139 | def GetHStrainShift(refl,SGData,phfx,parmDict): |
| 2140 | laue = SGData['SGLaue'] |
| 2141 | uniq = SGData['SGUniq'] |
| 2142 | h,k,l = refl[:3] |
| 2143 | if laue in ['m3','m3m']: |
| 2144 | Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \ |
| 2145 | refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2 |
| 2146 | elif laue in ['6/m','6/mmm','3m1','31m','3']: |
| 2147 | Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2 |
| 2148 | elif laue in ['3R','3mR']: |
| 2149 | Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l) |
| 2150 | elif laue in ['4/m','4/mmm']: |
| 2151 | Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2 |
| 2152 | elif laue in ['mmm']: |
| 2153 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2 |
| 2154 | elif laue in ['2/m']: |
| 2155 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2 |
| 2156 | if uniq == 'a': |
| 2157 | Dij += parmDict[phfx+'D23']*k*l |
| 2158 | elif uniq == 'b': |
| 2159 | Dij += parmDict[phfx+'D13']*h*l |
| 2160 | elif uniq == 'c': |
| 2161 | Dij += parmDict[phfx+'D12']*h*k |
| 2162 | else: |
| 2163 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \ |
| 2164 | parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l |
| 2165 | return Dij*refl[4]**2*tand(refl[5]/2.0) |
| 2166 | |
| 2167 | def GetHStrainShiftDerv(refl,SGData,phfx): |
| 2168 | laue = SGData['SGLaue'] |
| 2169 | uniq = SGData['SGUniq'] |
| 2170 | h,k,l = refl[:3] |
| 2171 | if laue in ['m3','m3m']: |
| 2172 | dDijDict = {phfx+'D11':h**2+k**2+l**2, |
| 2173 | phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2} |
| 2174 | elif laue in ['6/m','6/mmm','3m1','31m','3']: |
| 2175 | dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2} |
| 2176 | elif laue in ['3R','3mR']: |
| 2177 | dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l} |
| 2178 | elif laue in ['4/m','4/mmm']: |
| 2179 | dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2} |
| 2180 | elif laue in ['mmm']: |
| 2181 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2} |
| 2182 | elif laue in ['2/m']: |
| 2183 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2} |
| 2184 | if uniq == 'a': |
| 2185 | dDijDict[phfx+'D23'] = k*l |
| 2186 | elif uniq == 'b': |
| 2187 | dDijDict[phfx+'D13'] = h*l |
| 2188 | elif uniq == 'c': |
| 2189 | dDijDict[phfx+'D12'] = h*k |
| 2190 | names.append() |
| 2191 | else: |
| 2192 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2, |
| 2193 | phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l} |
| 2194 | for item in dDijDict: |
| 2195 | dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0) |
| 2196 | return dDijDict |
| 2197 | |
| 2198 | def GetFprime(controlDict,Histograms): |
| 2199 | FFtables = controlDict['FFtables'] |
| 2200 | if not FFtables: |
| 2201 | return |
| 2202 | histoList = Histograms.keys() |
| 2203 | histoList.sort() |
| 2204 | for histogram in histoList: |
| 2205 | if 'PWDR' in histogram[:4]: |
| 2206 | Histogram = Histograms[histogram] |
| 2207 | hId = Histogram['hId'] |
| 2208 | hfx = ':%d:'%(hId) |
| 2209 | keV = controlDict[hfx+'keV'] |
| 2210 | for El in FFtables: |
| 2211 | Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0]) |
| 2212 | FP,FPP,Mu = G2el.FPcalc(Orbs, keV) |
| 2213 | FFtables[El][hfx+'FP'] = FP |
| 2214 | FFtables[El][hfx+'FPP'] = FPP |
| 2215 | |
| 2216 | def GetFobsSq(Histograms,Phases,parmDict,calcControls): |
| 2217 | histoList = Histograms.keys() |
| 2218 | histoList.sort() |
| 2219 | for histogram in histoList: |
| 2220 | if 'PWDR' in histogram[:4]: |
| 2221 | Histogram = Histograms[histogram] |
| 2222 | hId = Histogram['hId'] |
| 2223 | hfx = ':%d:'%(hId) |
| 2224 | Limits = calcControls[hfx+'Limits'] |
| 2225 | shl = max(parmDict[hfx+'SH/L'],0.002) |
| 2226 | Ka2 = False |
| 2227 | kRatio = 0.0 |
| 2228 | if hfx+'Lam1' in parmDict.keys(): |
| 2229 | Ka2 = True |
| 2230 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
| 2231 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
| 2232 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 2233 | ymb = np.array(y-yb) |
| 2234 | ycmb = np.array(yc-yb) |
| 2235 | ratio = np.where(ycmb!=0.,ymb/ycmb,0.0) |
| 2236 | xB = np.searchsorted(x,Limits[0]) |
| 2237 | xF = np.searchsorted(x,Limits[1]) |
| 2238 | refLists = Histogram['Reflection Lists'] |
| 2239 | for phase in refLists: |
| 2240 | Phase = Phases[phase] |
| 2241 | pId = Phase['pId'] |
| 2242 | phfx = '%d:%d:'%(pId,hId) |
| 2243 | refList = refLists[phase] |
| 2244 | sumFo = 0.0 |
| 2245 | sumdF = 0.0 |
| 2246 | sumFosq = 0.0 |
| 2247 | sumdFsq = 0.0 |
| 2248 | for refl in refList: |
| 2249 | if 'C' in calcControls[hfx+'histType']: |
| 2250 | yp = np.zeros_like(yb) |
| 2251 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
| 2252 | iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin) |
| 2253 | iFin = np.searchsorted(x[xB:xF],refl[5]+fmax) |
| 2254 | iFin2 = iFin |
| 2255 | yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here |
| 2256 | if Ka2: |
| 2257 | pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) |
| 2258 | Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) |
| 2259 | iBeg2 = np.searchsorted(x,pos2-fmin) |
| 2260 | iFin2 = np.searchsorted(x,pos2+fmax) |
| 2261 | yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here |
| 2262 | refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0)) |
| 2263 | elif 'T' in calcControls[hfx+'histType']: |
| 2264 | print 'TOF Undefined at present' |
| 2265 | raise Exception #no TOF yet |
| 2266 | Fo = np.sqrt(np.abs(refl[8])) |
| 2267 | Fc = np.sqrt(np.abs(refl[9])) |
| 2268 | sumFo += Fo |
| 2269 | sumFosq += refl[8]**2 |
| 2270 | sumdF += np.abs(Fo-Fc) |
| 2271 | sumdFsq += (refl[8]-refl[9])**2 |
| 2272 | Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.) |
| 2273 | Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.) |
| 2274 | Histogram[phfx+'Nref'] = len(refList) |
| 2275 | |
| 2276 | def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): |
| 2277 | |
| 2278 | def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): |
| 2279 | U = parmDict[hfx+'U'] |
| 2280 | V = parmDict[hfx+'V'] |
| 2281 | W = parmDict[hfx+'W'] |
| 2282 | X = parmDict[hfx+'X'] |
| 2283 | Y = parmDict[hfx+'Y'] |
| 2284 | tanPos = tand(refl[5]/2.0) |
| 2285 | sig = U*tanPos**2+V*tanPos+W #save peak sigma |
| 2286 | sig = max(0.001,sig) |
| 2287 | gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma |
| 2288 | gam = max(0.001,gam) |
| 2289 | return sig,gam |
| 2290 | |
| 2291 | hId = Histogram['hId'] |
| 2292 | hfx = ':%d:'%(hId) |
| 2293 | bakType = calcControls[hfx+'bakType'] |
| 2294 | yb = G2pwd.getBackground(hfx,parmDict,bakType,x) |
| 2295 | yc = np.zeros_like(yb) |
| 2296 | |
| 2297 | if 'C' in calcControls[hfx+'histType']: |
| 2298 | shl = max(parmDict[hfx+'SH/L'],0.002) |
| 2299 | Ka2 = False |
| 2300 | if hfx+'Lam1' in parmDict.keys(): |
| 2301 | wave = parmDict[hfx+'Lam1'] |
| 2302 | Ka2 = True |
| 2303 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
| 2304 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
| 2305 | else: |
| 2306 | wave = parmDict[hfx+'Lam'] |
| 2307 | else: |
| 2308 | print 'TOF Undefined at present' |
| 2309 | raise ValueError |
| 2310 | for phase in Histogram['Reflection Lists']: |
| 2311 | refList = Histogram['Reflection Lists'][phase] |
| 2312 | Phase = Phases[phase] |
| 2313 | pId = Phase['pId'] |
| 2314 | pfx = '%d::'%(pId) |
| 2315 | phfx = '%d:%d:'%(pId,hId) |
| 2316 | hfx = ':%d:'%(hId) |
| 2317 | SGData = Phase['General']['SGData'] |
| 2318 | A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] |
| 2319 | G,g = G2lat.A2Gmat(A) #recip & real metric tensors |
| 2320 | GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies |
| 2321 | Vst = np.sqrt(nl.det(G)) #V* |
| 2322 | if 'Pawley' not in Phase['General']['Type']: |
| 2323 | refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict) |
| 2324 | for refl in refList: |
| 2325 | if 'C' in calcControls[hfx+'histType']: |
| 2326 | h,k,l = refl[:3] |
| 2327 | refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position |
| 2328 | Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.)) #Lorentz correction |
| 2329 | refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift |
| 2330 | refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam |
| 2331 | GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13] |
| 2332 | refl[13] *= Vst*Lorenz |
| 2333 | if 'Pawley' in Phase['General']['Type']: |
| 2334 | try: |
| 2335 | refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) |
| 2336 | except KeyError: |
| 2337 | # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) |
| 2338 | continue |
| 2339 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
| 2340 | iBeg = np.searchsorted(x,refl[5]-fmin) |
| 2341 | iFin = np.searchsorted(x,refl[5]+fmax) |
| 2342 | if not iBeg+iFin: #peak below low limit - skip peak |
| 2343 | continue |
| 2344 | elif not iBeg-iFin: #peak above high limit - done |
| 2345 | break |
| 2346 | yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here |
| 2347 | if Ka2: |
| 2348 | pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) |
| 2349 | Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) |
| 2350 | iBeg = np.searchsorted(x,pos2-fmin) |
| 2351 | iFin = np.searchsorted(x,pos2+fmax) |
| 2352 | if not iBeg+iFin: #peak below low limit - skip peak |
| 2353 | continue |
| 2354 | elif not iBeg-iFin: #peak above high limit - done |
| 2355 | return yc,yb |
| 2356 | yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here |
| 2357 | elif 'T' in calcControls[hfx+'histType']: |
| 2358 | print 'TOF Undefined at present' |
| 2359 | raise Exception #no TOF yet |
| 2360 | return yc,yb |
| 2361 | |
| 2362 | def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): |
| 2363 | |
| 2364 | def cellVaryDerv(pfx,SGData,dpdA): |
| 2365 | if SGData['SGLaue'] in ['-1',]: |
| 2366 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]], |
| 2367 | [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]] |
| 2368 | elif SGData['SGLaue'] in ['2/m',]: |
| 2369 | if SGData['SGUniq'] == 'a': |
| 2370 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]] |
| 2371 | elif SGData['SGUniq'] == 'b': |
| 2372 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]] |
| 2373 | else: |
| 2374 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]] |
| 2375 | elif SGData['SGLaue'] in ['mmm',]: |
| 2376 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]] |
| 2377 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
| 2378 | # return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]] |
| 2379 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
| 2380 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
| 2381 | # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]] |
| 2382 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
| 2383 | elif SGData['SGLaue'] in ['3R', '3mR']: |
| 2384 | return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]] |
| 2385 | elif SGData['SGLaue'] in ['m3m','m3']: |
| 2386 | # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]] |
| 2387 | return [[pfx+'A0',dpdA[0]]] |
| 2388 | # create a list of dependent variables and set up a dictionary to hold their derivatives |
| 2389 | dependentVars = G2mv.GetDependentVars() |
| 2390 | depDerivDict = {} |
| 2391 | for j in dependentVars: |
| 2392 | depDerivDict[j] = np.zeros(shape=(len(x))) |
| 2393 | #print 'dependent vars',dependentVars |
| 2394 | lenX = len(x) |
| 2395 | hId = Histogram['hId'] |
| 2396 | hfx = ':%d:'%(hId) |
| 2397 | bakType = calcControls[hfx+'bakType'] |
| 2398 | dMdv = np.zeros(shape=(len(varylist),len(x))) |
| 2399 | dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x) |
| 2400 | if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints |
| 2401 | bBpos =varylist.index(hfx+'Back:0') |
| 2402 | dMdv[bBpos:bBpos+len(dMdb)] = dMdb |
| 2403 | names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU'] |
| 2404 | for name in varylist: |
| 2405 | if 'Debye' in name: |
| 2406 | id = int(name.split(':')[-1]) |
| 2407 | parm = name[:int(name.rindex(':'))] |
| 2408 | ip = names.index(parm) |
| 2409 | dMdv[varylist.index(name)] = dMddb[3*id+ip] |
| 2410 | if 'C' in calcControls[hfx+'histType']: |
| 2411 | dx = x[1]-x[0] |
| 2412 | shl = max(parmDict[hfx+'SH/L'],0.002) |
| 2413 | Ka2 = False |
| 2414 | if hfx+'Lam1' in parmDict.keys(): |
| 2415 | wave = parmDict[hfx+'Lam1'] |
| 2416 | Ka2 = True |
| 2417 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
| 2418 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
| 2419 | else: |
| 2420 | wave = parmDict[hfx+'Lam'] |
| 2421 | else: |
| 2422 | print 'TOF Undefined at present' |
| 2423 | raise ValueError |
| 2424 | for phase in Histogram['Reflection Lists']: |
| 2425 | refList = Histogram['Reflection Lists'][phase] |
| 2426 | Phase = Phases[phase] |
| 2427 | SGData = Phase['General']['SGData'] |
| 2428 | pId = Phase['pId'] |
| 2429 | pfx = '%d::'%(pId) |
| 2430 | phfx = '%d:%d:'%(pId,hId) |
| 2431 | A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] |
| 2432 | G,g = G2lat.A2Gmat(A) #recip & real metric tensors |
| 2433 | GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies |
| 2434 | if 'Pawley' not in Phase['General']['Type']: |
| 2435 | dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict) |
| 2436 | for iref,refl in enumerate(refList): |
| 2437 | if 'C' in calcControls[hfx+'histType']: #CW powder |
| 2438 | h,k,l = refl[:3] |
| 2439 | dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) |
| 2440 | if 'Pawley' in Phase['General']['Type']: |
| 2441 | try: |
| 2442 | refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) |
| 2443 | except KeyError: |
| 2444 | # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) |
| 2445 | continue |
| 2446 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
| 2447 | iBeg = np.searchsorted(x,refl[5]-fmin) |
| 2448 | iFin = np.searchsorted(x,refl[5]+fmax) |
| 2449 | if not iBeg+iFin: #peak below low limit - skip peak |
| 2450 | continue |
| 2451 | elif not iBeg-iFin: #peak above high limit - done |
| 2452 | break |
| 2453 | pos = refl[5] |
| 2454 | tanth = tand(pos/2.0) |
| 2455 | costh = cosd(pos/2.0) |
| 2456 | lenBF = iFin-iBeg |
| 2457 | dMdpk = np.zeros(shape=(6,lenBF)) |
| 2458 | dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) |
| 2459 | for i in range(1,5): |
| 2460 | dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i] |
| 2461 | dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0] |
| 2462 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])} |
| 2463 | if Ka2: |
| 2464 | pos2 = refl[5]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) |
| 2465 | kdelt = int((pos2-refl[5])/dx) |
| 2466 | iBeg2 = min(lenX,iBeg+kdelt) |
| 2467 | iFin2 = min(lenX,iFin+kdelt) |
| 2468 | if iBeg2-iFin2: |
| 2469 | lenBF2 = iFin2-iBeg2 |
| 2470 | dMdpk2 = np.zeros(shape=(6,lenBF2)) |
| 2471 | dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) |
| 2472 | for i in range(1,5): |
| 2473 | dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i] |
| 2474 | dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0] |
| 2475 | dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0] |
| 2476 | dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} |
| 2477 | if 'Pawley' in Phase['General']['Type']: |
| 2478 | try: |
| 2479 | idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])) |
| 2480 | dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9] |
| 2481 | if Ka2: |
| 2482 | dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9] |
| 2483 | # Assuming Pawley variables not in constraints |
| 2484 | except ValueError: |
| 2485 | pass |
| 2486 | dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) |
| 2487 | names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], |
| 2488 | hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], |
| 2489 | hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], |
| 2490 | hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], |
| 2491 | hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], |
| 2492 | hfx+'DisplaceY':[dpdY,'pos'],} |
| 2493 | for name in names: |
| 2494 | item = names[name] |
| 2495 | if name in varylist: |
| 2496 | dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]] |
| 2497 | if Ka2: |
| 2498 | dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]] |
| 2499 | elif name in dependentVars: |
| 2500 | if Ka2: |
| 2501 | depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]] |
| 2502 | depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]] |
| 2503 | for iPO in dIdPO: |
| 2504 | if iPO in varylist: |
| 2505 | dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] |
| 2506 | if Ka2: |
| 2507 | dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] |
| 2508 | elif iPO in dependentVars: |
| 2509 | depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int'] |
| 2510 | if Ka2: |
| 2511 | depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int'] |
| 2512 | for i,name in enumerate(['omega','chi','phi']): |
| 2513 | aname = pfx+'SH '+name |
| 2514 | if aname in varylist: |
| 2515 | dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int'] |
| 2516 | if Ka2: |
| 2517 | dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] |
| 2518 | elif aname in dependentVars: |
| 2519 | depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int'] |
| 2520 | if Ka2: |
| 2521 | depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int'] |
| 2522 | for iSH in dFdODF: |
| 2523 | if iSH in varylist: |
| 2524 | dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] |
| 2525 | if Ka2: |
| 2526 | dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] |
| 2527 | elif iSH in dependentVars: |
| 2528 | depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int'] |
| 2529 | if Ka2: |
| 2530 | depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int'] |
| 2531 | cellDervNames = cellVaryDerv(pfx,SGData,dpdA) |
| 2532 | for name,dpdA in cellDervNames: |
| 2533 | if name in varylist: |
| 2534 | dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos'] |
| 2535 | if Ka2: |
| 2536 | dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos'] |
| 2537 | elif name in dependentVars: |
| 2538 | depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos'] |
| 2539 | if Ka2: |
| 2540 | depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos'] |
| 2541 | dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) |
| 2542 | for name in dDijDict: |
| 2543 | if name in varylist: |
| 2544 | dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos'] |
| 2545 | if Ka2: |
| 2546 | dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] |
| 2547 | elif name in dependentVars: |
| 2548 | depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos'] |
| 2549 | if Ka2: |
| 2550 | depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos'] |
| 2551 | gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) |
| 2552 | for name in gamDict: |
| 2553 | if name in varylist: |
| 2554 | dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam'] |
| 2555 | if Ka2: |
| 2556 | dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] |
| 2557 | elif name in dependentVars: |
| 2558 | depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam'] |
| 2559 | if Ka2: |
| 2560 | depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam'] |
| 2561 | |
| 2562 | elif 'T' in calcControls[hfx+'histType']: |
| 2563 | print 'TOF Undefined at present' |
| 2564 | raise Exception #no TOF yet |
| 2565 | #do atom derivatives - for F,X & U so far |
| 2566 | corr = dervDict['int']/refl[9] |
| 2567 | if Ka2: |
| 2568 | corr2 = dervDict2['int']/refl[9] |
| 2569 | for name in varylist+dependentVars: |
| 2570 | try: |
| 2571 | aname = name.split(pfx)[1][:2] |
| 2572 | if aname not in ['Af','dA','AU']: continue # skip anything not an atom param |
| 2573 | except IndexError: |
| 2574 | continue |
| 2575 | if name in varylist: |
| 2576 | dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr |
| 2577 | if Ka2: |
| 2578 | dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 |
| 2579 | elif name in dependentVars: |
| 2580 | depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr |
| 2581 | if Ka2: |
| 2582 | depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2 |
| 2583 | # now process derivatives in constraints |
| 2584 | G2mv.Dict2Deriv(varylist,depDerivDict,dMdv) |
| 2585 | return dMdv |
| 2586 | |
| 2587 | def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 2588 | parmdict.update(zip(varylist,values)) |
| 2589 | G2mv.Dict2Map(parmdict,varylist) |
| 2590 | Histograms,Phases = HistoPhases |
| 2591 | nvar = len(varylist) |
| 2592 | dMdv = np.empty(0) |
| 2593 | histoList = Histograms.keys() |
| 2594 | histoList.sort() |
| 2595 | for histogram in histoList: |
| 2596 | if 'PWDR' in histogram[:4]: |
| 2597 | Histogram = Histograms[histogram] |
| 2598 | hId = Histogram['hId'] |
| 2599 | hfx = ':%d:'%(hId) |
| 2600 | Limits = calcControls[hfx+'Limits'] |
| 2601 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 2602 | xB = np.searchsorted(x,Limits[0]) |
| 2603 | xF = np.searchsorted(x,Limits[1]) |
| 2604 | dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], |
| 2605 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
| 2606 | if len(dMdv): |
| 2607 | dMdv = np.concatenate((dMdv.T,dMdvh.T)).T |
| 2608 | else: |
| 2609 | dMdv = dMdvh |
| 2610 | return dMdv |
| 2611 | |
| 2612 | def ComputePowderHessian(args): |
| 2613 | Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args |
| 2614 | hId = Histogram['hId'] |
| 2615 | hfx = ':%d:'%(hId) |
| 2616 | Limits = calcControls[hfx+'Limits'] |
| 2617 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 2618 | dy = y-yc |
| 2619 | xB = np.searchsorted(x,Limits[0]) |
| 2620 | xF = np.searchsorted(x,Limits[1]) |
| 2621 | dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv( |
| 2622 | parmdict,x[xB:xF], |
| 2623 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
| 2624 | Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) |
| 2625 | Hess = np.inner(dMdvh,dMdvh) |
| 2626 | return Vec,Hess |
| 2627 | |
| 2628 | #def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 2629 | # parmdict.update(zip(varylist,values)) |
| 2630 | # G2mv.Dict2Map(parmdict,varylist) |
| 2631 | # Histograms,Phases = HistoPhases |
| 2632 | # nvar = len(varylist) |
| 2633 | # Hess = np.empty(0) |
| 2634 | # histoList = Histograms.keys() |
| 2635 | # histoList.sort() |
| 2636 | # for histogram in histoList: |
| 2637 | # if 'PWDR' in histogram[:4]: |
| 2638 | # Histogram = Histograms[histogram] |
| 2639 | # hId = Histogram['hId'] |
| 2640 | # hfx = ':%d:'%(hId) |
| 2641 | # Limits = calcControls[hfx+'Limits'] |
| 2642 | # x,y,w,yc,yb,yd = Histogram['Data'] |
| 2643 | # dy = y-yc |
| 2644 | # xB = np.searchsorted(x,Limits[0]) |
| 2645 | # xF = np.searchsorted(x,Limits[1]) |
| 2646 | # dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], |
| 2647 | # varylist,Histogram,Phases,calcControls,pawleyLookup) |
| 2648 | # if dlg: |
| 2649 | # dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0] |
| 2650 | # if len(Hess): |
| 2651 | # Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) |
| 2652 | # Hess += np.inner(dMdvh,dMdvh) |
| 2653 | # else: |
| 2654 | # Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) |
| 2655 | # Hess = np.inner(dMdvh,dMdvh) |
| 2656 | # return Vec,Hess |
| 2657 | |
| 2658 | def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 2659 | parmdict.update(zip(varylist,values)) |
| 2660 | G2mv.Dict2Map(parmdict,varylist) |
| 2661 | Histograms,Phases = HistoPhases |
| 2662 | nvar = len(varylist) |
| 2663 | HessSum = np.zeros((nvar,nvar)) |
| 2664 | VecSum = np.zeros(nvar) |
| 2665 | histoList = Histograms.keys() |
| 2666 | histoList.sort() |
| 2667 | argList = [] |
| 2668 | MaxProcess = calcControls['max Hprocess'] |
| 2669 | for histogram in histoList: |
| 2670 | if 'PWDR' in histogram[:4]: |
| 2671 | Histogram = Histograms[histogram] |
| 2672 | argList.append([ |
| 2673 | Histogram,parmdict,varylist,Phases, |
| 2674 | calcControls,pawleyLookup]) |
| 2675 | if MaxProcess > 1: |
| 2676 | mpPool = mp.Pool(processes=MaxProcess) |
| 2677 | results = mpPool.map(ComputePowderHessian,argList) |
| 2678 | for Vec,Hess in results: |
| 2679 | VecSum += Vec |
| 2680 | HessSum += Hess |
| 2681 | else: |
| 2682 | for arg in argList: |
| 2683 | Vec,Hess = ComputePowderHessian(arg) |
| 2684 | VecSum += Vec |
| 2685 | HessSum += Hess |
| 2686 | return VecSum,HessSum |
| 2687 | |
| 2688 | def ComputePowderProfile(args): |
| 2689 | Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args |
| 2690 | hId = Histogram['hId'] |
| 2691 | hfx = ':%d:'%(hId) |
| 2692 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 2693 | Limits = calcControls[hfx+'Limits'] |
| 2694 | xB = np.searchsorted(x,Limits[0]) |
| 2695 | xF = np.searchsorted(x,Limits[1]) |
| 2696 | yc,yb = getPowderProfile(parmdict,x[xB:xF], |
| 2697 | varylist,Histogram,Phases,calcControls, |
| 2698 | pawleyLookup) |
| 2699 | return xB,xF,yc,yb,Histogram['Reflection Lists'] |
| 2700 | |
| 2701 | #def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 2702 | # parmdict.update(zip(varylist,values)) |
| 2703 | # Values2Dict(parmdict, varylist, values) |
| 2704 | # G2mv.Dict2Map(parmdict,varylist) |
| 2705 | # Histograms,Phases = HistoPhases |
| 2706 | # M = np.empty(0) |
| 2707 | # sumwYo = 0 |
| 2708 | # Nobs = 0 |
| 2709 | # histoList = Histograms.keys() |
| 2710 | # histoList.sort() |
| 2711 | # for histogram in histoList: |
| 2712 | # if 'PWDR' in histogram[:4]: |
| 2713 | # Histogram = Histograms[histogram] |
| 2714 | # hId = Histogram['hId'] |
| 2715 | # hfx = ':%d:'%(hId) |
| 2716 | # Limits = calcControls[hfx+'Limits'] |
| 2717 | # x,y,w,yc,yb,yd = Histogram['Data'] |
| 2718 | # yc *= 0.0 #zero full calcd profiles |
| 2719 | # yb *= 0.0 |
| 2720 | # yd *= 0.0 |
| 2721 | # xB = np.searchsorted(x,Limits[0]) |
| 2722 | # xF = np.searchsorted(x,Limits[1]) |
| 2723 | # Histogram['Nobs'] = xF-xB |
| 2724 | # Nobs += Histogram['Nobs'] |
| 2725 | # Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) |
| 2726 | # sumwYo += Histogram['sumwYo'] |
| 2727 | # yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF], |
| 2728 | # varylist,Histogram,Phases,calcControls,pawleyLookup) |
| 2729 | # yc[xB:xF] += yb[xB:xF] |
| 2730 | # yd[xB:xF] = y[xB:xF]-yc[xB:xF] |
| 2731 | # Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) |
| 2732 | # wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF]) |
| 2733 | # Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) |
| 2734 | # if dlg: |
| 2735 | # dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0] |
| 2736 | # M = np.concatenate((M,wdy)) |
| 2737 | # Histograms['sumwYo'] = sumwYo |
| 2738 | # Histograms['Nobs'] = Nobs |
| 2739 | # Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.) |
| 2740 | # if dlg: |
| 2741 | # GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0] |
| 2742 | # if not GoOn: |
| 2743 | # parmDict['saved values'] = values |
| 2744 | # raise Exception #Abort!! |
| 2745 | # return M |
| 2746 | # |
| 2747 | def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 2748 | parmdict.update(zip(varylist,values)) |
| 2749 | Values2Dict(parmdict, varylist, values) |
| 2750 | G2mv.Dict2Map(parmdict,varylist) |
| 2751 | Histograms,Phases = HistoPhases |
| 2752 | M = np.empty(0) |
| 2753 | sumwYo = 0 |
| 2754 | Nobs = 0 |
| 2755 | histoList = Histograms.keys() |
| 2756 | histoList.sort() |
| 2757 | argList = [] |
| 2758 | MaxProcess = calcControls['max Hprocess'] |
| 2759 | for histogram in histoList: |
| 2760 | if 'PWDR' in histogram[:4]: |
| 2761 | Histogram = Histograms[histogram] |
| 2762 | argList.append( |
| 2763 | [Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup] |
| 2764 | ) |
| 2765 | if MaxProcess > 1: |
| 2766 | mpPool = mp.Pool(processes=MaxProcess) |
| 2767 | results = mpPool.map(ComputePowderProfile,argList) |
| 2768 | for arg,res in zip(argList,results): |
| 2769 | xB,xF,ycSect,ybSect,RL = res |
| 2770 | Histogram = arg[0] |
| 2771 | Histogram['Reflection Lists'] = RL |
| 2772 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 2773 | yc *= 0.0 #zero full calcd profiles |
| 2774 | yb *= 0.0 |
| 2775 | yd *= 0.0 |
| 2776 | Histogram['Nobs'] = xF-xB |
| 2777 | Nobs += Histogram['Nobs'] |
| 2778 | Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) |
| 2779 | sumwYo += Histogram['sumwYo'] |
| 2780 | |
| 2781 | yc[xB:xF] = ycSect |
| 2782 | yb[xB:xF] = ybSect |
| 2783 | yc[xB:xF] += yb[xB:xF] |
| 2784 | yd[xB:xF] = y[xB:xF]-yc[xB:xF] |
| 2785 | Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) |
| 2786 | wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF]) |
| 2787 | Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) |
| 2788 | M = np.concatenate((M,wdy)) |
| 2789 | else: |
| 2790 | for arg in argList: |
| 2791 | xB,xF,ycSect,ybSect,RL = ComputePowderProfile(arg) |
| 2792 | Histogram = arg[0] |
| 2793 | hId = Histogram['hId'] |
| 2794 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 2795 | yc *= 0.0 #zero full calcd profiles |
| 2796 | yb *= 0.0 |
| 2797 | yd *= 0.0 |
| 2798 | Histogram['Nobs'] = xF-xB |
| 2799 | Nobs += Histogram['Nobs'] |
| 2800 | Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) |
| 2801 | sumwYo += Histogram['sumwYo'] |
| 2802 | |
| 2803 | yc[xB:xF] = ycSect |
| 2804 | yb[xB:xF] = ybSect |
| 2805 | yc[xB:xF] += yb[xB:xF] |
| 2806 | yd[xB:xF] = y[xB:xF]-yc[xB:xF] |
| 2807 | Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) |
| 2808 | wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF]) |
| 2809 | Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) |
| 2810 | if dlg: |
| 2811 | dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0] |
| 2812 | M = np.concatenate((M,wdy)) |
| 2813 | |
| 2814 | Histograms['sumwYo'] = sumwYo |
| 2815 | Histograms['Nobs'] = Nobs |
| 2816 | Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.) |
| 2817 | if dlg: |
| 2818 | GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0] |
| 2819 | if not GoOn: |
| 2820 | parmDict['saved values'] = values |
| 2821 | raise Exception #Abort!! |
| 2822 | return M |
| 2823 | |
| 2824 | def Refine(GPXfile,dlg): |
| 2825 | import pytexture as ptx |
| 2826 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
| 2827 | |
| 2828 | ShowBanner() |
| 2829 | varyList = [] |
| 2830 | parmDict = {} |
| 2831 | G2mv.InitVars() |
| 2832 | Controls = GetControls(GPXfile) |
| 2833 | ShowControls(Controls) |
| 2834 | calcControls = {} |
| 2835 | calcControls.update(Controls) |
| 2836 | constrDict,constrFlag,fixedList = GetConstraints(GPXfile) |
| 2837 | Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) |
| 2838 | if not Phases: |
| 2839 | print ' *** ERROR - you have no histograms to refine! ***' |
| 2840 | print ' *** Refine aborted ***' |
| 2841 | raise Exception |
| 2842 | if not Histograms: |
| 2843 | print ' *** ERROR - you have no data to refine with! ***' |
| 2844 | print ' *** Refine aborted ***' |
| 2845 | raise Exception |
| 2846 | Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases) |
| 2847 | calcControls['Natoms'] = Natoms |
| 2848 | calcControls['FFtables'] = FFtables |
| 2849 | calcControls['BLtables'] = BLtables |
| 2850 | hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms) |
| 2851 | calcControls.update(controlDict) |
| 2852 | histVary,histDict,controlDict = GetHistogramData(Histograms) |
| 2853 | calcControls.update(controlDict) |
| 2854 | varyList = phaseVary+hapVary+histVary |
| 2855 | parmDict.update(phaseDict) |
| 2856 | parmDict.update(hapDict) |
| 2857 | parmDict.update(histDict) |
| 2858 | GetFprime(calcControls,Histograms) |
| 2859 | # do constraint processing |
| 2860 | try: |
| 2861 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
| 2862 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList) |
| 2863 | except: |
| 2864 | print ' *** ERROR - your constraints are internally inconsistent ***' |
| 2865 | raise Exception(' *** Refine aborted ***') |
| 2866 | # check to see which generated parameters are fully varied |
| 2867 | msg = G2mv.SetVaryFlags(varyList) |
| 2868 | if msg: |
| 2869 | print ' *** ERROR - you have not set the refine flags for constraints consistently! ***' |
| 2870 | print msg |
| 2871 | raise Exception(' *** Refine aborted ***') |
| 2872 | G2mv.Map2Dict(parmDict,varyList) |
| 2873 | # print G2mv.VarRemapShow(varyList) |
538 | | |
539 | | data = cPickle.load(infile) |
540 | | |
541 | | except EOFError: |
542 | | |
543 | | break |
544 | | |
545 | | datum = data[0] |
546 | | |
547 | | # print 'read: ',datum[0] |
548 | | |
549 | | if datum[0] == 'Phases': |
550 | | |
551 | | for iphase in range(len(data)): |
552 | | |
553 | | if data[iphase][0] in Phases: |
554 | | |
555 | | phaseName = data[iphase][0] |
556 | | |
557 | | data[iphase][1].update(Phases[phaseName]) |
558 | | |
559 | | elif datum[0] == 'Covariance': |
560 | | |
561 | | data[0][1] = CovData |
562 | | |
563 | | try: |
564 | | |
565 | | histogram = Histograms[datum[0]] |
566 | | |
567 | | # print 'found ',datum[0] |
568 | | |
569 | | data[0][1][1] = histogram['Data'] |
570 | | |
571 | | for datus in data[1:]: |
572 | | |
573 | | # print ' read: ',datus[0] |
574 | | |
575 | | if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: |
576 | | |
577 | | datus[1] = histogram[datus[0]] |
578 | | |
579 | | except KeyError: |
580 | | |
581 | | pass |
582 | | |
583 | | |
584 | | |
585 | | cPickle.dump(data,outfile,1) |
586 | | |
587 | | infile.close() |
588 | | |
589 | | outfile.close() |
590 | | |
591 | | print 'GPX file save successful' |
592 | | |
593 | | |
594 | | |
595 | | def SetSeqResult(GPXfile,Histograms,SeqResult): |
596 | | |
597 | | GPXback = GPXBackup(GPXfile) |
598 | | |
599 | | print '\n',135*'-' |
600 | | |
601 | | print 'Read from file:',GPXback |
602 | | |
603 | | print 'Save to file :',GPXfile |
604 | | |
605 | | infile = open(GPXback,'rb') |
606 | | |
607 | | outfile = open(GPXfile,'wb') |
608 | | |
609 | | while True: |
610 | | |
611 | | try: |
612 | | |
613 | | data = cPickle.load(infile) |
614 | | |
615 | | except EOFError: |
616 | | |
617 | | break |
618 | | |
619 | | datum = data[0] |
620 | | |
621 | | if datum[0] == 'Sequental results': |
622 | | |
623 | | data[0][1] = SeqResult |
624 | | |
625 | | try: |
626 | | |
627 | | histogram = Histograms[datum[0]] |
628 | | |
629 | | data[0][1][1] = histogram['Data'] |
630 | | |
631 | | for datus in data[1:]: |
632 | | |
633 | | if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: |
634 | | |
635 | | datus[1] = histogram[datus[0]] |
636 | | |
637 | | except KeyError: |
638 | | |
639 | | pass |
640 | | |
641 | | |
642 | | |
643 | | cPickle.dump(data,outfile,1) |
644 | | |
645 | | infile.close() |
646 | | |
647 | | outfile.close() |
648 | | |
649 | | print 'GPX file save successful' |
650 | | |
651 | | |
652 | | |
653 | | def GetPWDRdata(GPXfile,PWDRname): |
654 | | |
655 | | ''' Returns powder data from GSASII gpx file |
656 | | |
657 | | input: |
658 | | |
659 | | GPXfile = .gpx full file name |
660 | | |
661 | | PWDRname = powder histogram name as obtained from GetHistogramNames |
662 | | |
663 | | return: |
664 | | |
665 | | PWDRdata = powder data dictionary with: |
666 | | |
667 | | Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters |
668 | | |
669 | | |
670 | | |
671 | | ''' |
672 | | |
673 | | file = open(GPXfile,'rb') |
674 | | |
675 | | PWDRdata = {} |
676 | | |
677 | | while True: |
678 | | |
679 | | try: |
680 | | |
681 | | data = cPickle.load(file) |
682 | | |
683 | | except EOFError: |
684 | | |
685 | | break |
686 | | |
687 | | datum = data[0] |
688 | | |
689 | | if datum[0] == PWDRname: |
690 | | |
691 | | PWDRdata['Data'] = datum[1][1] #powder data arrays |
692 | | |
693 | | PWDRdata[data[2][0]] = data[2][1] #Limits |
694 | | |
695 | | PWDRdata[data[3][0]] = data[3][1] #Background |
696 | | |
697 | | PWDRdata[data[4][0]] = data[4][1] #Instrument parameters |
698 | | |
699 | | PWDRdata[data[5][0]] = data[5][1] #Sample parameters |
700 | | |
| 2909 | covMatrix = result[1]*GOF |
| 2910 | sig = np.sqrt(np.diag(covMatrix)) |
| 2911 | if np.any(np.isnan(sig)): |
| 2912 | print '*** Least squares aborted - some invalid esds possible ***' |
| 2913 | # table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig))) |
| 2914 | # for item in table: print item,table[item] #useful debug - are things shifting? |
| 2915 | break #refinement succeeded - finish up! |
| 2916 | except TypeError: #result[1] is None on singular matrix |
| 2917 | print '**** Refinement failed - singular matrix ****' |
| 2918 | if 'Hessian' in Controls['deriv type']: |
| 2919 | for i in result[2]['psing'].reverse(): |
| 2920 | print 'Removing parameter: ',varyList[i] |
| 2921 | del(varyList[i]) |
| 2922 | else: |
| 2923 | Ipvt = result[2]['ipvt'] |
| 2924 | for i,ipvt in enumerate(Ipvt): |
| 2925 | if not np.sum(result[2]['fjac'],axis=1)[i]: |
| 2926 | print 'Removing parameter: ',varyList[ipvt-1] |
| 2927 | del(varyList[ipvt-1]) |
| 2928 | break |
| 2929 | |
| 2930 | # print 'dependentParmList: ',G2mv.dependentParmList |
| 2931 | # print 'arrayList: ',G2mv.arrayList |
| 2932 | # print 'invarrayList: ',G2mv.invarrayList |
| 2933 | # print 'indParmList: ',G2mv.indParmList |
| 2934 | # print 'fixedDict: ',G2mv.fixedDict |
| 2935 | # print 'test1' |
| 2936 | GetFobsSq(Histograms,Phases,parmDict,calcControls) |
| 2937 | # print 'test2' |
| 2938 | sigDict = dict(zip(varyList,sig)) |
| 2939 | newCellDict = GetNewCellParms(parmDict,varyList) |
| 2940 | newAtomDict = ApplyXYZshifts(parmDict,varyList) |
| 2941 | covData = {'variables':result[0],'varyList':varyList,'sig':sig, |
| 2942 | 'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict} |
| 2943 | # add the uncertainties into the esd dictionary (sigDict) |
| 2944 | sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict)) |
| 2945 | SetPhaseData(parmDict,sigDict,Phases,covData) |
| 2946 | SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms) |
| 2947 | SetHistogramData(parmDict,sigDict,Histograms) |
| 2948 | G2mv.PrintIndependentVars(parmDict,varyList,sigDict) |
| 2949 | SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData) |
| 2950 | |
| 2951 | #for testing purposes!!! |
| 2952 | # import cPickle |
| 2953 | # file = open('structTestdata.dat','wb') |
| 2954 | # cPickle.dump(parmDict,file,1) |
| 2955 | # cPickle.dump(varyList,file,1) |
| 2956 | # for histogram in Histograms: |
| 2957 | # if 'PWDR' in histogram[:4]: |
| 2958 | # Histogram = Histograms[histogram] |
| 2959 | # cPickle.dump(Histogram,file,1) |
| 2960 | # cPickle.dump(Phases,file,1) |
| 2961 | # cPickle.dump(calcControls,file,1) |
| 2962 | # cPickle.dump(pawleyLookup,file,1) |
| 2963 | # file.close() |
| 2964 | |
| 2965 | if dlg: |
| 2966 | return Rwp |
| 2967 | |
| 2968 | def SeqRefine(GPXfile,dlg): |
| 2969 | import pytexture as ptx |
| 2970 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
| 2971 | |
| 2972 | ShowBanner() |
| 2973 | print ' Sequential Refinement' |
| 2974 | G2mv.InitVars() |
| 2975 | Controls = GetControls(GPXfile) |
| 2976 | ShowControls(Controls) |
| 2977 | constrDict,constrFlag,fixedList = GetConstraints(GPXfile) |
| 2978 | Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) |
| 2979 | if not Phases: |
| 2980 | print ' *** ERROR - you have no histograms to refine! ***' |
| 2981 | print ' *** Refine aborted ***' |
| 2982 | raise Exception |
| 2983 | if not Histograms: |
| 2984 | print ' *** ERROR - you have no data to refine with! ***' |
| 2985 | print ' *** Refine aborted ***' |
| 2986 | raise Exception |
| 2987 | Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False) |
| 2988 | if 'Seq Data' in Controls: |
| 2989 | histNames = Controls['Seq Data'] |
| 2990 | else: |
| 2991 | histNames = GetHistogramNames(GPXfile,['PWDR',]) |
| 2992 | if 'Reverse Seq' in Controls: |
| 2993 | if Controls['Reverse Seq']: |
| 2994 | histNames.reverse() |
| 2995 | SeqResult = {'histNames':histNames} |
| 2996 | makeBack = True |
| 2997 | for ihst,histogram in enumerate(histNames): |
| 2998 | ifPrint = False |
| 2999 | if dlg: |
| 3000 | dlg.SetTitle('Residual for histogram '+str(ihst)) |
| 3001 | calcControls = {} |
| 3002 | calcControls['Natoms'] = Natoms |
| 3003 | calcControls['FFtables'] = FFtables |
| 3004 | calcControls['BLtables'] = BLtables |
| 3005 | varyList = [] |
| 3006 | parmDict = {} |
| 3007 | Histo = {histogram:Histograms[histogram],} |
| 3008 | hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False) |
| 3009 | calcControls.update(controlDict) |
| 3010 | histVary,histDict,controlDict = GetHistogramData(Histo,False) |
| 3011 | calcControls.update(controlDict) |
| 3012 | varyList = phaseVary+hapVary+histVary |
| 3013 | if not ihst: |
| 3014 | saveVaryList = varyList[:] |
| 3015 | for i,item in enumerate(saveVaryList): |
| 3016 | items = item.split(':') |
| 3017 | if items[1]: |
| 3018 | items[1] = '' |
| 3019 | item = ':'.join(items) |
| 3020 | saveVaryList[i] = item |
| 3021 | SeqResult['varyList'] = saveVaryList |
| 3022 | else: |
| 3023 | newVaryList = varyList[:] |
| 3024 | for i,item in enumerate(newVaryList): |
| 3025 | items = item.split(':') |
| 3026 | if items[1]: |
| 3027 | items[1] = '' |
| 3028 | item = ':'.join(items) |
| 3029 | newVaryList[i] = item |
| 3030 | if newVaryList != SeqResult['varyList']: |
| 3031 | print newVaryList |
| 3032 | print SeqResult['varyList'] |
| 3033 | print '**** ERROR - variable list for this histogram does not match previous' |
| 3034 | raise Exception |
| 3035 | parmDict.update(phaseDict) |
| 3036 | parmDict.update(hapDict) |
| 3037 | parmDict.update(histDict) |
| 3038 | GetFprime(calcControls,Histo) |
| 3039 | constrDict,constrFlag,fixedList = G2mv.InputParse([]) #constraints go here? |
| 3040 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
| 3041 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList) |
| 3042 | G2mv.Map2Dict(parmDict,varyList) |
| 3043 | |
| 3044 | while True: |
| 3045 | begin = time.time() |
| 3046 | values = np.array(Dict2Values(parmDict, varyList)) |
| 3047 | Ftol = Controls['min dM/M'] |
| 3048 | Factor = Controls['shift factor'] |
| 3049 | maxCyc = Controls['max cyc'] |
| 3050 | |
| 3051 | if 'Jacobian' in Controls['deriv type']: |
| 3052 | result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True, |
| 3053 | ftol=Ftol,col_deriv=True,factor=Factor, |
| 3054 | args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
| 3055 | ncyc = int(result[2]['nfev']/2) |
| 3056 | elif 'Hessian' in Controls['deriv type']: |
| 3057 | result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc, |
| 3058 | args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
| 3059 | ncyc = result[2]['num cyc']+1 |
| 3060 | else: #'numeric' |
| 3061 | result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor, |
| 3062 | args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
| 3063 | ncyc = int(result[2]['nfev']/len(varyList)) |
| 3064 | |
| 3065 | |
| 3066 | |
| 3067 | runtime = time.time()-begin |
| 3068 | chisq = np.sum(result[2]['fvec']**2) |
| 3069 | Values2Dict(parmDict, varyList, result[0]) |
| 3070 | G2mv.Dict2Map(parmDict,varyList) |
| 3071 | |
| 3072 | Rwp = np.sqrt(chisq/Histo['sumwYo'])*100. #to % |
| 3073 | GOF = chisq/(Histo['Nobs']-len(varyList)) |
| 3074 | print '\n Refinement results for histogram: v'+histogram |
| 3075 | print 135*'-' |
| 3076 | print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList) |
| 3077 | print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc) |
| 3078 | print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF) |
702 | | |
703 | | PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing |
704 | | |
705 | | except IndexError: |
706 | | |
707 | | PWDRdata['Reflection lists'] = {} |
708 | | |
709 | | file.close() |
710 | | |
711 | | return PWDRdata |
712 | | |
713 | | |
714 | | |
715 | | def GetHKLFdata(GPXfile,HKLFname): |
716 | | |
717 | | ''' Returns single crystal data from GSASII gpx file |
718 | | |
719 | | input: |
720 | | |
721 | | GPXfile = .gpx full file name |
722 | | |
723 | | HKLFname = single crystal histogram name as obtained from GetHistogramNames |
724 | | |
725 | | return: |
726 | | |
727 | | HKLFdata = single crystal data list of reflections: for each reflection: |
728 | | |
729 | | HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase] |
730 | | |
731 | | ''' |
732 | | |
733 | | file = open(GPXfile,'rb') |
734 | | |
735 | | HKLFdata = [] |
736 | | |
737 | | while True: |
738 | | |
739 | | try: |
740 | | |
741 | | data = cPickle.load(file) |
742 | | |
743 | | except EOFError: |
744 | | |
745 | | break |
746 | | |
747 | | datum = data[0] |
748 | | |
749 | | if datum[0] == HKLFname: |
750 | | |
751 | | HKLFdata = datum[1:][0] |
752 | | |
753 | | file.close() |
754 | | |
755 | | return HKLFdata |
756 | | |
757 | | |
758 | | |
759 | | def ShowBanner(): |
760 | | |
761 | | print 80*'*' |
762 | | |
763 | | print ' General Structure Analysis System-II Crystal Structure Refinement' |
764 | | |
765 | | print ' by Robert B. Von Dreele, Argonne National Laboratory(C), 2010' |
766 | | |
767 | | print ' This product includes software developed by the UChicago Argonne, LLC,' |
768 | | |
769 | | print ' as Operator of Argonne National Laboratory.' |
770 | | |
771 | | print 80*'*','\n' |
772 | | |
773 | | |
774 | | |
775 | | def ShowControls(Controls): |
776 | | |
777 | | print ' Least squares controls:' |
778 | | |
779 | | print ' Derivative type: ',Controls['deriv type'] |
780 | | |
781 | | print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M']) |
782 | | |
783 | | print ' Initial shift factor: ','%.3f'%(Controls['shift factor']) |
784 | | |
785 | | |
786 | | |
787 | | def GetFFtable(General): |
788 | | |
789 | | ''' returns a dictionary of form factor data for atom types found in General |
790 | | |
791 | | input: |
792 | | |
793 | | General = dictionary of phase info.; includes AtomTypes |
794 | | |
795 | | return: |
796 | | |
797 | | FFtable = dictionary of form factor data; key is atom type |
798 | | |
799 | | ''' |
800 | | |
801 | | atomTypes = General['AtomTypes'] |
802 | | |
803 | | FFtable = {} |
804 | | |
805 | | for El in atomTypes: |
806 | | |
807 | | FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0]) |
808 | | |
809 | | for item in FFs: |
810 | | |
811 | | if item['Symbol'] == El.upper(): |
812 | | |
813 | | FFtable[El] = item |
814 | | |
815 | | return FFtable |
816 | | |
817 | | |
818 | | |
819 | | def GetBLtable(General): |
820 | | |
821 | | ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General |
822 | | |
823 | | input: |
824 | | |
825 | | General = dictionary of phase info.; includes AtomTypes & Isotopes |
826 | | |
827 | | return: |
828 | | |
829 | | BLtable = dictionary of scattering length data; key is atom type |
830 | | |
831 | | ''' |
832 | | |
833 | | atomTypes = General['AtomTypes'] |
834 | | |
835 | | BLtable = {} |
836 | | |
837 | | isotopes = General['Isotopes'] |
838 | | |
839 | | isotope = General['Isotope'] |
840 | | |
841 | | for El in atomTypes: |
842 | | |
843 | | BLtable[El] = [isotope[El],isotopes[El][isotope[El]]] |
844 | | |
845 | | return BLtable |
846 | | |
847 | | |
848 | | |
849 | | def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary): |
850 | | |
851 | | if SGLaue in ['-1','2/m','mmm']: |
852 | | |
853 | | return #no Pawley symmetry required constraints |
854 | | |
855 | | for i,varyI in enumerate(pawleyVary): |
856 | | |
857 | | refI = int(varyI.split(':')[-1]) |
858 | | |
859 | | ih,ik,il = PawleyRef[refI][:3] |
860 | | |
861 | | for varyJ in pawleyVary[0:i]: |
862 | | |
863 | | refJ = int(varyJ.split(':')[-1]) |
864 | | |
865 | | jh,jk,jl = PawleyRef[refJ][:3] |
866 | | |
867 | | if SGLaue in ['4/m','4/mmm']: |
868 | | |
869 | | isum = ih**2+ik**2 |
870 | | |
871 | | jsum = jh**2+jk**2 |
872 | | |
873 | | if abs(il) == abs(jl) and isum == jsum: |
874 | | |
875 | | G2mv.StoreEquivalence(varyJ,(varyI,)) |
876 | | |
877 | | elif SGLaue in ['3R','3mR']: |
878 | | |
879 | | isum = ih**2+ik**2+il**2 |
880 | | |
881 | | jsum = jh**2+jk**2*jl**2 |
882 | | |
883 | | isum2 = ih*ik+ih*il+ik*il |
884 | | |
885 | | jsum2 = jh*jk+jh*jl+jk*jl |
886 | | |
887 | | if isum == jsum and isum2 == jsum2: |
888 | | |
889 | | G2mv.StoreEquivalence(varyJ,(varyI,)) |
890 | | |
891 | | elif SGLaue in ['3','3m1','31m','6/m','6/mmm']: |
892 | | |
893 | | isum = ih**2+ik**2+ih*ik |
894 | | |
895 | | jsum = jh**2+jk**2+jh*jk |
896 | | |
897 | | if abs(il) == abs(jl) and isum == jsum: |
898 | | |
899 | | G2mv.StoreEquivalence(varyJ,(varyI,)) |
900 | | |
901 | | elif SGLaue in ['m3','m3m']: |
902 | | |
903 | | isum = ih**2+ik**2+il**2 |
904 | | |
905 | | jsum = jh**2+jk**2+jl**2 |
906 | | |
907 | | if isum == jsum: |
908 | | |
909 | | G2mv.StoreEquivalence(varyJ,(varyI,)) |
910 | | |
911 | | |
912 | | |
913 | | def cellVary(pfx,SGData): |
914 | | |
915 | | if SGData['SGLaue'] in ['-1',]: |
916 | | |
917 | | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
918 | | |
919 | | elif SGData['SGLaue'] in ['2/m',]: |
920 | | |
921 | | if SGData['SGUniq'] == 'a': |
922 | | |
923 | | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3'] |
924 | | |
925 | | elif SGData['SGUniq'] == 'b': |
926 | | |
927 | | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4'] |
928 | | |
929 | | else: |
930 | | |
931 | | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5'] |
932 | | |
933 | | elif SGData['SGLaue'] in ['mmm',]: |
934 | | |
935 | | return [pfx+'A0',pfx+'A1',pfx+'A2'] |
936 | | |
937 | | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
938 | | |
939 | | return [pfx+'A0',pfx+'A2'] |
940 | | |
941 | | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
942 | | |
943 | | return [pfx+'A0',pfx+'A2'] |
944 | | |
945 | | elif SGData['SGLaue'] in ['3R', '3mR']: |
946 | | |
947 | | return [pfx+'A0',pfx+'A3'] |
948 | | |
949 | | elif SGData['SGLaue'] in ['m3m','m3']: |
950 | | |
951 | | return [pfx+'A0',] |
952 | | |
953 | | |
954 | | |
955 | | def GetPhaseData(PhaseData,Print=True): |
| 3080 | covMatrix = result[1]*GOF |
| 3081 | sig = np.sqrt(np.diag(covMatrix)) |
| 3082 | if np.any(np.isnan(sig)): |
| 3083 | print '*** Least squares aborted - some invalid esds possible ***' |
| 3084 | ifPrint = True |
| 3085 | break #refinement succeeded - finish up! |
| 3086 | except TypeError: #result[1] is None on singular matrix |
| 3087 | print '**** Refinement failed - singular matrix ****' |
| 3088 | if 'Hessian' in Controls['deriv type']: |
| 3089 | for i in result[2]['psing'].reverse(): |
| 3090 | print 'Removing parameter: ',varyList[i] |
| 3091 | del(varyList[i]) |
| 3092 | else: |
| 3093 | Ipvt = result[2]['ipvt'] |
| 3094 | for i,ipvt in enumerate(Ipvt): |
| 3095 | if not np.sum(result[2]['fjac'],axis=1)[i]: |
| 3096 | print 'Removing parameter: ',varyList[ipvt-1] |
| 3097 | del(varyList[ipvt-1]) |
| 3098 | break |
| 3099 | |
| 3100 | GetFobsSq(Histo,Phases,parmDict,calcControls) |
| 3101 | sigDict = dict(zip(varyList,sig)) |
| 3102 | newCellDict = GetNewCellParms(parmDict,varyList) |
| 3103 | newAtomDict = ApplyXYZshifts(parmDict,varyList) |
| 3104 | covData = {'variables':result[0],'varyList':varyList,'sig':sig, |
| 3105 | 'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict} |
| 3106 | SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint) |
| 3107 | SetHistogramData(parmDict,sigDict,Histo,ifPrint) |
| 3108 | SeqResult[histogram] = covData |
| 3109 | SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack) |
| 3110 | makeBack = False |
| 3111 | SetSeqResult(GPXfile,Histograms,SeqResult) |
| 3112 | |
| 3113 | def DistAngle(DisAglCtls,DisAglData): |
| 3114 | import numpy.ma as ma |
| 3115 | |
| 3116 | def ShowBanner(name): |
| 3117 | print 80*'*' |
| 3118 | print ' Interatomic Distances and Angles for phase '+name |
| 3119 | print 80*'*','\n' |
| 3120 | |
| 3121 | ShowBanner(DisAglCtls['Name']) |
| 3122 | SGData = DisAglData['SGData'] |
| 3123 | SGtext = G2spc.SGPrint(SGData) |
| 3124 | for line in SGtext: print line |
| 3125 | Cell = DisAglData['Cell'] |
| 3126 | |
| 3127 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
| 3128 | covData = {} |
| 3129 | if 'covData' in DisAglData: |
| 3130 | covData = DisAglData['covData'] |
| 3131 | covMatrix = covData['covMatrix'] |
| 3132 | varyList = covData['varyList'] |
| 3133 | pfx = str(DisAglData['pId'])+'::' |
| 3134 | A = G2lat.cell2A(Cell[:6]) |
| 3135 | cellSig = getCellEsd(pfx,SGData,A,covData) |
| 3136 | names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = '] |
| 3137 | valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)] |
| 3138 | line = '\n Unit cell:' |
| 3139 | for name,vals in zip(names,valEsd): |
| 3140 | line += name+vals |
| 3141 | print line |
| 3142 | else: |
| 3143 | print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \ |
| 3144 | ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \ |
| 3145 | '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6]) |
| 3146 | Factor = DisAglCtls['Factors'] |
| 3147 | Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii']))) |
| 3148 | Units = np.array([ #is there a nicer way to make this? |
| 3149 | [-1,-1,-1],[-1,-1,0],[-1,-1,1],[-1,0,-1],[-1,0,0],[-1,0,1],[-1,1,-1],[-1,1,0],[-1,1,1], |
| 3150 | [0,-1,-1],[0,-1,0],[0,-1,1],[0,0,-1],[0,0,0],[0,0,1],[0,1,-1],[0,1,0],[0,1,1], |
| 3151 | [1,-1,-1],[1,-1,0],[1,-1,1],[1,0,-1],[1,0,0],[1,0,1],[1,1,-1],[1,1,0],[1,1,1]]) |
| 3152 | origAtoms = DisAglData['OrigAtoms'] |
| 3153 | targAtoms = DisAglData['TargAtoms'] |
| 3154 | for Oatom in origAtoms: |
| 3155 | OxyzNames = '' |
| 3156 | IndBlist = [] |
| 3157 | Dist = [] |
| 3158 | Vect = [] |
| 3159 | VectA = [] |
| 3160 | angles = [] |
| 3161 | for Tatom in targAtoms: |
| 3162 | Xvcov = [] |
| 3163 | TxyzNames = '' |
| 3164 | if 'covData' in DisAglData: |
| 3165 | OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])] |
| 3166 | TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])] |
| 3167 | Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix) |
| 3168 | result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False) |
| 3169 | BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0] |
| 3170 | AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1] |
| 3171 | for Txyz,Top,Tunit in result: |
| 3172 | Dx = (Txyz-np.array(Oatom[3:6]))+Units |
| 3173 | dx = np.inner(Amat,Dx) |
| 3174 | dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5) |
| 3175 | IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.)) |
| 3176 | if np.any(IndB): |
| 3177 | for indb in IndB: |
| 3178 | for i in range(len(indb)): |
| 3179 | if str(dx.T[indb][i]) not in IndBlist: |
| 3180 | IndBlist.append(str(dx.T[indb][i])) |
| 3181 | unit = Units[indb][i] |
| 3182 | tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2]) |
| 3183 | pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData) |
| 3184 | sig = 0.0 |
| 3185 | if len(Xvcov): |
| 3186 | sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx))) |
| 3187 | Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig]) |
| 3188 | if (Dist[-1][-1]-AsumR) <= 0.: |
| 3189 | Vect.append(dx.T[indb][i]/Dist[-1][-2]) |
| 3190 | VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top]) |
| 3191 | else: |
| 3192 | Vect.append([0.,0.,0.]) |
| 3193 | VectA.append([]) |
| 3194 | Vect = np.array(Vect) |
| 3195 | angles = np.zeros((len(Vect),len(Vect))) |
| 3196 | angsig = np.zeros((len(Vect),len(Vect))) |
| 3197 | for i,veca in enumerate(Vect): |
| 3198 | if np.any(veca): |
| 3199 | for j,vecb in enumerate(Vect): |
| 3200 | if np.any(vecb): |
| 3201 | angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData) |
| 3202 | line = '' |
| 3203 | for i,x in enumerate(Oatom[3:6]): |
| 3204 | if len(Xvcov): |
| 3205 | line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True)) |
| 3206 | else: |
| 3207 | line += '%12.5f'%(x) |
| 3208 | print '\n Distances & angles for ',Oatom[1],' at ',line |
| 3209 | print 80*'*' |
| 3210 | line = '' |
| 3211 | for dist in Dist[:-1]: |
| 3212 | line += '%12s'%(dist[1].center(12)) |
| 3213 | print ' To cell +(sym. op.) dist. ',line |
| 3214 | for i,dist in enumerate(Dist): |
| 3215 | line = '' |
| 3216 | for j,angle in enumerate(angles[i][0:i]): |
| 3217 | sig = angsig[i][j] |
| 3218 | if angle: |
| 3219 | if sig: |
| 3220 | line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12)) |
| 3221 | else: |
| 3222 | val = '%.3f'%(angle) |
| 3223 | line += '%12s'%(val.center(12)) |
| 3224 | else: |
| 3225 | line += 12*' ' |
| 3226 | if dist[5]: #sig exists! |
| 3227 | val = G2mth.ValEsd(dist[4],dist[5]) |
| 3228 | else: |
| 3229 | val = '%8.4f'%(dist[4]) |
| 3230 | print ' %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line |
| 3231 | |
| 3232 | def DisAglTor(DATData): |
| 3233 | SGData = DATData['SGData'] |
| 3234 | Cell = DATData['Cell'] |
| 3235 | |
| 3236 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
| 3237 | covData = {} |
| 3238 | pfx = '' |
| 3239 | if 'covData' in DATData: |
| 3240 | covData = DATData['covData'] |
| 3241 | covMatrix = covData['covMatrix'] |
| 3242 | varyList = covData['varyList'] |
| 3243 | pfx = str(DATData['pId'])+'::' |
| 3244 | Datoms = [] |
| 3245 | Oatoms = [] |
| 3246 | for i,atom in enumerate(DATData['Datoms']): |
| 3247 | symop = atom[-1].split('+') |
| 3248 | if len(symop) == 1: |
| 3249 | symop.append('0,0,0') |
| 3250 | symop[0] = int(symop[0]) |
| 3251 | symop[1] = eval(symop[1]) |
| 3252 | atom.append(symop) |
| 3253 | Datoms.append(atom) |
| 3254 | oatom = DATData['Oatoms'][i] |
| 3255 | names = ['','',''] |
| 3256 | if pfx: |
| 3257 | names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])] |
| 3258 | oatom += [names,] |
| 3259 | Oatoms.append(oatom) |
| 3260 | atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms] |
| 3261 | if DATData['Natoms'] == 4: #torsion |
| 3262 | Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData) |
| 3263 | print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig) |
| 3264 | print ' NB: Atom sequence determined by selection order' |
| 3265 | return # done with torsion |
| 3266 | elif DATData['Natoms'] == 3: #angle |
| 3267 | Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData) |
| 3268 | print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig) |
| 3269 | print ' NB: Atom sequence determined by selection order' |
| 3270 | else: #2 atoms - distance |
| 3271 | Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData) |
| 3272 | print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig) |
| 3273 | |
| 3274 | def BestPlane(PlaneData): |
| 3275 | |
| 3276 | def ShowBanner(name): |
| 3277 | print 80*'*' |
| 3278 | print ' Best plane result for phase '+name |
| 3279 | print 80*'*','\n' |
| 3280 | |
| 3281 | ShowBanner(PlaneData['Name']) |
| 3282 | |
| 3283 | Cell = PlaneData['Cell'] |
| 3284 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
| 3285 | Atoms = PlaneData['Atoms'] |
| 3286 | sumXYZ = np.zeros(3) |
| 3287 | XYZ = [] |
| 3288 | Natoms = len(Atoms) |
| 3289 | for atom in Atoms: |
| 3290 | xyz = np.array(atom[3:6]) |
| 3291 | XYZ.append(xyz) |
| 3292 | sumXYZ += xyz |
| 3293 | sumXYZ /= Natoms |
| 3294 | XYZ = np.array(XYZ)-sumXYZ |
| 3295 | XYZ = np.inner(Amat,XYZ).T |
| 3296 | Zmat = np.zeros((3,3)) |
| 3297 | for i,xyz in enumerate(XYZ): |
| 3298 | Zmat += np.outer(xyz.T,xyz) |
| 3299 | print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]) |
| 3300 | Evec,Emat = nl.eig(Zmat) |
| 3301 | Evec = np.sqrt(Evec)/(Natoms-3) |
| 3302 | Order = np.argsort(Evec) |
| 3303 | XYZ = np.inner(XYZ,Emat.T).T |
| 3304 | XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T |
| 3305 | print ' Atoms in Cartesian best plane coordinates:' |
| 3306 | print ' Name X Y Z' |
| 3307 | for i,xyz in enumerate(XYZ): |
| 3308 | print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]) |
| 3309 | print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]) |
958 | | |
959 | | def PrintFFtable(FFtable): |
960 | | |
961 | | print '\n X-ray scattering factors:' |
962 | | |
963 | | print ' Symbol fa fb fc' |
964 | | |
965 | | print 99*'-' |
966 | | |
967 | | for Ename in FFtable: |
968 | | |
969 | | ffdata = FFtable[Ename] |
970 | | |
971 | | fa = ffdata['fa'] |
972 | | |
973 | | fb = ffdata['fb'] |
974 | | |
975 | | print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \ |
976 | | |
977 | | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc']) |
978 | | |
979 | | |
980 | | |
981 | | def PrintBLtable(BLtable): |
982 | | |
983 | | print '\n Neutron scattering factors:' |
984 | | |
985 | | print ' Symbol isotope mass b resonant terms' |
986 | | |
987 | | print 99*'-' |
988 | | |
989 | | for Ename in BLtable: |
990 | | |
991 | | bldata = BLtable[Ename] |
992 | | |
993 | | isotope = bldata[0] |
994 | | |
995 | | mass = bldata[1][0] |
996 | | |
997 | | blen = bldata[1][1] |
998 | | |
999 | | bres = [] |
1000 | | |
1001 | | if len(bldata[1]) > 2: |
1002 | | |
1003 | | bres = bldata[1][2:] |
1004 | | |
1005 | | line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen) |
1006 | | |
1007 | | for item in bres: |
1008 | | |
1009 | | line += '%10.5g'%(item) |
1010 | | |
1011 | | print line |
1012 | | |
1013 | | |
1014 | | |
1015 | | def PrintAtoms(General,Atoms): |
1016 | | |
1017 | | print '\n Atoms:' |
1018 | | |
1019 | | line = ' name type refine? x y z '+ \ |
1020 | | |
1021 | | ' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23' |
1022 | | |
1023 | | if General['Type'] == 'magnetic': |
1024 | | |
1025 | | line += ' Mx My Mz' |
1026 | | |
1027 | | elif General['Type'] == 'macromolecular': |
1028 | | |
1029 | | line = ' res no residue chain '+line |
1030 | | |
1031 | | print line |
1032 | | |
1033 | | if General['Type'] == 'nuclear': |
1034 | | |
1035 | | print 135*'-' |
1036 | | |
1037 | | for i,at in enumerate(Atoms): |
1038 | | |
1039 | | line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \ |
1040 | | |
1041 | | '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9]) |
1042 | | |
1043 | | if at[9] == 'I': |
1044 | | |
1045 | | line += '%8.4f'%(at[10])+48*' ' |
1046 | | |
1047 | | else: |
1048 | | |
1049 | | line += 8*' ' |
1050 | | |
1051 | | for j in range(6): |
1052 | | |
1053 | | line += '%8.4f'%(at[11+j]) |
1054 | | |
1055 | | print line |
1056 | | |
1057 | | |
1058 | | |
1059 | | def PrintTexture(textureData): |
1060 | | |
1061 | | topstr = '\n Spherical harmonics texture: Order:' + \ |
1062 | | |
1063 | | str(textureData['Order']) |
1064 | | |
1065 | | if textureData['Order']: |
1066 | | |
1067 | | print topstr+' Refine? '+str(textureData['SH Coeff'][0]) |
1068 | | |
1069 | | else: |
1070 | | |
1071 | | print topstr |
1072 | | |
1073 | | return |
1074 | | |
1075 | | names = ['omega','chi','phi'] |
1076 | | |
1077 | | line = '\n' |
1078 | | |
1079 | | for name in names: |
1080 | | |
1081 | | line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0]) |
1082 | | |
1083 | | print line |
1084 | | |
1085 | | print '\n Texture coefficients:' |
1086 | | |
1087 | | ptlbls = ' names :' |
1088 | | |
1089 | | ptstr = ' values:' |
1090 | | |
1091 | | SHcoeff = textureData['SH Coeff'][1] |
1092 | | |
1093 | | for item in SHcoeff: |
1094 | | |
1095 | | ptlbls += '%12s'%(item) |
1096 | | |
1097 | | ptstr += '%12.4f'%(SHcoeff[item]) |
1098 | | |
1099 | | print ptlbls |
1100 | | |
1101 | | print ptstr |
1102 | | |
1103 | | |
1104 | | |
1105 | | if Print: print ' Phases:' |
1106 | | |
1107 | | phaseVary = [] |
1108 | | |
1109 | | phaseDict = {} |
1110 | | |
1111 | | phaseConstr = {} |
1112 | | |
1113 | | pawleyLookup = {} |
1114 | | |
1115 | | FFtables = {} #scattering factors - xrays |
1116 | | |
1117 | | BLtables = {} # neutrons |
1118 | | |
1119 | | Natoms = {} |
1120 | | |
1121 | | AtMults = {} |
1122 | | |
1123 | | AtIA = {} |
1124 | | |
1125 | | shModels = ['cylindrical','none','shear - 2/m','rolling - mmm'] |
1126 | | |
1127 | | SamSym = dict(zip(shModels,['0','-1','2/m','mmm'])) |
1128 | | |
1129 | | for name in PhaseData: |
1130 | | |
1131 | | General = PhaseData[name]['General'] |
1132 | | |
1133 | | pId = PhaseData[name]['pId'] |
1134 | | |
1135 | | pfx = str(pId)+'::' |
1136 | | |
1137 | | FFtable = GetFFtable(General) |
1138 | | |
1139 | | BLtable = GetBLtable(General) |
1140 | | |
1141 | | FFtables.update(FFtable) |
1142 | | |
1143 | | BLtables.update(BLtable) |
1144 | | |
1145 | | Atoms = PhaseData[name]['Atoms'] |
1146 | | |
1147 | | try: |
1148 | | |
1149 | | PawleyRef = PhaseData[name]['Pawley ref'] |
1150 | | |
1151 | | except KeyError: |
1152 | | |
1153 | | PawleyRef = [] |
1154 | | |
1155 | | SGData = General['SGData'] |
1156 | | |
1157 | | SGtext = G2spc.SGPrint(SGData) |
1158 | | |
1159 | | cell = General['Cell'] |
1160 | | |
1161 | | A = G2lat.cell2A(cell[1:7]) |
1162 | | |
1163 | | phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]}) |
1164 | | |
1165 | | if cell[0]: |
1166 | | |
1167 | | phaseVary += cellVary(pfx,SGData) |
1168 | | |
1169 | | Natoms[pfx] = 0 |
1170 | | |
1171 | | if Atoms: |
1172 | | |
1173 | | if General['Type'] == 'nuclear': |
1174 | | |
1175 | | Natoms[pfx] = len(Atoms) |
1176 | | |
1177 | | for i,at in enumerate(Atoms): |
1178 | | |
1179 | | phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8], |
1180 | | |
1181 | | pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5], |
1182 | | |
1183 | | pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z |
1184 | | |
1185 | | pfx+'AI/A:'+str(i):at[9],}) |
1186 | | |
1187 | | if at[9] == 'I': |
1188 | | |
1189 | | phaseDict[pfx+'AUiso:'+str(i)] = at[10] |
1190 | | |
1191 | | else: |
1192 | | |
1193 | | phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13], |
1194 | | |
1195 | | pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]}) |
1196 | | |
1197 | | if 'F' in at[2]: |
1198 | | |
1199 | | phaseVary.append(pfx+'Afrac:'+str(i)) |
1200 | | |
1201 | | if 'X' in at[2]: |
1202 | | |
1203 | | xId,xCoef = G2spc.GetCSxinel(at[7]) |
1204 | | |
1205 | | delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)] |
1206 | | |
1207 | | for j in range(3): |
1208 | | |
1209 | | if xId[j] > 0: |
1210 | | |
1211 | | phaseVary.append(delnames[j]) |
1212 | | |
1213 | | for k in range(j): |
1214 | | |
1215 | | if xId[j] == xId[k]: |
1216 | | |
1217 | | G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) |
1218 | | |
1219 | | if 'U' in at[2]: |
1220 | | |
1221 | | if at[9] == 'I': |
1222 | | |
1223 | | phaseVary.append(pfx+'AUiso:'+str(i)) |
1224 | | |
1225 | | else: |
1226 | | |
1227 | | uId,uCoef = G2spc.GetCSuinel(at[7])[:2] |
1228 | | |
1229 | | names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i), |
1230 | | |
1231 | | pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)] |
1232 | | |
1233 | | for j in range(6): |
1234 | | |
1235 | | if uId[j] > 0: |
1236 | | |
1237 | | phaseVary.append(names[j]) |
1238 | | |
1239 | | for k in range(j): |
1240 | | |
1241 | | if uId[j] == uId[k]: |
1242 | | |
1243 | | G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),)) |
1244 | | |
1245 | | # elif General['Type'] == 'magnetic': |
1246 | | |
1247 | | # elif General['Type'] == 'macromolecular': |
1248 | | |
1249 | | |
1250 | | |
1251 | | |
1252 | | |
1253 | | if 'SH Texture' in General: |
1254 | | |
1255 | | textureData = General['SH Texture'] |
1256 | | |
1257 | | phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']] |
1258 | | |
1259 | | phaseDict[pfx+'SHorder'] = textureData['Order'] |
1260 | | |
1261 | | for name in ['omega','chi','phi']: |
1262 | | |
1263 | | phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1] |
1264 | | |
1265 | | if textureData['Sample '+name][0]: |
1266 | | |
1267 | | phaseVary.append(pfx+'SH '+name) |
1268 | | |
1269 | | for name in textureData['SH Coeff'][1]: |
1270 | | |
1271 | | phaseDict[pfx+name] = textureData['SH Coeff'][1][name] |
1272 | | |
1273 | | if textureData['SH Coeff'][0]: |
1274 | | |
1275 | | phaseVary.append(pfx+name) |
1276 | | |
1277 | | |
1278 | | |
1279 | | if Print: |
1280 | | |
1281 | | print '\n Phase name: ',General['Name'] |
1282 | | |
1283 | | print 135*'-' |
1284 | | |
1285 | | PrintFFtable(FFtable) |
1286 | | |
1287 | | PrintBLtable(BLtable) |
1288 | | |
1289 | | print '' |
1290 | | |
1291 | | for line in SGtext: print line |
1292 | | |
1293 | | PrintAtoms(General,Atoms) |
1294 | | |
1295 | | print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \ |
1296 | | |
1297 | | ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \ |
1298 | | |
1299 | | '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0] |
1300 | | |
1301 | | if 'SH Texture' in General: |
1302 | | |
1303 | | PrintTexture(textureData) |
1304 | | |
1305 | | |
1306 | | |
1307 | | elif PawleyRef: |
1308 | | |
1309 | | pawleyVary = [] |
1310 | | |
1311 | | for i,refl in enumerate(PawleyRef): |
1312 | | |
1313 | | phaseDict[pfx+'PWLref:'+str(i)] = refl[6] |
1314 | | |
1315 | | pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i |
1316 | | |
1317 | | if refl[5]: |
1318 | | |
1319 | | pawleyVary.append(pfx+'PWLref:'+str(i)) |
1320 | | |
1321 | | GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence |
1322 | | |
1323 | | phaseVary += pawleyVary |
1324 | | |
1325 | | |
1326 | | |
1327 | | return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables |
1328 | | |
1329 | | |
1330 | | |
1331 | | def cellFill(pfx,SGData,parmDict,sigDict): |
1332 | | |
1333 | | if SGData['SGLaue'] in ['-1',]: |
1334 | | |
1335 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
1336 | | |
1337 | | parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']] |
1338 | | |
1339 | | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
1340 | | |
1341 | | sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']] |
1342 | | |
1343 | | elif SGData['SGLaue'] in ['2/m',]: |
1344 | | |
1345 | | if SGData['SGUniq'] == 'a': |
1346 | | |
1347 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
1348 | | |
1349 | | parmDict[pfx+'A3'],0,0] |
1350 | | |
1351 | | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
1352 | | |
1353 | | sigDict[pfx+'A3'],0,0] |
1354 | | |
1355 | | elif SGData['SGUniq'] == 'b': |
1356 | | |
1357 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
1358 | | |
1359 | | 0,parmDict[pfx+'A4'],0] |
1360 | | |
1361 | | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
1362 | | |
1363 | | 0,sigDict[pfx+'A4'],0] |
1364 | | |
1365 | | else: |
1366 | | |
1367 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
1368 | | |
1369 | | 0,0,parmDict[pfx+'A5']] |
1370 | | |
1371 | | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
1372 | | |
1373 | | 0,0,sigDict[pfx+'A5']] |
1374 | | |
1375 | | elif SGData['SGLaue'] in ['mmm',]: |
1376 | | |
1377 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0] |
1378 | | |
1379 | | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0] |
1380 | | |
1381 | | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
1382 | | |
1383 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0] |
1384 | | |
1385 | | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
1386 | | |
1387 | | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
1388 | | |
1389 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'], |
1390 | | |
1391 | | parmDict[pfx+'A0'],0,0] |
1392 | | |
1393 | | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
1394 | | |
1395 | | elif SGData['SGLaue'] in ['3R', '3mR']: |
1396 | | |
1397 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'], |
1398 | | |
1399 | | parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']] |
1400 | | |
1401 | | sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0] |
1402 | | |
1403 | | elif SGData['SGLaue'] in ['m3m','m3']: |
1404 | | |
1405 | | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0] |
1406 | | |
1407 | | sigA = [sigDict[pfx+'A0'],0,0,0,0,0] |
1408 | | |
1409 | | return A,sigA |
1410 | | |
1411 | | |
1412 | | |
1413 | | def getCellEsd(pfx,SGData,A,covData): |
1414 | | |
1415 | | dpr = 180./np.pi |
1416 | | |
1417 | | rVsq = G2lat.calc_rVsq(A) |
1418 | | |
1419 | | G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors |
1420 | | |
1421 | | cell = np.array(G2lat.Gmat2cell(g)) #real cell |
1422 | | |
1423 | | cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell |
1424 | | |
1425 | | scos = cosd(cellst[3:6]) |
1426 | | |
1427 | | ssin = sind(cellst[3:6]) |
1428 | | |
1429 | | scot = scos/ssin |
1430 | | |
1431 | | rcos = cosd(cell[3:6]) |
1432 | | |
1433 | | rsin = sind(cell[3:6]) |
1434 | | |
1435 | | rcot = rcos/rsin |
1436 | | |
1437 | | RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
1438 | | |
1439 | | varyList = covData['varyList'] |
1440 | | |
1441 | | covMatrix = covData['covMatrix'] |
1442 | | |
1443 | | vcov = G2mth.getVCov(RMnames,varyList,covMatrix) |
1444 | | |
1445 | | Ax = np.array(A) |
1446 | | |
1447 | | Ax[3:] /= 2. |
1448 | | |
1449 | | drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2, |
1450 | | |
1451 | | Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]]) |
1452 | | |
1453 | | srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T)) |
1454 | | |
1455 | | Vol = 1/np.sqrt(rVsq) |
1456 | | |
1457 | | sigVol = Vol**3*np.sqrt(srcvlsq)/2. |
1458 | | |
1459 | | R123 = Ax[0]*Ax[1]*Ax[2] |
1460 | | |
1461 | | dsasdg = np.zeros((3,6)) |
1462 | | |
1463 | | dadg = np.zeros((6,6)) |
1464 | | |
1465 | | for i0 in range(3): #0 1 2 |
1466 | | |
1467 | | i1 = (i0+1)%3 #1 2 0 |
1468 | | |
1469 | | i2 = (i1+1)%3 #2 0 1 |
1470 | | |
1471 | | i3 = 5-i2 #3 5 4 |
1472 | | |
1473 | | i4 = 5-i1 #4 3 5 |
1474 | | |
1475 | | i5 = 5-i0 #5 4 3 |
1476 | | |
1477 | | dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1] |
1478 | | |
1479 | | dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2] |
1480 | | |
1481 | | dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2]) |
1482 | | |
1483 | | denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2) |
1484 | | |
1485 | | denom = np.sqrt(denmsq) |
1486 | | |
1487 | | dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2) |
1488 | | |
1489 | | dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2) |
1490 | | |
1491 | | dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2) |
1492 | | |
1493 | | dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2) |
1494 | | |
1495 | | dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4]) |
1496 | | |
1497 | | dadg[i5][i5] = -Ax[i0]/denom |
1498 | | |
1499 | | for i0 in range(3): |
1500 | | |
1501 | | i1 = (i0+1)%3 |
1502 | | |
1503 | | i2 = (i1+1)%3 |
1504 | | |
1505 | | i3 = 5-i2 |
1506 | | |
1507 | | for ij in range(6): |
1508 | | |
1509 | | dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1]) |
1510 | | |
1511 | | if ij == i0: |
1512 | | |
1513 | | dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0] |
1514 | | |
1515 | | dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr |
1516 | | |
1517 | | sigMat = np.inner(dadg,np.inner(vcov,dadg.T)) |
1518 | | |
1519 | | var = np.diag(sigMat) |
1520 | | |
1521 | | CS = np.where(var>0.,np.sqrt(var),0.) |
1522 | | |
1523 | | cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order |
1524 | | |
1525 | | return cellSig |
1526 | | |
1527 | | |
1528 | | |
1529 | | def SetPhaseData(parmDict,sigDict,Phases,covData): |
1530 | | |
1531 | | |
1532 | | |
1533 | | def PrintAtomsAndSig(General,Atoms,atomsSig): |
1534 | | |
1535 | | print '\n Atoms:' |
1536 | | |
1537 | | line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23' |
1538 | | |
1539 | | if General['Type'] == 'magnetic': |
1540 | | |
1541 | | line += ' Mx My Mz' |
1542 | | |
1543 | | elif General['Type'] == 'macromolecular': |
1544 | | |
1545 | | line = ' res no residue chain '+line |
1546 | | |
1547 | | print line |
1548 | | |
1549 | | if General['Type'] == 'nuclear': |
1550 | | |
1551 | | print 135*'-' |
1552 | | |
1553 | | fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f', |
1554 | | |
1555 | | 11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'} |
1556 | | |
1557 | | noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']} |
1558 | | |
1559 | | for i,at in enumerate(Atoms): |
1560 | | |
1561 | | name = fmt[0]%(at[0])+fmt[1]%(at[1])+':' |
1562 | | |
1563 | | valstr = ' values:' |
1564 | | |
1565 | | sigstr = ' sig :' |
1566 | | |
1567 | | for ind in [3,4,5,6]: |
1568 | | |
1569 | | sigind = str(i)+':'+str(ind) |
1570 | | |
1571 | | valstr += fmt[ind]%(at[ind]) |
1572 | | |
1573 | | if sigind in atomsSig: |
1574 | | |
1575 | | sigstr += fmt[ind]%(atomsSig[sigind]) |
1576 | | |
1577 | | else: |
1578 | | |
1579 | | sigstr += noFXsig[ind][1]%(noFXsig[ind][0]) |
1580 | | |
1581 | | if at[9] == 'I': |
1582 | | |
1583 | | valstr += fmt[10]%(at[10]) |
1584 | | |
1585 | | if str(i)+':10' in atomsSig: |
1586 | | |
1587 | | sigstr += fmt[10]%(atomsSig[str(i)+':10']) |
1588 | | |
1589 | | else: |
1590 | | |
1591 | | sigstr += 8*' ' |
1592 | | |
1593 | | else: |
1594 | | |
1595 | | valstr += 8*' ' |
1596 | | |
1597 | | sigstr += 8*' ' |
1598 | | |
1599 | | for ind in [11,12,13,14,15,16]: |
1600 | | |
1601 | | sigind = str(i)+':'+str(ind) |
1602 | | |
1603 | | valstr += fmt[ind]%(at[ind]) |
1604 | | |
1605 | | if sigind in atomsSig: |
1606 | | |
1607 | | sigstr += fmt[ind]%(atomsSig[sigind]) |
1608 | | |
1609 | | else: |
1610 | | |
1611 | | sigstr += 8*' ' |
1612 | | |
1613 | | print name |
1614 | | |
1615 | | print valstr |
1616 | | |
1617 | | print sigstr |
1618 | | |
1619 | | |
1620 | | |
1621 | | def PrintSHtextureAndSig(textureData,SHtextureSig): |
1622 | | |
1623 | | print '\n Spherical harmonics texture: Order:' + str(textureData['Order']) |
1624 | | |
1625 | | names = ['omega','chi','phi'] |
1626 | | |
1627 | | namstr = ' names :' |
1628 | | |
1629 | | ptstr = ' values:' |
1630 | | |
1631 | | sigstr = ' esds :' |
1632 | | |
1633 | | for name in names: |
1634 | | |
1635 | | namstr += '%12s'%(name) |
1636 | | |
1637 | | ptstr += '%12.3f'%(textureData['Sample '+name][1]) |
1638 | | |
1639 | | if 'Sample '+name in SHtextureSig: |
1640 | | |
1641 | | sigstr += '%12.3f'%(SHtextureSig['Sample '+name]) |
1642 | | |
1643 | | else: |
1644 | | |
1645 | | sigstr += 12*' ' |
1646 | | |
1647 | | print namstr |
1648 | | |
1649 | | print ptstr |
1650 | | |
1651 | | print sigstr |
1652 | | |
1653 | | print '\n Texture coefficients:' |
1654 | | |
1655 | | namstr = ' names :' |
1656 | | |
1657 | | ptstr = ' values:' |
1658 | | |
1659 | | sigstr = ' esds :' |
1660 | | |
1661 | | SHcoeff = textureData['SH Coeff'][1] |
1662 | | |
1663 | | for name in SHcoeff: |
1664 | | |
1665 | | namstr += '%12s'%(name) |
1666 | | |
1667 | | ptstr += '%12.3f'%(SHcoeff[name]) |
1668 | | |
1669 | | if name in SHtextureSig: |
1670 | | |
1671 | | sigstr += '%12.3f'%(SHtextureSig[name]) |
1672 | | |
1673 | | else: |
1674 | | |
1675 | | sigstr += 12*' ' |
1676 | | |
1677 | | print namstr |
1678 | | |
1679 | | print ptstr |
1680 | | |
1681 | | print sigstr |
1682 | | |
1683 | | |
1684 | | |
1685 | | |
1686 | | |
1687 | | print '\n Phases:' |
1688 | | |
1689 | | for phase in Phases: |
1690 | | |
1691 | | print ' Result for phase: ',phase |
1692 | | |
1693 | | Phase = Phases[phase] |
1694 | | |
1695 | | General = Phase['General'] |
1696 | | |
1697 | | SGData = General['SGData'] |
1698 | | |
1699 | | Atoms = Phase['Atoms'] |
1700 | | |
1701 | | cell = General['Cell'] |
1702 | | |
1703 | | pId = Phase['pId'] |
1704 | | |
1705 | | pfx = str(pId)+'::' |
1706 | | |
1707 | | if cell[0]: |
1708 | | |
1709 | | A,sigA = cellFill(pfx,SGData,parmDict,sigDict) |
1710 | | |
1711 | | cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol |
1712 | | |
1713 | | print ' Reciprocal metric tensor: ' |
1714 | | |
1715 | | ptfmt = "%15.9f" |
1716 | | |
1717 | | names = ['A11','A22','A33','A12','A13','A23'] |
1718 | | |
1719 | | namstr = ' names :' |
1720 | | |
1721 | | ptstr = ' values:' |
1722 | | |
1723 | | sigstr = ' esds :' |
1724 | | |
1725 | | for name,a,siga in zip(names,A,sigA): |
1726 | | |
1727 | | namstr += '%15s'%(name) |
1728 | | |
1729 | | ptstr += ptfmt%(a) |
1730 | | |
1731 | | if siga: |
1732 | | |
1733 | | sigstr += ptfmt%(siga) |
1734 | | |
1735 | | else: |
1736 | | |
1737 | | sigstr += 15*' ' |
1738 | | |
1739 | | print namstr |
1740 | | |
1741 | | print ptstr |
1742 | | |
1743 | | print sigstr |
1744 | | |
1745 | | cell[1:7] = G2lat.A2cell(A) |
1746 | | |
1747 | | cell[7] = G2lat.calc_V(A) |
1748 | | |
1749 | | print ' New unit cell:' |
1750 | | |
1751 | | ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"] |
1752 | | |
1753 | | names = ['a','b','c','alpha','beta','gamma','Volume'] |
1754 | | |
1755 | | namstr = ' names :' |
1756 | | |
1757 | | ptstr = ' values:' |
1758 | | |
1759 | | sigstr = ' esds :' |
1760 | | |
1761 | | for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig): |
1762 | | |
1763 | | namstr += '%12s'%(name) |
1764 | | |
1765 | | ptstr += fmt%(a) |
1766 | | |
1767 | | if siga: |
1768 | | |
1769 | | sigstr += fmt%(siga) |
1770 | | |
1771 | | else: |
1772 | | |
1773 | | sigstr += 12*' ' |
1774 | | |
1775 | | print namstr |
1776 | | |
1777 | | print ptstr |
1778 | | |
1779 | | print sigstr |
1780 | | |
1781 | | |
1782 | | |
1783 | | if 'Pawley' in Phase['General']['Type']: |
1784 | | |
1785 | | pawleyRef = Phase['Pawley ref'] |
1786 | | |
1787 | | for i,refl in enumerate(pawleyRef): |
1788 | | |
1789 | | key = pfx+'PWLref:'+str(i) |
1790 | | |
1791 | | refl[6] = abs(parmDict[key]) #suppress negative Fsq |
1792 | | |
1793 | | if key in sigDict: |
1794 | | |
1795 | | refl[7] = sigDict[key] |
1796 | | |
1797 | | else: |
1798 | | |
1799 | | refl[7] = 0 |
1800 | | |
1801 | | else: |
1802 | | |
1803 | | atomsSig = {} |
1804 | | |
1805 | | if General['Type'] == 'nuclear': |
1806 | | |
1807 | | for i,at in enumerate(Atoms): |
1808 | | |
1809 | | names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i), |
1810 | | |
1811 | | 10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i), |
1812 | | |
1813 | | 14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)} |
1814 | | |
1815 | | for ind in [3,4,5,6]: |
1816 | | |
1817 | | at[ind] = parmDict[names[ind]] |
1818 | | |
1819 | | if ind in [3,4,5]: |
1820 | | |
1821 | | name = names[ind].replace('A','dA') |
1822 | | |
1823 | | else: |
1824 | | |
1825 | | name = names[ind] |
1826 | | |
1827 | | if name in sigDict: |
1828 | | |
1829 | | atomsSig[str(i)+':'+str(ind)] = sigDict[name] |
1830 | | |
1831 | | if at[9] == 'I': |
1832 | | |
1833 | | at[10] = parmDict[names[10]] |
1834 | | |
1835 | | if names[10] in sigDict: |
1836 | | |
1837 | | atomsSig[str(i)+':10'] = sigDict[names[10]] |
1838 | | |
1839 | | else: |
1840 | | |
1841 | | for ind in [11,12,13,14,15,16]: |
1842 | | |
1843 | | at[ind] = parmDict[names[ind]] |
1844 | | |
1845 | | if names[ind] in sigDict: |
1846 | | |
1847 | | atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]] |
1848 | | |
1849 | | PrintAtomsAndSig(General,Atoms,atomsSig) |
1850 | | |
1851 | | |
1852 | | |
1853 | | if 'SH Texture' in General: |
1854 | | |
1855 | | textureData = General['SH Texture'] |
1856 | | |
1857 | | if textureData['Order']: |
1858 | | |
1859 | | SHtextureSig = {} |
1860 | | |
1861 | | for name in ['omega','chi','phi']: |
1862 | | |
1863 | | aname = pfx+'SH '+name |
1864 | | |
1865 | | textureData['Sample '+name][1] = parmDict[aname] |
1866 | | |
1867 | | if aname in sigDict: |
1868 | | |
1869 | | SHtextureSig['Sample '+name] = sigDict[aname] |
1870 | | |
1871 | | for name in textureData['SH Coeff'][1]: |
1872 | | |
1873 | | aname = pfx+name |
1874 | | |
1875 | | textureData['SH Coeff'][1][name] = parmDict[aname] |
1876 | | |
1877 | | if aname in sigDict: |
1878 | | |
1879 | | SHtextureSig[name] = sigDict[aname] |
1880 | | |
1881 | | PrintSHtextureAndSig(textureData,SHtextureSig) |
1882 | | |
1883 | | |
1884 | | |
1885 | | def GetHistogramPhaseData(Phases,Histograms,Print=True): |
1886 | | |
1887 | | |
1888 | | |
1889 | | def PrintSize(hapData): |
1890 | | |
1891 | | if hapData[0] in ['isotropic','uniaxial']: |
1892 | | |
1893 | | line = '\n Size model : %9s'%(hapData[0]) |
1894 | | |
1895 | | line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
1896 | | |
1897 | | if hapData[0] == 'uniaxial': |
1898 | | |
1899 | | line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
1900 | | |
1901 | | print line |
1902 | | |
1903 | | else: |
1904 | | |
1905 | | print '\n Size model : %s'%(hapData[0]) |
1906 | | |
1907 | | Snames = ['S11','S22','S33','S12','S13','S23'] |
1908 | | |
1909 | | ptlbls = ' names :' |
1910 | | |
1911 | | ptstr = ' values:' |
1912 | | |
1913 | | varstr = ' refine:' |
1914 | | |
1915 | | for i,name in enumerate(Snames): |
1916 | | |
1917 | | ptlbls += '%12s' % (name) |
1918 | | |
1919 | | ptstr += '%12.6f' % (hapData[4][i]) |
1920 | | |
1921 | | varstr += '%12s' % (str(hapData[5][i])) |
1922 | | |
1923 | | print ptlbls |
1924 | | |
1925 | | print ptstr |
1926 | | |
1927 | | print varstr |
1928 | | |
1929 | | |
1930 | | |
1931 | | def PrintMuStrain(hapData,SGData): |
1932 | | |
1933 | | if hapData[0] in ['isotropic','uniaxial']: |
1934 | | |
1935 | | line = '\n Mustrain model: %9s'%(hapData[0]) |
1936 | | |
1937 | | line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
1938 | | |
1939 | | if hapData[0] == 'uniaxial': |
1940 | | |
1941 | | line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
1942 | | |
1943 | | print line |
1944 | | |
1945 | | else: |
1946 | | |
1947 | | print '\n Mustrain model: %s'%(hapData[0]) |
1948 | | |
1949 | | Snames = G2spc.MustrainNames(SGData) |
1950 | | |
1951 | | ptlbls = ' names :' |
1952 | | |
1953 | | ptstr = ' values:' |
1954 | | |
1955 | | varstr = ' refine:' |
1956 | | |
1957 | | for i,name in enumerate(Snames): |
1958 | | |
1959 | | ptlbls += '%12s' % (name) |
1960 | | |
1961 | | ptstr += '%12.6f' % (hapData[4][i]) |
1962 | | |
1963 | | varstr += '%12s' % (str(hapData[5][i])) |
1964 | | |
1965 | | print ptlbls |
1966 | | |
1967 | | print ptstr |
1968 | | |
1969 | | print varstr |
1970 | | |
1971 | | |
1972 | | |
1973 | | def PrintHStrain(hapData,SGData): |
1974 | | |
1975 | | print '\n Hydrostatic/elastic strain: ' |
1976 | | |
1977 | | Hsnames = G2spc.HStrainNames(SGData) |
1978 | | |
1979 | | ptlbls = ' names :' |
1980 | | |
1981 | | ptstr = ' values:' |
1982 | | |
1983 | | varstr = ' refine:' |
1984 | | |
1985 | | for i,name in enumerate(Hsnames): |
1986 | | |
1987 | | ptlbls += '%12s' % (name) |
1988 | | |
1989 | | ptstr += '%12.6f' % (hapData[0][i]) |
1990 | | |
1991 | | varstr += '%12s' % (str(hapData[1][i])) |
1992 | | |
1993 | | print ptlbls |
1994 | | |
1995 | | print ptstr |
1996 | | |
1997 | | print varstr |
1998 | | |
1999 | | |
2000 | | |
2001 | | def PrintSHPO(hapData): |
2002 | | |
2003 | | print '\n Spherical harmonics preferred orientation: Order:' + \ |
2004 | | |
2005 | | str(hapData[4])+' Refine? '+str(hapData[2]) |
2006 | | |
2007 | | ptlbls = ' names :' |
2008 | | |
2009 | | ptstr = ' values:' |
2010 | | |
2011 | | for item in hapData[5]: |
2012 | | |
2013 | | ptlbls += '%12s'%(item) |
2014 | | |
2015 | | ptstr += '%12.3f'%(hapData[5][item]) |
2016 | | |
2017 | | print ptlbls |
2018 | | |
2019 | | print ptstr |
2020 | | |
2021 | | |
2022 | | |
2023 | | hapDict = {} |
2024 | | |
2025 | | hapVary = [] |
2026 | | |
2027 | | controlDict = {} |
2028 | | |
2029 | | poType = {} |
2030 | | |
2031 | | poAxes = {} |
2032 | | |
2033 | | spAxes = {} |
2034 | | |
2035 | | spType = {} |
2036 | | |
2037 | | |
2038 | | |
2039 | | for phase in Phases: |
2040 | | |
2041 | | HistoPhase = Phases[phase]['Histograms'] |
2042 | | |
2043 | | SGData = Phases[phase]['General']['SGData'] |
2044 | | |
2045 | | cell = Phases[phase]['General']['Cell'][1:7] |
2046 | | |
2047 | | A = G2lat.cell2A(cell) |
2048 | | |
2049 | | pId = Phases[phase]['pId'] |
2050 | | |
2051 | | histoList = HistoPhase.keys() |
2052 | | |
2053 | | histoList.sort() |
2054 | | |
2055 | | for histogram in histoList: |
2056 | | |
2057 | | try: |
2058 | | |
2059 | | Histogram = Histograms[histogram] |
2060 | | |
2061 | | except KeyError: |
2062 | | |
2063 | | #skip if histogram not included e.g. in a sequential refinement |
2064 | | |
2065 | | continue |
2066 | | |
2067 | | hapData = HistoPhase[histogram] |
2068 | | |
2069 | | hId = Histogram['hId'] |
2070 | | |
2071 | | limits = Histogram['Limits'][1] |
2072 | | |
2073 | | inst = Histogram['Instrument Parameters'] |
2074 | | |
2075 | | inst = dict(zip(inst[3],inst[1])) |
2076 | | |
2077 | | Zero = inst['Zero'] |
2078 | | |
2079 | | if 'C' in inst['Type']: |
2080 | | |
2081 | | try: |
2082 | | |
2083 | | wave = inst['Lam'] |
2084 | | |
2085 | | except KeyError: |
2086 | | |
2087 | | wave = inst['Lam1'] |
2088 | | |
2089 | | dmin = wave/(2.0*sind(limits[1]/2.0)) |
2090 | | |
2091 | | pfx = str(pId)+':'+str(hId)+':' |
2092 | | |
2093 | | for item in ['Scale','Extinction']: |
2094 | | |
2095 | | hapDict[pfx+item] = hapData[item][0] |
2096 | | |
2097 | | if hapData[item][1]: |
2098 | | |
2099 | | hapVary.append(pfx+item) |
2100 | | |
2101 | | names = G2spc.HStrainNames(SGData) |
2102 | | |
2103 | | for i,name in enumerate(names): |
2104 | | |
2105 | | hapDict[pfx+name] = hapData['HStrain'][0][i] |
2106 | | |
2107 | | if hapData['HStrain'][1][i]: |
2108 | | |
2109 | | hapVary.append(pfx+name) |
2110 | | |
2111 | | controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0] |
2112 | | |
2113 | | if hapData['Pref.Ori.'][0] == 'MD': |
2114 | | |
2115 | | hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1] |
2116 | | |
2117 | | controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3] |
2118 | | |
2119 | | if hapData['Pref.Ori.'][2]: |
2120 | | |
2121 | | hapVary.append(pfx+'MD') |
2122 | | |
2123 | | else: #'SH' spherical harmonics |
2124 | | |
2125 | | controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4] |
2126 | | |
2127 | | controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5]) |
2128 | | |
2129 | | for item in hapData['Pref.Ori.'][5]: |
2130 | | |
2131 | | hapDict[pfx+item] = hapData['Pref.Ori.'][5][item] |
2132 | | |
2133 | | if hapData['Pref.Ori.'][2]: |
2134 | | |
2135 | | hapVary.append(pfx+item) |
2136 | | |
2137 | | for item in ['Mustrain','Size']: |
2138 | | |
2139 | | controlDict[pfx+item+'Type'] = hapData[item][0] |
2140 | | |
2141 | | if hapData[item][0] in ['isotropic','uniaxial']: |
2142 | | |
2143 | | hapDict[pfx+item+':i'] = hapData[item][1][0] |
2144 | | |
2145 | | if hapData[item][2][0]: |
2146 | | |
2147 | | hapVary.append(pfx+item+':i') |
2148 | | |
2149 | | if hapData[item][0] == 'uniaxial': |
2150 | | |
2151 | | controlDict[pfx+item+'Axis'] = hapData[item][3] |
2152 | | |
2153 | | hapDict[pfx+item+':a'] = hapData[item][1][1] |
2154 | | |
2155 | | if hapData[item][2][1]: |
2156 | | |
2157 | | hapVary.append(pfx+item+':a') |
2158 | | |
2159 | | else: #generalized for mustrain or ellipsoidal for size |
2160 | | |
2161 | | if item == 'Mustrain': |
2162 | | |
2163 | | names = G2spc.MustrainNames(SGData) |
2164 | | |
2165 | | pwrs = [] |
2166 | | |
2167 | | for name in names: |
2168 | | |
2169 | | h,k,l = name[1:] |
2170 | | |
2171 | | pwrs.append([int(h),int(k),int(l)]) |
2172 | | |
2173 | | controlDict[pfx+'MuPwrs'] = pwrs |
2174 | | |
2175 | | for i in range(len(hapData[item][4])): |
2176 | | |
2177 | | sfx = ':'+str(i) |
2178 | | |
2179 | | hapDict[pfx+item+sfx] = hapData[item][4][i] |
2180 | | |
2181 | | if hapData[item][5][i]: |
2182 | | |
2183 | | hapVary.append(pfx+item+sfx) |
2184 | | |
2185 | | |
2186 | | |
2187 | | if Print: |
2188 | | |
2189 | | print '\n Phase: ',phase,' in histogram: ',histogram |
2190 | | |
2191 | | print 135*'-' |
2192 | | |
2193 | | print ' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1] |
2194 | | |
2195 | | print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1] |
2196 | | |
2197 | | if hapData['Pref.Ori.'][0] == 'MD': |
2198 | | |
2199 | | Ax = hapData['Pref.Ori.'][3] |
2200 | | |
2201 | | print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \ |
2202 | | |
2203 | | ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2]) |
2204 | | |
2205 | | else: #'SH' for spherical harmonics |
2206 | | |
2207 | | PrintSHPO(hapData['Pref.Ori.']) |
2208 | | |
2209 | | PrintSize(hapData['Size']) |
2210 | | |
2211 | | PrintMuStrain(hapData['Mustrain'],SGData) |
2212 | | |
2213 | | PrintHStrain(hapData['HStrain'],SGData) |
2214 | | |
2215 | | HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A)) |
2216 | | |
2217 | | refList = [] |
2218 | | |
2219 | | for h,k,l,d in HKLd: |
2220 | | |
2221 | | ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData) |
2222 | | |
2223 | | if ext: |
2224 | | |
2225 | | continue |
2226 | | |
2227 | | if 'C' in inst['Type']: |
2228 | | |
2229 | | pos = 2.0*asind(wave/(2.0*d))+Zero |
2230 | | |
2231 | | if limits[0] < pos < limits[1]: |
2232 | | |
2233 | | refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0]) |
2234 | | |
2235 | | else: |
2236 | | |
2237 | | raise ValueError |
2238 | | |
2239 | | Histogram['Reflection Lists'][phase] = refList |
2240 | | |
2241 | | return hapVary,hapDict,controlDict |
2242 | | |
2243 | | |
2244 | | |
2245 | | def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True): |
2246 | | |
2247 | | |
2248 | | |
2249 | | def PrintSizeAndSig(hapData,sizeSig): |
2250 | | |
2251 | | line = '\n Size model: %9s'%(hapData[0]) |
2252 | | |
2253 | | refine = False |
2254 | | |
2255 | | if hapData[0] in ['isotropic','uniaxial']: |
2256 | | |
2257 | | line += ' equatorial:%12.3f'%(hapData[1][0]) |
2258 | | |
2259 | | if sizeSig[0][0]: |
2260 | | |
2261 | | line += ', sig: %8.3f'%(sizeSig[0][0]) |
2262 | | |
2263 | | refine = True |
2264 | | |
2265 | | if hapData[0] == 'uniaxial': |
2266 | | |
2267 | | line += ' axial:%12.3f'%(hapData[1][1]) |
2268 | | |
2269 | | if sizeSig[0][1]: |
2270 | | |
2271 | | refine = True |
2272 | | |
2273 | | line += ', sig: %8.3f'%(sizeSig[0][1]) |
2274 | | |
2275 | | if refine: |
2276 | | |
2277 | | print line |
2278 | | |
2279 | | else: |
2280 | | |
2281 | | Snames = ['S11','S22','S33','S12','S13','S23'] |
2282 | | |
2283 | | ptlbls = ' name :' |
2284 | | |
2285 | | ptstr = ' value :' |
2286 | | |
2287 | | sigstr = ' sig :' |
2288 | | |
2289 | | for i,name in enumerate(Snames): |
2290 | | |
2291 | | ptlbls += '%12s' % (name) |
2292 | | |
2293 | | ptstr += '%12.6f' % (hapData[4][i]) |
2294 | | |
2295 | | if sizeSig[1][i]: |
2296 | | |
2297 | | refine = True |
2298 | | |
2299 | | sigstr += '%12.6f' % (sizeSig[1][i]) |
2300 | | |
2301 | | else: |
2302 | | |
2303 | | sigstr += 12*' ' |
2304 | | |
2305 | | if refine: |
2306 | | |
2307 | | print line |
2308 | | |
2309 | | print ptlbls |
2310 | | |
2311 | | print ptstr |
2312 | | |
2313 | | print sigstr |
2314 | | |
2315 | | |
2316 | | |
2317 | | def PrintMuStrainAndSig(hapData,mustrainSig,SGData): |
2318 | | |
2319 | | line = '\n Mustrain model: %9s'%(hapData[0]) |
2320 | | |
2321 | | refine = False |
2322 | | |
2323 | | if hapData[0] in ['isotropic','uniaxial']: |
2324 | | |
2325 | | line += ' equatorial:%12.1f'%(hapData[1][0]) |
2326 | | |
2327 | | if mustrainSig[0][0]: |
2328 | | |
2329 | | line += ', sig: %8.1f'%(mustrainSig[0][0]) |
2330 | | |
2331 | | refine = True |
2332 | | |
2333 | | if hapData[0] == 'uniaxial': |
2334 | | |
2335 | | line += ' axial:%12.1f'%(hapData[1][1]) |
2336 | | |
2337 | | if mustrainSig[0][1]: |
2338 | | |
2339 | | line += ', sig: %8.1f'%(mustrainSig[0][1]) |
2340 | | |
2341 | | if refine: |
2342 | | |
2343 | | print line |
2344 | | |
2345 | | else: |
2346 | | |
2347 | | Snames = G2spc.MustrainNames(SGData) |
2348 | | |
2349 | | ptlbls = ' name :' |
2350 | | |
2351 | | ptstr = ' value :' |
2352 | | |
2353 | | sigstr = ' sig :' |
2354 | | |
2355 | | for i,name in enumerate(Snames): |
2356 | | |
2357 | | ptlbls += '%12s' % (name) |
2358 | | |
2359 | | ptstr += '%12.6f' % (hapData[4][i]) |
2360 | | |
2361 | | if mustrainSig[1][i]: |
2362 | | |
2363 | | refine = True |
2364 | | |
2365 | | sigstr += '%12.6f' % (mustrainSig[1][i]) |
2366 | | |
2367 | | else: |
2368 | | |
2369 | | sigstr += 12*' ' |
2370 | | |
2371 | | if refine: |
2372 | | |
2373 | | print line |
2374 | | |
2375 | | print ptlbls |
2376 | | |
2377 | | print ptstr |
2378 | | |
2379 | | print sigstr |
2380 | | |
2381 | | |
2382 | | |
2383 | | def PrintHStrainAndSig(hapData,strainSig,SGData): |
2384 | | |
2385 | | Hsnames = G2spc.HStrainNames(SGData) |
2386 | | |
2387 | | ptlbls = ' name :' |
2388 | | |
2389 | | ptstr = ' value :' |
2390 | | |
2391 | | sigstr = ' sig :' |
2392 | | |
2393 | | refine = False |
2394 | | |
2395 | | for i,name in enumerate(Hsnames): |
2396 | | |
2397 | | ptlbls += '%12s' % (name) |
2398 | | |
2399 | | ptstr += '%12.6g' % (hapData[0][i]) |
2400 | | |
2401 | | if name in strainSig: |
2402 | | |
2403 | | refine = True |
2404 | | |
2405 | | sigstr += '%12.6g' % (strainSig[name]) |
2406 | | |
2407 | | else: |
2408 | | |
2409 | | sigstr += 12*' ' |
2410 | | |
2411 | | if refine: |
2412 | | |
2413 | | print '\n Hydrostatic/elastic strain: ' |
2414 | | |
2415 | | print ptlbls |
2416 | | |
2417 | | print ptstr |
2418 | | |
2419 | | print sigstr |
2420 | | |
2421 | | |
2422 | | |
2423 | | def PrintSHPOAndSig(hapData,POsig): |
2424 | | |
2425 | | print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4]) |
2426 | | |
2427 | | ptlbls = ' names :' |
2428 | | |
2429 | | ptstr = ' values:' |
2430 | | |
2431 | | sigstr = ' sig :' |
2432 | | |
2433 | | for item in hapData[5]: |
2434 | | |
2435 | | ptlbls += '%12s'%(item) |
2436 | | |
2437 | | ptstr += '%12.3f'%(hapData[5][item]) |
2438 | | |
2439 | | if item in POsig: |
2440 | | |
2441 | | sigstr += '%12.3f'%(POsig[item]) |
2442 | | |
2443 | | else: |
2444 | | |
2445 | | sigstr += 12*' ' |
2446 | | |
2447 | | print ptlbls |
2448 | | |
2449 | | print ptstr |
2450 | | |
2451 | | print sigstr |
2452 | | |
2453 | | |
2454 | | |
2455 | | for phase in Phases: |
2456 | | |
2457 | | HistoPhase = Phases[phase]['Histograms'] |
2458 | | |
2459 | | SGData = Phases[phase]['General']['SGData'] |
2460 | | |
2461 | | pId = Phases[phase]['pId'] |
2462 | | |
2463 | | histoList = HistoPhase.keys() |
2464 | | |
2465 | | histoList.sort() |
2466 | | |
2467 | | for histogram in histoList: |
2468 | | |
2469 | | try: |
2470 | | |
2471 | | Histogram = Histograms[histogram] |
2472 | | |
2473 | | except KeyError: |
2474 | | |
2475 | | #skip if histogram not included e.g. in a sequential refinement |
2476 | | |
2477 | | continue |
2478 | | |
2479 | | print '\n Phase: ',phase,' in histogram: ',histogram |
2480 | | |
2481 | | print 130*'-' |
2482 | | |
2483 | | hapData = HistoPhase[histogram] |
2484 | | |
2485 | | hId = Histogram['hId'] |
2486 | | |
2487 | | pfx = str(pId)+':'+str(hId)+':' |
2488 | | |
2489 | | print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \ |
2490 | | |
2491 | | %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref']) |
2492 | | |
2493 | | |
2494 | | |
2495 | | PhFrExtPOSig = {} |
2496 | | |
2497 | | for item in ['Scale','Extinction']: |
2498 | | |
2499 | | hapData[item][0] = parmDict[pfx+item] |
2500 | | |
2501 | | if pfx+item in sigDict: |
2502 | | |
2503 | | PhFrExtPOSig[item] = sigDict[pfx+item] |
2504 | | |
2505 | | if hapData['Pref.Ori.'][0] == 'MD': |
2506 | | |
2507 | | hapData['Pref.Ori.'][1] = parmDict[pfx+'MD'] |
2508 | | |
2509 | | if pfx+'MD' in sigDict: |
2510 | | |
2511 | | PhFrExtPOSig['MD'] = sigDict[pfx+'MD'] |
2512 | | |
2513 | | else: #'SH' spherical harmonics |
2514 | | |
2515 | | for item in hapData['Pref.Ori.'][5]: |
2516 | | |
2517 | | hapData['Pref.Ori.'][5][item] = parmDict[pfx+item] |
2518 | | |
2519 | | if pfx+item in sigDict: |
2520 | | |
2521 | | PhFrExtPOSig[item] = sigDict[pfx+item] |
2522 | | |
2523 | | if Print: |
2524 | | |
2525 | | if 'Scale' in PhFrExtPOSig: |
2526 | | |
2527 | | print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale']) |
2528 | | |
2529 | | if 'Extinction' in PhFrExtPOSig: |
2530 | | |
2531 | | print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction']) |
2532 | | |
2533 | | if hapData['Pref.Ori.'][0] == 'MD': |
2534 | | |
2535 | | if 'MD' in PhFrExtPOSig: |
2536 | | |
2537 | | print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD']) |
2538 | | |
2539 | | else: |
2540 | | |
2541 | | PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig) |
2542 | | |
2543 | | SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]], |
2544 | | |
2545 | | 'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]], |
2546 | | |
2547 | | 'HStrain':{}} |
2548 | | |
2549 | | for item in ['Mustrain','Size']: |
2550 | | |
2551 | | if hapData[item][0] in ['isotropic','uniaxial']: |
2552 | | |
2553 | | hapData[item][1][0] = parmDict[pfx+item+':i'] |
2554 | | |
2555 | | if item == 'Size': |
2556 | | |
2557 | | hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0])) |
2558 | | |
2559 | | if pfx+item+':i' in sigDict: |
2560 | | |
2561 | | SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i'] |
2562 | | |
2563 | | if hapData[item][0] == 'uniaxial': |
2564 | | |
2565 | | hapData[item][1][1] = parmDict[pfx+item+':a'] |
2566 | | |
2567 | | if item == 'Size': |
2568 | | |
2569 | | hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1])) |
2570 | | |
2571 | | if pfx+item+':a' in sigDict: |
2572 | | |
2573 | | SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a'] |
2574 | | |
2575 | | else: #generalized for mustrain or ellipsoidal for size |
2576 | | |
2577 | | for i in range(len(hapData[item][4])): |
2578 | | |
2579 | | sfx = ':'+str(i) |
2580 | | |
2581 | | hapData[item][4][i] = parmDict[pfx+item+sfx] |
2582 | | |
2583 | | if pfx+item+sfx in sigDict: |
2584 | | |
2585 | | SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx] |
2586 | | |
2587 | | names = G2spc.HStrainNames(SGData) |
2588 | | |
2589 | | for i,name in enumerate(names): |
2590 | | |
2591 | | hapData['HStrain'][0][i] = parmDict[pfx+name] |
2592 | | |
2593 | | if pfx+name in sigDict: |
2594 | | |
2595 | | SizeMuStrSig['HStrain'][name] = sigDict[pfx+name] |
2596 | | |
2597 | | if Print: |
2598 | | |
2599 | | PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size']) |
2600 | | |
2601 | | PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData) |
2602 | | |
2603 | | PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData) |
2604 | | |
2605 | | |
2606 | | |
2607 | | def GetHistogramData(Histograms,Print=True): |
2608 | | |
2609 | | |
2610 | | |
2611 | | def GetBackgroundParms(hId,Background): |
2612 | | |
2613 | | Back = Background[0] |
2614 | | |
2615 | | Debye = Background[1] |
2616 | | |
2617 | | bakType,bakFlag = Back[:2] |
2618 | | |
2619 | | backVals = Back[3:] |
2620 | | |
2621 | | backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))] |
2622 | | |
2623 | | backDict = dict(zip(backNames,backVals)) |
2624 | | |
2625 | | backVary = [] |
2626 | | |
2627 | | if bakFlag: |
2628 | | |
2629 | | backVary = backNames |
2630 | | |
2631 | | backDict[':'+str(hId)+':nDebye'] = Debye['nDebye'] |
2632 | | |
2633 | | debyeDict = {} |
2634 | | |
2635 | | debyeList = [] |
2636 | | |
2637 | | for i in range(Debye['nDebye']): |
2638 | | |
2639 | | debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)] |
2640 | | |
2641 | | debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2]))) |
2642 | | |
2643 | | debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2]) |
2644 | | |
2645 | | debyeVary = [] |
2646 | | |
2647 | | for item in debyeList: |
2648 | | |
2649 | | if item[1]: |
2650 | | |
2651 | | debyeVary.append(item[0]) |
2652 | | |
2653 | | backDict.update(debyeDict) |
2654 | | |
2655 | | backVary += debyeVary |
2656 | | |
2657 | | return bakType,backDict,backVary |
2658 | | |
2659 | | |
2660 | | |
2661 | | def GetInstParms(hId,Inst): |
2662 | | |
2663 | | insVals,insFlags,insNames = Inst[1:4] |
2664 | | |
2665 | | dataType = insVals[0] |
2666 | | |
2667 | | instDict = {} |
2668 | | |
2669 | | insVary = [] |
2670 | | |
2671 | | pfx = ':'+str(hId)+':' |
2672 | | |
2673 | | for i,flag in enumerate(insFlags): |
2674 | | |
2675 | | insName = pfx+insNames[i] |
2676 | | |
2677 | | instDict[insName] = insVals[i] |
2678 | | |
2679 | | if flag: |
2680 | | |
2681 | | insVary.append(insName) |
2682 | | |
2683 | | instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001) |
2684 | | |
2685 | | instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001) |
2686 | | |
2687 | | instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005) |
2688 | | |
2689 | | return dataType,instDict,insVary |
2690 | | |
2691 | | |
2692 | | |
2693 | | def GetSampleParms(hId,Sample): |
2694 | | |
2695 | | sampVary = [] |
2696 | | |
2697 | | hfx = ':'+str(hId)+':' |
2698 | | |
2699 | | sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'], |
2700 | | |
2701 | | hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']} |
2702 | | |
2703 | | Type = Sample['Type'] |
2704 | | |
2705 | | if 'Bragg' in Type: #Bragg-Brentano |
2706 | | |
2707 | | for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering? |
2708 | | |
2709 | | sampDict[hfx+item] = Sample[item][0] |
2710 | | |
2711 | | if Sample[item][1]: |
2712 | | |
2713 | | sampVary.append(hfx+item) |
2714 | | |
2715 | | elif 'Debye' in Type: #Debye-Scherrer |
2716 | | |
2717 | | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
2718 | | |
2719 | | sampDict[hfx+item] = Sample[item][0] |
2720 | | |
2721 | | if Sample[item][1]: |
2722 | | |
2723 | | sampVary.append(hfx+item) |
2724 | | |
2725 | | return Type,sampDict,sampVary |
2726 | | |
2727 | | |
2728 | | |
2729 | | def PrintBackground(Background): |
2730 | | |
2731 | | Back = Background[0] |
2732 | | |
2733 | | Debye = Background[1] |
2734 | | |
2735 | | print '\n Background function: ',Back[0],' Refine?',bool(Back[1]) |
2736 | | |
2737 | | line = ' Coefficients: ' |
2738 | | |
2739 | | for i,back in enumerate(Back[3:]): |
2740 | | |
2741 | | line += '%10.3f'%(back) |
2742 | | |
2743 | | if i and not i%10: |
2744 | | |
2745 | | line += '\n'+15*' ' |
2746 | | |
2747 | | print line |
2748 | | |
2749 | | if Debye['nDebye']: |
2750 | | |
2751 | | print '\n Debye diffuse scattering coefficients' |
2752 | | |
2753 | | parms = ['DebyeA','DebyeR','DebyeU'] |
2754 | | |
2755 | | line = ' names :' |
2756 | | |
2757 | | for parm in parms: |
2758 | | |
2759 | | line += '%16s'%(parm) |
2760 | | |
2761 | | print line |
2762 | | |
2763 | | for j,term in enumerate(Debye['debyeTerms']): |
2764 | | |
2765 | | line = ' term'+'%2d'%(j)+':' |
2766 | | |
2767 | | for i in range(3): |
2768 | | |
2769 | | line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1])) |
2770 | | |
2771 | | print line |
2772 | | |
2773 | | |
2774 | | |
2775 | | def PrintInstParms(Inst): |
2776 | | |
2777 | | print '\n Instrument Parameters:' |
2778 | | |
2779 | | ptlbls = ' name :' |
2780 | | |
2781 | | ptstr = ' value :' |
2782 | | |
2783 | | varstr = ' refine:' |
2784 | | |
2785 | | instNames = Inst[3][1:] |
2786 | | |
2787 | | for i,name in enumerate(instNames): |
2788 | | |
2789 | | ptlbls += '%12s' % (name) |
2790 | | |
2791 | | ptstr += '%12.6f' % (Inst[1][i+1]) |
2792 | | |
2793 | | if name in ['Lam1','Lam2','Azimuth']: |
2794 | | |
2795 | | varstr += 12*' ' |
2796 | | |
2797 | | else: |
2798 | | |
2799 | | varstr += '%12s' % (str(bool(Inst[2][i+1]))) |
2800 | | |
2801 | | print ptlbls |
2802 | | |
2803 | | print ptstr |
2804 | | |
2805 | | print varstr |
2806 | | |
2807 | | |
2808 | | |
2809 | | def PrintSampleParms(Sample): |
2810 | | |
2811 | | print '\n Sample Parameters:' |
2812 | | |
2813 | | print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \ |
2814 | | |
2815 | | (Sample['Omega'],Sample['Chi'],Sample['Phi']) |
2816 | | |
2817 | | ptlbls = ' name :' |
2818 | | |
2819 | | ptstr = ' value :' |
2820 | | |
2821 | | varstr = ' refine:' |
2822 | | |
2823 | | if 'Bragg' in Sample['Type']: |
2824 | | |
2825 | | for item in ['Scale','Shift','Transparency']: |
2826 | | |
2827 | | ptlbls += '%14s'%(item) |
2828 | | |
2829 | | ptstr += '%14.4f'%(Sample[item][0]) |
2830 | | |
2831 | | varstr += '%14s'%(str(bool(Sample[item][1]))) |
2832 | | |
2833 | | |
2834 | | |
2835 | | elif 'Debye' in Type: #Debye-Scherrer |
2836 | | |
2837 | | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
2838 | | |
2839 | | ptlbls += '%14s'%(item) |
2840 | | |
2841 | | ptstr += '%14.4f'%(Sample[item][0]) |
2842 | | |
2843 | | varstr += '%14s'%(str(bool(Sample[item][1]))) |
2844 | | |
2845 | | |
2846 | | |
2847 | | print ptlbls |
2848 | | |
2849 | | print ptstr |
2850 | | |
2851 | | print varstr |
2852 | | |
2853 | | |
2854 | | |
2855 | | |
2856 | | |
2857 | | histDict = {} |
2858 | | |
2859 | | histVary = [] |
2860 | | |
2861 | | controlDict = {} |
2862 | | |
2863 | | histoList = Histograms.keys() |
2864 | | |
2865 | | histoList.sort() |
2866 | | |
2867 | | for histogram in histoList: |
2868 | | |
2869 | | Histogram = Histograms[histogram] |
2870 | | |
2871 | | hId = Histogram['hId'] |
2872 | | |
2873 | | pfx = ':'+str(hId)+':' |
2874 | | |
2875 | | controlDict[pfx+'Limits'] = Histogram['Limits'][1] |
2876 | | |
2877 | | |
2878 | | |
2879 | | Background = Histogram['Background'] |
2880 | | |
2881 | | Type,bakDict,bakVary = GetBackgroundParms(hId,Background) |
2882 | | |
2883 | | controlDict[pfx+'bakType'] = Type |
2884 | | |
2885 | | histDict.update(bakDict) |
2886 | | |
2887 | | histVary += bakVary |
2888 | | |
2889 | | |
2890 | | |
2891 | | Inst = Histogram['Instrument Parameters'] |
2892 | | |
2893 | | Type,instDict,insVary = GetInstParms(hId,Inst) |
2894 | | |
2895 | | controlDict[pfx+'histType'] = Type |
2896 | | |
2897 | | if pfx+'Lam1' in instDict: |
2898 | | |
2899 | | controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1'] |
2900 | | |
2901 | | else: |
2902 | | |
2903 | | controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam'] |
2904 | | |
2905 | | histDict.update(instDict) |
2906 | | |
2907 | | histVary += insVary |
2908 | | |
2909 | | |
2910 | | |
2911 | | Sample = Histogram['Sample Parameters'] |
2912 | | |
2913 | | Type,sampDict,sampVary = GetSampleParms(hId,Sample) |
2914 | | |
2915 | | controlDict[pfx+'instType'] = Type |
2916 | | |
2917 | | histDict.update(sampDict) |
2918 | | |
2919 | | histVary += sampVary |
2920 | | |
2921 | | |
2922 | | |
2923 | | if Print: |
2924 | | |
2925 | | print '\n Histogram: ',histogram,' histogram Id: ',hId |
2926 | | |
2927 | | print 135*'-' |
2928 | | |
2929 | | Units = {'C':' deg','T':' msec'} |
2930 | | |
2931 | | units = Units[controlDict[pfx+'histType'][2]] |
2932 | | |
2933 | | Limits = controlDict[pfx+'Limits'] |
2934 | | |
2935 | | print ' Instrument type: ',Sample['Type'] |
2936 | | |
2937 | | print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units) |
2938 | | |
2939 | | PrintSampleParms(Sample) |
2940 | | |
2941 | | PrintInstParms(Inst) |
2942 | | |
2943 | | PrintBackground(Background) |
2944 | | |
2945 | | |
2946 | | |
2947 | | return histVary,histDict,controlDict |
2948 | | |
2949 | | |
2950 | | |
2951 | | def SetHistogramData(parmDict,sigDict,Histograms,Print=True): |
2952 | | |
2953 | | |
2954 | | |
2955 | | def SetBackgroundParms(pfx,Background,parmDict,sigDict): |
2956 | | |
2957 | | Back = Background[0] |
2958 | | |
2959 | | Debye = Background[1] |
2960 | | |
2961 | | lenBack = len(Back[3:]) |
2962 | | |
2963 | | backSig = [0 for i in range(lenBack+3*Debye['nDebye'])] |
2964 | | |
2965 | | for i in range(lenBack): |
2966 | | |
2967 | | Back[3+i] = parmDict[pfx+'Back:'+str(i)] |
2968 | | |
2969 | | if pfx+'Back:'+str(i) in sigDict: |
2970 | | |
2971 | | backSig[i] = sigDict[pfx+'Back:'+str(i)] |
2972 | | |
2973 | | if Debye['nDebye']: |
2974 | | |
2975 | | for i in range(Debye['nDebye']): |
2976 | | |
2977 | | names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)] |
2978 | | |
2979 | | for j,name in enumerate(names): |
2980 | | |
2981 | | Debye['debyeTerms'][i][2*j] = parmDict[name] |
2982 | | |
2983 | | if name in sigDict: |