239 | | def getBackupName(GPXfile,makeBack): #recovered old one |
240 | | GPXpath,GPXname = ospath.split(GPXfile) |
241 | | if GPXpath == '': GPXpath = '.' |
242 | | Name = ospath.splitext(GPXname)[0] |
243 | | files = os.listdir(GPXpath) |
244 | | last = 0 |
245 | | for name in files: |
246 | | name = name.split('.') |
247 | | if len(name) == 3 and name[0] == Name and 'bak' in name[1]: |
248 | | if makeBack: |
249 | | last = max(last,int(name[1].strip('bak'))+1) |
| 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 | |
| 507 | input: |
| 508 | |
| 509 | GPXfile = .gpx full file name |
| 510 | |
| 511 | Histograms = dictionary of histograms as {name:data,...} |
| 512 | |
| 513 | Phases = dictionary of phases that use histograms |
| 514 | |
| 515 | CovData = dictionary of refined variables, varyList, & covariance matrix |
| 516 | |
| 517 | makeBack = True if new backup of .gpx file is to be made; else use the last one made |
| 518 | |
| 519 | ''' |
| 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') |
| 534 | |
| 535 | while True: |
| 536 | |
| 537 | try: |
| 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 | |
| 701 | try: |
| 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): |
| 956 | |
| 957 | |
| 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 | |
255 | | def GPXBackup(GPXfile,makeBack=True): |
256 | | import distutils.file_util as dfu |
257 | | GPXback = getBackupName(GPXfile,makeBack) |
258 | | dfu.copy_file(GPXfile,GPXback) |
259 | | return GPXback |
260 | | |
261 | | def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True): |
262 | | ''' Updates gpxfile from all histograms that are found in any phase |
263 | | and any phase that used a histogram |
| 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: |
| 2984 | |
| 2985 | backSig[lenBack+3*i+j] = sigDict[name] |
| 2986 | |
| 2987 | return backSig |
| 2988 | |
| 2989 | |
| 2990 | |
| 2991 | def SetInstParms(pfx,Inst,parmDict,sigDict): |
| 2992 | |
| 2993 | insVals,insFlags,insNames = Inst[1:4] |
| 2994 | |
| 2995 | instSig = [0 for i in range(len(insVals))] |
| 2996 | |
| 2997 | for i,flag in enumerate(insFlags): |
| 2998 | |
| 2999 | insName = pfx+insNames[i] |
| 3000 | |
| 3001 | insVals[i] = parmDict[insName] |
| 3002 | |
| 3003 | if insName in sigDict: |
| 3004 | |
| 3005 | instSig[i] = sigDict[insName] |
| 3006 | |
| 3007 | return instSig |
| 3008 | |
| 3009 | |
| 3010 | |
| 3011 | def SetSampleParms(pfx,Sample,parmDict,sigDict): |
| 3012 | |
| 3013 | if 'Bragg' in Sample['Type']: #Bragg-Brentano |
| 3014 | |
| 3015 | sampSig = [0 for i in range(3)] |
| 3016 | |
| 3017 | for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering? |
| 3018 | |
| 3019 | Sample[item][0] = parmDict[pfx+item] |
| 3020 | |
| 3021 | if pfx+item in sigDict: |
| 3022 | |
| 3023 | sampSig[i] = sigDict[pfx+item] |
| 3024 | |
| 3025 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
| 3026 | |
| 3027 | sampSig = [0 for i in range(4)] |
| 3028 | |
| 3029 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
| 3030 | |
| 3031 | Sample[item][0] = parmDict[pfx+item] |
| 3032 | |
| 3033 | if pfx+item in sigDict: |
| 3034 | |
| 3035 | sampSig[i] = sigDict[pfx+item] |
| 3036 | |
| 3037 | return sampSig |
| 3038 | |
| 3039 | |
| 3040 | |
| 3041 | def PrintBackgroundSig(Background,backSig): |
| 3042 | |
| 3043 | Back = Background[0] |
| 3044 | |
| 3045 | Debye = Background[1] |
| 3046 | |
| 3047 | lenBack = len(Back[3:]) |
| 3048 | |
| 3049 | valstr = ' value : ' |
| 3050 | |
| 3051 | sigstr = ' sig : ' |
| 3052 | |
| 3053 | refine = False |
| 3054 | |
| 3055 | for i,back in enumerate(Back[3:]): |
| 3056 | |
| 3057 | valstr += '%10.4g'%(back) |
| 3058 | |
| 3059 | if Back[1]: |
| 3060 | |
| 3061 | refine = True |
| 3062 | |
| 3063 | sigstr += '%10.4g'%(backSig[i]) |
| 3064 | |
| 3065 | else: |
| 3066 | |
| 3067 | sigstr += 10*' ' |
| 3068 | |
| 3069 | if refine: |
| 3070 | |
| 3071 | print '\n Background function: ',Back[0] |
| 3072 | |
| 3073 | print valstr |
| 3074 | |
| 3075 | print sigstr |
| 3076 | |
| 3077 | if Debye['nDebye']: |
| 3078 | |
| 3079 | ifAny = False |
| 3080 | |
| 3081 | ptfmt = "%12.5f" |
| 3082 | |
| 3083 | names = ' names :' |
| 3084 | |
| 3085 | ptstr = ' values:' |
| 3086 | |
| 3087 | sigstr = ' esds :' |
| 3088 | |
| 3089 | for item in sigDict: |
| 3090 | |
| 3091 | if 'Debye' in item: |
| 3092 | |
| 3093 | ifAny = True |
| 3094 | |
| 3095 | names += '%12s'%(item) |
| 3096 | |
| 3097 | ptstr += ptfmt%(parmDict[item]) |
| 3098 | |
| 3099 | sigstr += ptfmt%(sigDict[item]) |
| 3100 | |
| 3101 | if ifAny: |
| 3102 | |
| 3103 | print '\n Debye diffuse scattering coefficients' |
| 3104 | |
| 3105 | print names |
| 3106 | |
| 3107 | print ptstr |
| 3108 | |
| 3109 | print sigstr |
| 3110 | |
| 3111 | |
| 3112 | |
| 3113 | def PrintInstParmsSig(Inst,instSig): |
| 3114 | |
| 3115 | ptlbls = ' names :' |
| 3116 | |
| 3117 | ptstr = ' value :' |
| 3118 | |
| 3119 | sigstr = ' sig :' |
| 3120 | |
| 3121 | instNames = Inst[3][1:] |
| 3122 | |
| 3123 | refine = False |
| 3124 | |
| 3125 | for i,name in enumerate(instNames): |
| 3126 | |
| 3127 | ptlbls += '%12s' % (name) |
| 3128 | |
| 3129 | ptstr += '%12.6f' % (Inst[1][i+1]) |
| 3130 | |
| 3131 | if instSig[i+1]: |
| 3132 | |
| 3133 | refine = True |
| 3134 | |
| 3135 | sigstr += '%12.6f' % (instSig[i+1]) |
| 3136 | |
| 3137 | else: |
| 3138 | |
| 3139 | sigstr += 12*' ' |
| 3140 | |
| 3141 | if refine: |
| 3142 | |
| 3143 | print '\n Instrument Parameters:' |
| 3144 | |
| 3145 | print ptlbls |
| 3146 | |
| 3147 | print ptstr |
| 3148 | |
| 3149 | print sigstr |
| 3150 | |
| 3151 | |
| 3152 | |
| 3153 | def PrintSampleParmsSig(Sample,sampleSig): |
| 3154 | |
| 3155 | ptlbls = ' names :' |
| 3156 | |
| 3157 | ptstr = ' values:' |
| 3158 | |
| 3159 | sigstr = ' sig :' |
| 3160 | |
| 3161 | refine = False |
| 3162 | |
| 3163 | if 'Bragg' in Sample['Type']: |
| 3164 | |
| 3165 | for i,item in enumerate(['Scale','Shift','Transparency']): |
| 3166 | |
| 3167 | ptlbls += '%14s'%(item) |
| 3168 | |
| 3169 | ptstr += '%14.4f'%(Sample[item][0]) |
| 3170 | |
| 3171 | if sampleSig[i]: |
| 3172 | |
| 3173 | refine = True |
| 3174 | |
| 3175 | sigstr += '%14.4f'%(sampleSig[i]) |
| 3176 | |
| 3177 | else: |
| 3178 | |
| 3179 | sigstr += 14*' ' |
| 3180 | |
| 3181 | |
| 3182 | |
| 3183 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
| 3184 | |
| 3185 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
| 3186 | |
| 3187 | ptlbls += '%14s'%(item) |
| 3188 | |
| 3189 | ptstr += '%14.4f'%(Sample[item][0]) |
| 3190 | |
| 3191 | if sampleSig[i]: |
| 3192 | |
| 3193 | refine = True |
| 3194 | |
| 3195 | sigstr += '%14.4f'%(sampleSig[i]) |
| 3196 | |
| 3197 | else: |
| 3198 | |
| 3199 | sigstr += 14*' ' |
| 3200 | |
| 3201 | |
| 3202 | |
| 3203 | if refine: |
| 3204 | |
| 3205 | print '\n Sample Parameters:' |
| 3206 | |
| 3207 | print ptlbls |
| 3208 | |
| 3209 | print ptstr |
| 3210 | |
| 3211 | print sigstr |
| 3212 | |
| 3213 | |
| 3214 | |
| 3215 | histoList = Histograms.keys() |
| 3216 | |
| 3217 | histoList.sort() |
| 3218 | |
| 3219 | for histogram in histoList: |
| 3220 | |
| 3221 | if 'PWDR' in histogram: |
| 3222 | |
| 3223 | Histogram = Histograms[histogram] |
| 3224 | |
| 3225 | hId = Histogram['hId'] |
| 3226 | |
| 3227 | pfx = ':'+str(hId)+':' |
| 3228 | |
| 3229 | Background = Histogram['Background'] |
| 3230 | |
| 3231 | backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict) |
| 3232 | |
| 3233 | |
| 3234 | |
| 3235 | Inst = Histogram['Instrument Parameters'] |
| 3236 | |
| 3237 | instSig = SetInstParms(pfx,Inst,parmDict,sigDict) |
| 3238 | |
| 3239 | |
| 3240 | |
| 3241 | Sample = Histogram['Sample Parameters'] |
| 3242 | |
| 3243 | sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict) |
| 3244 | |
| 3245 | |
| 3246 | |
| 3247 | print '\n Histogram: ',histogram,' histogram Id: ',hId |
| 3248 | |
| 3249 | print 135*'-' |
| 3250 | |
| 3251 | print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs']) |
| 3252 | |
| 3253 | if Print: |
| 3254 | |
| 3255 | print ' Instrument type: ',Sample['Type'] |
| 3256 | |
| 3257 | PrintSampleParmsSig(Sample,sampSig) |
| 3258 | |
| 3259 | PrintInstParmsSig(Inst,instSig) |
| 3260 | |
| 3261 | PrintBackgroundSig(Background,backSig) |
| 3262 | |
| 3263 | |
| 3264 | |
| 3265 | def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict): |
| 3266 | |
| 3267 | Natoms = calcControls['Natoms'][pfx] |
| 3268 | |
| 3269 | Tdata = Natoms*[' ',] |
| 3270 | |
| 3271 | Mdata = np.zeros(Natoms) |
| 3272 | |
| 3273 | IAdata = Natoms*[' ',] |
| 3274 | |
| 3275 | Fdata = np.zeros(Natoms) |
| 3276 | |
| 3277 | FFdata = [] |
| 3278 | |
| 3279 | BLdata = [] |
| 3280 | |
| 3281 | Xdata = np.zeros((3,Natoms)) |
| 3282 | |
| 3283 | dXdata = np.zeros((3,Natoms)) |
| 3284 | |
| 3285 | Uisodata = np.zeros(Natoms) |
| 3286 | |
| 3287 | Uijdata = np.zeros((6,Natoms)) |
| 3288 | |
| 3289 | keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata, |
| 3290 | |
| 3291 | 'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2], |
| 3292 | |
| 3293 | 'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata, |
| 3294 | |
| 3295 | 'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2], |
| 3296 | |
| 3297 | 'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]} |
| 3298 | |
| 3299 | for iatm in range(Natoms): |
| 3300 | |
| 3301 | for key in keys: |
| 3302 | |
| 3303 | parm = pfx+key+str(iatm) |
| 3304 | |
| 3305 | if parm in parmDict: |
| 3306 | |
| 3307 | keys[key][iatm] = parmDict[parm] |
| 3308 | |
| 3309 | FFdata.append(FFtables[Tdata[iatm]]) |
| 3310 | |
| 3311 | BLdata.append(BLtables[Tdata[iatm]][1]) |
| 3312 | |
| 3313 | return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata |
| 3314 | |
| 3315 | |
| 3316 | |
| 3317 | def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict): |
| 3318 | |
| 3319 | ''' Compute structure factors for all h,k,l for phase |
| 3320 | |
271 | | |
272 | | GPXback = GPXBackup(GPXfile,makeBack) |
273 | | print '\n',135*'-' |
274 | | print 'Read from file:',GPXback |
275 | | print 'Save to file :',GPXfile |
276 | | infile = open(GPXback,'rb') |
277 | | outfile = open(GPXfile,'wb') |
| 3664 | |
| 3665 | newAtomDict = {} |
| 3666 | |
| 3667 | for item in parmDict: |
| 3668 | |
| 3669 | if 'dA' in item: |
| 3670 | |
| 3671 | parm = ''.join(item.split('d')) |
| 3672 | |
| 3673 | parmDict[parm] += parmDict[item] |
| 3674 | |
| 3675 | newAtomDict[item] = [parm,parmDict[parm]] |
| 3676 | |
| 3677 | return newAtomDict |
| 3678 | |
| 3679 | |
| 3680 | |
| 3681 | def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict): |
| 3682 | |
| 3683 | IFCoup = 'Bragg' in calcControls[hfx+'instType'] |
| 3684 | |
| 3685 | odfCor = 1.0 |
| 3686 | |
| 3687 | H = refl[:3] |
| 3688 | |
| 3689 | cell = G2lat.Gmat2cell(g) |
| 3690 | |
| 3691 | Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']] |
| 3692 | |
| 3693 | Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] |
| 3694 | |
| 3695 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 3696 | |
| 3697 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs. |
| 3698 | |
| 3699 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) |
| 3700 | |
| 3701 | for item in SHnames: |
| 3702 | |
| 3703 | L,M,N = eval(item.strip('C')) |
| 3704 | |
| 3705 | Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta) |
| 3706 | |
| 3707 | Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam) |
| 3708 | |
| 3709 | Lnorm = G2lat.Lnorm(L) |
| 3710 | |
| 3711 | odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl |
| 3712 | |
| 3713 | return odfCor |
| 3714 | |
| 3715 | |
| 3716 | |
| 3717 | def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict): |
| 3718 | |
| 3719 | FORPI = 12.5663706143592 |
| 3720 | |
| 3721 | IFCoup = 'Bragg' in calcControls[hfx+'instType'] |
| 3722 | |
| 3723 | odfCor = 1.0 |
| 3724 | |
| 3725 | dFdODF = {} |
| 3726 | |
| 3727 | dFdSA = [0,0,0] |
| 3728 | |
| 3729 | H = refl[:3] |
| 3730 | |
| 3731 | cell = G2lat.Gmat2cell(g) |
| 3732 | |
| 3733 | Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']] |
| 3734 | |
| 3735 | Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] |
| 3736 | |
| 3737 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 3738 | |
| 3739 | psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) |
| 3740 | |
| 3741 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) |
| 3742 | |
| 3743 | for item in SHnames: |
| 3744 | |
| 3745 | L,M,N = eval(item.strip('C')) |
| 3746 | |
| 3747 | Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta) |
| 3748 | |
| 3749 | Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam) |
| 3750 | |
| 3751 | Lnorm = G2lat.Lnorm(L) |
| 3752 | |
| 3753 | odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl |
| 3754 | |
| 3755 | dFdODF[pfx+item] = Lnorm*Kcl*Ksl |
| 3756 | |
| 3757 | for i in range(3): |
| 3758 | |
| 3759 | dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i]) |
| 3760 | |
| 3761 | return odfCor,dFdODF,dFdSA |
| 3762 | |
| 3763 | |
| 3764 | |
| 3765 | def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict): |
| 3766 | |
| 3767 | odfCor = 1.0 |
| 3768 | |
| 3769 | H = refl[:3] |
| 3770 | |
| 3771 | cell = G2lat.Gmat2cell(g) |
| 3772 | |
| 3773 | Sangl = [0.,0.,0.] |
| 3774 | |
| 3775 | if 'Bragg' in calcControls[hfx+'instType']: |
| 3776 | |
| 3777 | Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']] |
| 3778 | |
| 3779 | IFCoup = True |
| 3780 | |
| 3781 | else: |
| 3782 | |
| 3783 | Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']] |
| 3784 | |
| 3785 | IFCoup = False |
| 3786 | |
| 3787 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 3788 | |
| 3789 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. |
| 3790 | |
| 3791 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) |
| 3792 | |
| 3793 | for item in SHnames: |
| 3794 | |
| 3795 | L,N = eval(item.strip('C')) |
| 3796 | |
| 3797 | Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) |
| 3798 | |
| 3799 | odfCor += parmDict[phfx+item]*Lnorm*Kcsl |
| 3800 | |
| 3801 | return odfCor |
| 3802 | |
| 3803 | |
| 3804 | |
| 3805 | def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict): |
| 3806 | |
| 3807 | FORPI = 12.5663706143592 |
| 3808 | |
| 3809 | odfCor = 1.0 |
| 3810 | |
| 3811 | dFdODF = {} |
| 3812 | |
| 3813 | H = refl[:3] |
| 3814 | |
| 3815 | cell = G2lat.Gmat2cell(g) |
| 3816 | |
| 3817 | Sangl = [0.,0.,0.] |
| 3818 | |
| 3819 | if 'Bragg' in calcControls[hfx+'instType']: |
| 3820 | |
| 3821 | Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']] |
| 3822 | |
| 3823 | IFCoup = True |
| 3824 | |
| 3825 | else: |
| 3826 | |
| 3827 | Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']] |
| 3828 | |
| 3829 | IFCoup = False |
| 3830 | |
| 3831 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
| 3832 | |
| 3833 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. |
| 3834 | |
| 3835 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) |
| 3836 | |
| 3837 | for item in SHnames: |
| 3838 | |
| 3839 | L,N = eval(item.strip('C')) |
| 3840 | |
| 3841 | Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) |
| 3842 | |
| 3843 | odfCor += parmDict[phfx+item]*Lnorm*Kcsl |
| 3844 | |
| 3845 | dFdODF[phfx+item] = Kcsl*Lnorm |
| 3846 | |
| 3847 | return odfCor,dFdODF |
| 3848 | |
| 3849 | |
| 3850 | |
| 3851 | def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): |
| 3852 | |
| 3853 | if calcControls[phfx+'poType'] == 'MD': |
| 3854 | |
| 3855 | MD = parmDict[phfx+'MD'] |
| 3856 | |
| 3857 | MDAxis = calcControls[phfx+'MDAxis'] |
| 3858 | |
| 3859 | sumMD = 0 |
| 3860 | |
| 3861 | for H in refl[11]: |
| 3862 | |
| 3863 | cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) |
| 3864 | |
| 3865 | A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) |
| 3866 | |
| 3867 | sumMD += A**3 |
| 3868 | |
| 3869 | POcorr = sumMD/len(refl[11]) |
| 3870 | |
| 3871 | else: #spherical harmonics |
| 3872 | |
| 3873 | POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict) |
| 3874 | |
| 3875 | return POcorr |
| 3876 | |
| 3877 | |
| 3878 | |
| 3879 | def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): |
| 3880 | |
| 3881 | POderv = {} |
| 3882 | |
| 3883 | if calcControls[phfx+'poType'] == 'MD': |
| 3884 | |
| 3885 | MD = parmDict[phfx+'MD'] |
| 3886 | |
| 3887 | MDAxis = calcControls[phfx+'MDAxis'] |
| 3888 | |
| 3889 | sumMD = 0 |
| 3890 | |
| 3891 | sumdMD = 0 |
| 3892 | |
| 3893 | for H in refl[11]: |
| 3894 | |
| 3895 | cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) |
| 3896 | |
| 3897 | A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) |
| 3898 | |
| 3899 | sumMD += A**3 |
| 3900 | |
| 3901 | sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2) |
| 3902 | |
| 3903 | POcorr = sumMD/len(refl[11]) |
| 3904 | |
| 3905 | POderv[phfx+'MD'] = sumdMD/len(refl[11]) |
| 3906 | |
| 3907 | else: #spherical harmonics |
| 3908 | |
| 3909 | POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict) |
| 3910 | |
| 3911 | return POcorr,POderv |
| 3912 | |
| 3913 | |
| 3914 | |
| 3915 | def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): |
| 3916 | |
| 3917 | Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity |
| 3918 | |
| 3919 | if 'X' in parmDict[hfx+'Type']: |
| 3920 | |
| 3921 | Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0] |
| 3922 | |
| 3923 | Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) |
| 3924 | |
| 3925 | if pfx+'SHorder' in parmDict: |
| 3926 | |
| 3927 | Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict) |
| 3928 | |
| 3929 | refl[13] = Icorr |
| 3930 | |
| 3931 | |
| 3932 | |
| 3933 | def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): |
| 3934 | |
| 3935 | dIdsh = 1./parmDict[hfx+'Scale'] |
| 3936 | |
| 3937 | dIdsp = 1./parmDict[phfx+'Scale'] |
| 3938 | |
| 3939 | if 'X' in parmDict[hfx+'Type']: |
| 3940 | |
| 3941 | pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth']) |
| 3942 | |
| 3943 | dIdPola /= pola |
| 3944 | |
| 3945 | else: #'N' |
| 3946 | |
| 3947 | dIdPola = 0.0 |
| 3948 | |
| 3949 | POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) |
| 3950 | |
| 3951 | for iPO in dIdPO: |
| 3952 | |
| 3953 | dIdPO[iPO] /= POcorr |
| 3954 | |
| 3955 | dFdODF = {} |
| 3956 | |
| 3957 | dFdSA = [0,0,0] |
| 3958 | |
| 3959 | if pfx+'SHorder' in parmDict: |
| 3960 | |
| 3961 | odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict) |
| 3962 | |
| 3963 | for iSH in dFdODF: |
| 3964 | |
| 3965 | dFdODF[iSH] /= odfCor |
| 3966 | |
| 3967 | for i in range(3): |
| 3968 | |
| 3969 | dFdSA[i] /= odfCor |
| 3970 | |
| 3971 | return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA |
| 3972 | |
| 3973 | |
| 3974 | |
| 3975 | def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict): |
| 3976 | |
| 3977 | costh = cosd(refl[5]/2.) |
| 3978 | |
| 3979 | #crystallite size |
| 3980 | |
| 3981 | if calcControls[phfx+'SizeType'] == 'isotropic': |
| 3982 | |
| 3983 | gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh) |
| 3984 | |
| 3985 | elif calcControls[phfx+'SizeType'] == 'uniaxial': |
| 3986 | |
| 3987 | H = np.array(refl[:3]) |
| 3988 | |
| 3989 | P = np.array(calcControls[phfx+'SizeAxis']) |
| 3990 | |
| 3991 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 3992 | |
| 3993 | gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh) |
| 3994 | |
| 3995 | gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2) |
| 3996 | |
| 3997 | else: #ellipsoidal crystallites |
| 3998 | |
| 3999 | Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] |
| 4000 | |
| 4001 | H = np.array(refl[:3]) |
| 4002 | |
| 4003 | lenR = G2pwd.ellipseSize(H,Sij,GB) |
| 4004 | |
| 4005 | gam = 1.8*wave/(np.pi*costh*lenR) |
| 4006 | |
| 4007 | #microstrain |
| 4008 | |
| 4009 | if calcControls[phfx+'MustrainType'] == 'isotropic': |
| 4010 | |
| 4011 | gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi |
| 4012 | |
| 4013 | elif calcControls[phfx+'MustrainType'] == 'uniaxial': |
| 4014 | |
| 4015 | H = np.array(refl[:3]) |
| 4016 | |
| 4017 | P = np.array(calcControls[phfx+'MustrainAxis']) |
| 4018 | |
| 4019 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 4020 | |
| 4021 | Si = parmDict[phfx+'Mustrain:i'] |
| 4022 | |
| 4023 | Sa = parmDict[phfx+'Mustrain:a'] |
| 4024 | |
| 4025 | gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) |
| 4026 | |
| 4027 | else: #generalized - P.W. Stephens model |
| 4028 | |
| 4029 | pwrs = calcControls[phfx+'MuPwrs'] |
| 4030 | |
| 4031 | sum = 0 |
| 4032 | |
| 4033 | for i,pwr in enumerate(pwrs): |
| 4034 | |
| 4035 | sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] |
| 4036 | |
| 4037 | gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum |
| 4038 | |
| 4039 | return gam |
| 4040 | |
| 4041 | |
| 4042 | |
| 4043 | def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict): |
| 4044 | |
| 4045 | gamDict = {} |
| 4046 | |
| 4047 | costh = cosd(refl[5]/2.) |
| 4048 | |
| 4049 | tanth = tand(refl[5]/2.) |
| 4050 | |
| 4051 | #crystallite size derivatives |
| 4052 | |
| 4053 | if calcControls[phfx+'SizeType'] == 'isotropic': |
| 4054 | |
| 4055 | gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh) |
| 4056 | |
| 4057 | elif calcControls[phfx+'SizeType'] == 'uniaxial': |
| 4058 | |
| 4059 | H = np.array(refl[:3]) |
| 4060 | |
| 4061 | P = np.array(calcControls[phfx+'SizeAxis']) |
| 4062 | |
| 4063 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 4064 | |
| 4065 | Si = parmDict[phfx+'Size:i'] |
| 4066 | |
| 4067 | Sa = parmDict[phfx+'Size:a'] |
| 4068 | |
| 4069 | gami = (1.8*wave/np.pi)/(Si*Sa) |
| 4070 | |
| 4071 | sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2) |
| 4072 | |
| 4073 | gam = gami*sqtrm/costh |
| 4074 | |
| 4075 | gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si |
| 4076 | |
| 4077 | gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa |
| 4078 | |
| 4079 | else: #ellipsoidal crystallites |
| 4080 | |
| 4081 | const = 1.8*wave/(np.pi*costh) |
| 4082 | |
| 4083 | Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] |
| 4084 | |
| 4085 | H = np.array(refl[:3]) |
| 4086 | |
| 4087 | R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) |
| 4088 | |
| 4089 | for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): |
| 4090 | |
| 4091 | gamDict[item] = -(const/R**2)*dRdS[i] |
| 4092 | |
| 4093 | #microstrain derivatives |
| 4094 | |
| 4095 | if calcControls[phfx+'MustrainType'] == 'isotropic': |
| 4096 | |
| 4097 | gamDict[phfx+'Mustrain:i'] = 0.018*tanth/np.pi |
| 4098 | |
| 4099 | elif calcControls[phfx+'MustrainType'] == 'uniaxial': |
| 4100 | |
| 4101 | H = np.array(refl[:3]) |
| 4102 | |
| 4103 | P = np.array(calcControls[phfx+'MustrainAxis']) |
| 4104 | |
| 4105 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
| 4106 | |
| 4107 | Si = parmDict[phfx+'Mustrain:i'] |
| 4108 | |
| 4109 | Sa = parmDict[phfx+'Mustrain:a'] |
| 4110 | |
| 4111 | gami = 0.018*Si*Sa*tanth/np.pi |
| 4112 | |
| 4113 | sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) |
| 4114 | |
| 4115 | gam = gami/sqtrm |
| 4116 | |
| 4117 | gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3 |
| 4118 | |
| 4119 | gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3 |
| 4120 | |
| 4121 | else: #generalized - P.W. Stephens model |
| 4122 | |
| 4123 | pwrs = calcControls[phfx+'MuPwrs'] |
| 4124 | |
| 4125 | const = 0.018*refl[4]**2*tanth |
| 4126 | |
| 4127 | for i,pwr in enumerate(pwrs): |
| 4128 | |
| 4129 | gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] |
| 4130 | |
| 4131 | return gamDict |
| 4132 | |
| 4133 | |
| 4134 | |
| 4135 | def GetReflPos(refl,wave,G,hfx,calcControls,parmDict): |
| 4136 | |
| 4137 | h,k,l = refl[:3] |
| 4138 | |
| 4139 | dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G) |
| 4140 | |
| 4141 | d = np.sqrt(dsq) |
| 4142 | |
| 4143 | |
| 4144 | |
| 4145 | refl[4] = d |
| 4146 | |
| 4147 | pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero'] |
| 4148 | |
| 4149 | const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns |
| 4150 | |
| 4151 | if 'Bragg' in calcControls[hfx+'instType']: |
| 4152 | |
| 4153 | pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \ |
| 4154 | |
| 4155 | parmDict[hfx+'Transparency']*sind(pos)*100.0) #trans(=1/mueff) in cm |
| 4156 | |
| 4157 | else: #Debye-Scherrer - simple but maybe not right |
| 4158 | |
| 4159 | pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos)) |
| 4160 | |
| 4161 | return pos |
| 4162 | |
| 4163 | |
| 4164 | |
| 4165 | def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict): |
| 4166 | |
| 4167 | dpr = 180./np.pi |
| 4168 | |
| 4169 | h,k,l = refl[:3] |
| 4170 | |
| 4171 | dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A) |
| 4172 | |
| 4173 | dst = np.sqrt(dstsq) |
| 4174 | |
| 4175 | pos = refl[5]-parmDict[hfx+'Zero'] |
| 4176 | |
| 4177 | const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0) |
| 4178 | |
| 4179 | dpdw = const*dst |
| 4180 | |
| 4181 | dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) |
| 4182 | |
| 4183 | dpdA *= const*wave/(2.0*dst) |
| 4184 | |
| 4185 | dpdZ = 1.0 |
| 4186 | |
| 4187 | const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns |
| 4188 | |
| 4189 | if 'Bragg' in calcControls[hfx+'instType']: |
| 4190 | |
| 4191 | dpdSh = -4.*const*cosd(pos/2.0) |
| 4192 | |
| 4193 | dpdTr = -const*sind(pos)*100.0 |
| 4194 | |
| 4195 | return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0. |
| 4196 | |
| 4197 | else: #Debye-Scherrer - simple but maybe not right |
| 4198 | |
| 4199 | dpdXd = -const*cosd(pos) |
| 4200 | |
| 4201 | dpdYd = -const*sind(pos) |
| 4202 | |
| 4203 | return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd |
| 4204 | |
| 4205 | |
| 4206 | |
| 4207 | def GetHStrainShift(refl,SGData,phfx,parmDict): |
| 4208 | |
| 4209 | laue = SGData['SGLaue'] |
| 4210 | |
| 4211 | uniq = SGData['SGUniq'] |
| 4212 | |
| 4213 | h,k,l = refl[:3] |
| 4214 | |
| 4215 | if laue in ['m3','m3m']: |
| 4216 | |
| 4217 | Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \ |
| 4218 | |
| 4219 | refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2 |
| 4220 | |
| 4221 | elif laue in ['6/m','6/mmm','3m1','31m','3']: |
| 4222 | |
| 4223 | Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2 |
| 4224 | |
| 4225 | elif laue in ['3R','3mR']: |
| 4226 | |
| 4227 | Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l) |
| 4228 | |
| 4229 | elif laue in ['4/m','4/mmm']: |
| 4230 | |
| 4231 | Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2 |
| 4232 | |
| 4233 | elif laue in ['mmm']: |
| 4234 | |
| 4235 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2 |
| 4236 | |
| 4237 | elif laue in ['2/m']: |
| 4238 | |
| 4239 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2 |
| 4240 | |
| 4241 | if uniq == 'a': |
| 4242 | |
| 4243 | Dij += parmDict[phfx+'D23']*k*l |
| 4244 | |
| 4245 | elif uniq == 'b': |
| 4246 | |
| 4247 | Dij += parmDict[phfx+'D13']*h*l |
| 4248 | |
| 4249 | elif uniq == 'c': |
| 4250 | |
| 4251 | Dij += parmDict[phfx+'D12']*h*k |
| 4252 | |
| 4253 | else: |
| 4254 | |
| 4255 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \ |
| 4256 | |
| 4257 | parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l |
| 4258 | |
| 4259 | return Dij*refl[4]**2*tand(refl[5]/2.0) |
| 4260 | |
| 4261 | |
| 4262 | |
| 4263 | def GetHStrainShiftDerv(refl,SGData,phfx): |
| 4264 | |
| 4265 | laue = SGData['SGLaue'] |
| 4266 | |
| 4267 | uniq = SGData['SGUniq'] |
| 4268 | |
| 4269 | h,k,l = refl[:3] |
| 4270 | |
| 4271 | if laue in ['m3','m3m']: |
| 4272 | |
| 4273 | dDijDict = {phfx+'D11':h**2+k**2+l**2, |
| 4274 | |
| 4275 | phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2} |
| 4276 | |
| 4277 | elif laue in ['6/m','6/mmm','3m1','31m','3']: |
| 4278 | |
| 4279 | dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2} |
| 4280 | |
| 4281 | elif laue in ['3R','3mR']: |
| 4282 | |
| 4283 | dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l} |
| 4284 | |
| 4285 | elif laue in ['4/m','4/mmm']: |
| 4286 | |
| 4287 | dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2} |
| 4288 | |
| 4289 | elif laue in ['mmm']: |
| 4290 | |
| 4291 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2} |
| 4292 | |
| 4293 | elif laue in ['2/m']: |
| 4294 | |
| 4295 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2} |
| 4296 | |
| 4297 | if uniq == 'a': |
| 4298 | |
| 4299 | dDijDict[phfx+'D23'] = k*l |
| 4300 | |
| 4301 | elif uniq == 'b': |
| 4302 | |
| 4303 | dDijDict[phfx+'D13'] = h*l |
| 4304 | |
| 4305 | elif uniq == 'c': |
| 4306 | |
| 4307 | dDijDict[phfx+'D12'] = h*k |
| 4308 | |
| 4309 | names.append() |
| 4310 | |
| 4311 | else: |
| 4312 | |
| 4313 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2, |
| 4314 | |
| 4315 | phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l} |
| 4316 | |
| 4317 | for item in dDijDict: |
| 4318 | |
| 4319 | dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0) |
| 4320 | |
| 4321 | return dDijDict |
| 4322 | |
| 4323 | |
| 4324 | |
| 4325 | def GetFprime(controlDict,Histograms): |
| 4326 | |
| 4327 | FFtables = controlDict['FFtables'] |
| 4328 | |
| 4329 | if not FFtables: |
| 4330 | |
| 4331 | return |
| 4332 | |
| 4333 | histoList = Histograms.keys() |
| 4334 | |
| 4335 | histoList.sort() |
| 4336 | |
| 4337 | for histogram in histoList: |
| 4338 | |
| 4339 | if 'PWDR' in histogram[:4]: |
| 4340 | |
| 4341 | Histogram = Histograms[histogram] |
| 4342 | |
| 4343 | hId = Histogram['hId'] |
| 4344 | |
| 4345 | hfx = ':%d:'%(hId) |
| 4346 | |
| 4347 | keV = controlDict[hfx+'keV'] |
| 4348 | |
| 4349 | for El in FFtables: |
| 4350 | |
| 4351 | Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0]) |
| 4352 | |
| 4353 | FP,FPP,Mu = G2el.FPcalc(Orbs, keV) |
| 4354 | |
| 4355 | FFtables[El][hfx+'FP'] = FP |
| 4356 | |
| 4357 | FFtables[El][hfx+'FPP'] = FPP |
| 4358 | |
| 4359 | |
| 4360 | |
| 4361 | def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): |
| 4362 | |
| 4363 | |
| 4364 | |
| 4365 | def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): |
| 4366 | |
| 4367 | U = parmDict[hfx+'U'] |
| 4368 | |
| 4369 | V = parmDict[hfx+'V'] |
| 4370 | |
| 4371 | W = parmDict[hfx+'W'] |
| 4372 | |
| 4373 | X = parmDict[hfx+'X'] |
| 4374 | |
| 4375 | Y = parmDict[hfx+'Y'] |
| 4376 | |
| 4377 | tanPos = tand(refl[5]/2.0) |
| 4378 | |
| 4379 | sig = U*tanPos**2+V*tanPos+W #save peak sigma |
| 4380 | |
| 4381 | sig = max(0.001,sig) |
| 4382 | |
| 4383 | gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma |
| 4384 | |
| 4385 | gam = max(0.001,gam) |
| 4386 | |
| 4387 | return sig,gam |
| 4388 | |
| 4389 | |
| 4390 | |
| 4391 | hId = Histogram['hId'] |
| 4392 | |
| 4393 | hfx = ':%d:'%(hId) |
| 4394 | |
| 4395 | bakType = calcControls[hfx+'bakType'] |
| 4396 | |
| 4397 | yb = G2pwd.getBackground(hfx,parmDict,bakType,x) |
| 4398 | |
| 4399 | yc = np.zeros_like(yb) |
| 4400 | |
| 4401 | |
| 4402 | |
| 4403 | if 'C' in calcControls[hfx+'histType']: |
| 4404 | |
| 4405 | shl = max(parmDict[hfx+'SH/L'],0.002) |
| 4406 | |
| 4407 | Ka2 = False |
| 4408 | |
| 4409 | if hfx+'Lam1' in parmDict.keys(): |
| 4410 | |
| 4411 | wave = parmDict[hfx+'Lam1'] |
| 4412 | |
| 4413 | Ka2 = True |
| 4414 | |
| 4415 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
| 4416 | |
| 4417 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
| 4418 | |
| 4419 | else: |
| 4420 | |
| 4421 | wave = parmDict[hfx+'Lam'] |
| 4422 | |
| 4423 | else: |
| 4424 | |
| 4425 | print 'TOF Undefined at present' |
| 4426 | |
| 4427 | raise ValueError |
| 4428 | |
| 4429 | for phase in Histogram['Reflection Lists']: |
| 4430 | |
| 4431 | refList = Histogram['Reflection Lists'][phase] |
| 4432 | |
| 4433 | Phase = Phases[phase] |
| 4434 | |
| 4435 | pId = Phase['pId'] |
| 4436 | |
| 4437 | pfx = '%d::'%(pId) |
| 4438 | |
| 4439 | phfx = '%d:%d:'%(pId,hId) |
| 4440 | |
| 4441 | hfx = ':%d:'%(hId) |
| 4442 | |
| 4443 | SGData = Phase['General']['SGData'] |
| 4444 | |
| 4445 | A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] |
| 4446 | |
| 4447 | G,g = G2lat.A2Gmat(A) #recip & real metric tensors |
| 4448 | |
| 4449 | GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies |
| 4450 | |
| 4451 | Vst = np.sqrt(nl.det(G)) #V* |
| 4452 | |
| 4453 | if 'Pawley' not in Phase['General']['Type']: |
| 4454 | |
| 4455 | refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict) |
| 4456 | |
| 4457 | for refl in refList: |
| 4458 | |
| 4459 | if 'C' in calcControls[hfx+'histType']: |
| 4460 | |
| 4461 | h,k,l = refl[:3] |
| 4462 | |
| 4463 | refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position |
| 4464 | |
| 4465 | Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.)) #Lorentz correction |
| 4466 | |
| 4467 | refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift |
| 4468 | |
| 4469 | refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam |
| 4470 | |
| 4471 | GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13] |
| 4472 | |
| 4473 | refl[13] *= Vst*Lorenz |
| 4474 | |
| 4475 | if 'Pawley' in Phase['General']['Type']: |
| 4476 | |
| 4477 | try: |
| 4478 | |
| 4479 | refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) |
| 4480 | |
| 4481 | except KeyError: |
| 4482 | |
| 4483 | # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) |
| 4484 | |
| 4485 | continue |
| 4486 | |
| 4487 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
| 4488 | |
| 4489 | iBeg = np.searchsorted(x,refl[5]-fmin) |
| 4490 | |
| 4491 | iFin = np.searchsorted(x,refl[5]+fmax) |
| 4492 | |
| 4493 | if not iBeg+iFin: #peak below low limit - skip peak |
| 4494 | |
| 4495 | continue |
| 4496 | |
| 4497 | elif not iBeg-iFin: #peak above high limit - done |
| 4498 | |
| 4499 | break |
| 4500 | |
| 4501 | yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here |
| 4502 | |
| 4503 | if Ka2: |
| 4504 | |
| 4505 | pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) |
| 4506 | |
| 4507 | Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) |
| 4508 | |
| 4509 | iBeg = np.searchsorted(x,pos2-fmin) |
| 4510 | |
| 4511 | iFin = np.searchsorted(x,pos2+fmax) |
| 4512 | |
| 4513 | if not iBeg+iFin: #peak below low limit - skip peak |
| 4514 | |
| 4515 | continue |
| 4516 | |
| 4517 | elif not iBeg-iFin: #peak above high limit - done |
| 4518 | |
| 4519 | return yc,yb |
| 4520 | |
| 4521 | yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here |
| 4522 | |
| 4523 | elif 'T' in calcControls[hfx+'histType']: |
| 4524 | |
| 4525 | print 'TOF Undefined at present' |
| 4526 | |
| 4527 | raise Exception #no TOF yet |
| 4528 | |
| 4529 | return yc,yb |
| 4530 | |
| 4531 | |
| 4532 | |
| 4533 | def GetFobsSq(Histograms,Phases,parmDict,calcControls): |
| 4534 | |
| 4535 | histoList = Histograms.keys() |
| 4536 | |
| 4537 | histoList.sort() |
| 4538 | |
| 4539 | for histogram in histoList: |
| 4540 | |
| 4541 | if 'PWDR' in histogram[:4]: |
| 4542 | |
| 4543 | Histogram = Histograms[histogram] |
| 4544 | |
| 4545 | hId = Histogram['hId'] |
| 4546 | |
| 4547 | hfx = ':%d:'%(hId) |
| 4548 | |
| 4549 | Limits = calcControls[hfx+'Limits'] |
| 4550 | |
| 4551 | shl = max(parmDict[hfx+'SH/L'],0.002) |
| 4552 | |
| 4553 | Ka2 = False |
| 4554 | |
| 4555 | kRatio = 0.0 |
| 4556 | |
| 4557 | if hfx+'Lam1' in parmDict.keys(): |
| 4558 | |
| 4559 | Ka2 = True |
| 4560 | |
| 4561 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
| 4562 | |
| 4563 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
| 4564 | |
| 4565 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 4566 | |
| 4567 | ymb = np.array(y-yb) |
| 4568 | |
| 4569 | ycmb = np.array(yc-yb) |
| 4570 | |
| 4571 | ratio = np.where(ycmb!=0.,ymb/ycmb,0.0) |
| 4572 | |
| 4573 | xB = np.searchsorted(x,Limits[0]) |
| 4574 | |
| 4575 | xF = np.searchsorted(x,Limits[1]) |
| 4576 | |
| 4577 | refLists = Histogram['Reflection Lists'] |
| 4578 | |
| 4579 | for phase in refLists: |
| 4580 | |
| 4581 | Phase = Phases[phase] |
| 4582 | |
| 4583 | pId = Phase['pId'] |
| 4584 | |
| 4585 | phfx = '%d:%d:'%(pId,hId) |
| 4586 | |
| 4587 | refList = refLists[phase] |
| 4588 | |
| 4589 | sumFo = 0.0 |
| 4590 | |
| 4591 | sumdF = 0.0 |
| 4592 | |
| 4593 | sumFosq = 0.0 |
| 4594 | |
| 4595 | sumdFsq = 0.0 |
| 4596 | |
| 4597 | for refl in refList: |
| 4598 | |
| 4599 | if 'C' in calcControls[hfx+'histType']: |
| 4600 | |
| 4601 | yp = np.zeros_like(yb) |
| 4602 | |
| 4603 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
| 4604 | |
| 4605 | iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin) |
| 4606 | |
| 4607 | iFin = np.searchsorted(x[xB:xF],refl[5]+fmax) |
| 4608 | |
| 4609 | iFin2 = iFin |
| 4610 | |
| 4611 | yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here |
| 4612 | |
| 4613 | if Ka2: |
| 4614 | |
| 4615 | pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) |
| 4616 | |
| 4617 | Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) |
| 4618 | |
| 4619 | iBeg2 = np.searchsorted(x,pos2-fmin) |
| 4620 | |
| 4621 | iFin2 = np.searchsorted(x,pos2+fmax) |
| 4622 | |
| 4623 | yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here |
| 4624 | |
| 4625 | refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0)) |
| 4626 | |
| 4627 | elif 'T' in calcControls[hfx+'histType']: |
| 4628 | |
| 4629 | print 'TOF Undefined at present' |
| 4630 | |
| 4631 | raise Exception #no TOF yet |
| 4632 | |
| 4633 | Fo = np.sqrt(np.abs(refl[8])) |
| 4634 | |
| 4635 | Fc = np.sqrt(np.abs(refl[9])) |
| 4636 | |
| 4637 | sumFo += Fo |
| 4638 | |
| 4639 | sumFosq += refl[8]**2 |
| 4640 | |
| 4641 | sumdF += np.abs(Fo-Fc) |
| 4642 | |
| 4643 | sumdFsq += (refl[8]-refl[9])**2 |
| 4644 | |
| 4645 | Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.) |
| 4646 | |
| 4647 | Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.) |
| 4648 | |
| 4649 | Histogram[phfx+'Nref'] = len(refList) |
| 4650 | |
| 4651 | |
| 4652 | |
| 4653 | def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): |
| 4654 | |
| 4655 | |
| 4656 | |
| 4657 | def cellVaryDerv(pfx,SGData,dpdA): |
| 4658 | |
| 4659 | if SGData['SGLaue'] in ['-1',]: |
| 4660 | |
| 4661 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]], |
| 4662 | |
| 4663 | [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]] |
| 4664 | |
| 4665 | elif SGData['SGLaue'] in ['2/m',]: |
| 4666 | |
| 4667 | if SGData['SGUniq'] == 'a': |
| 4668 | |
| 4669 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]] |
| 4670 | |
| 4671 | elif SGData['SGUniq'] == 'b': |
| 4672 | |
| 4673 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]] |
| 4674 | |
| 4675 | else: |
| 4676 | |
| 4677 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]] |
| 4678 | |
| 4679 | elif SGData['SGLaue'] in ['mmm',]: |
| 4680 | |
| 4681 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]] |
| 4682 | |
| 4683 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
| 4684 | |
| 4685 | # return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]] |
| 4686 | |
| 4687 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
| 4688 | |
| 4689 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
| 4690 | |
| 4691 | # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]] |
| 4692 | |
| 4693 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
| 4694 | |
| 4695 | elif SGData['SGLaue'] in ['3R', '3mR']: |
| 4696 | |
| 4697 | return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]] |
| 4698 | |
| 4699 | elif SGData['SGLaue'] in ['m3m','m3']: |
| 4700 | |
| 4701 | # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]] |
| 4702 | |
| 4703 | return [[pfx+'A0',dpdA[0]]] |
| 4704 | |
| 4705 | # create a list of dependent variables and set up a dictionary to hold their derivatives |
| 4706 | |
| 4707 | dependentVars = G2mv.GetDependentVars() |
| 4708 | |
| 4709 | depDerivDict = {} |
| 4710 | |
| 4711 | for j in dependentVars: |
| 4712 | |
| 4713 | depDerivDict[j] = np.zeros(shape=(len(x))) |
| 4714 | |
| 4715 | #print 'dependent vars',dependentVars |
| 4716 | |
| 4717 | lenX = len(x) |
| 4718 | |
| 4719 | hId = Histogram['hId'] |
| 4720 | |
| 4721 | hfx = ':%d:'%(hId) |
| 4722 | |
| 4723 | bakType = calcControls[hfx+'bakType'] |
| 4724 | |
| 4725 | dMdv = np.zeros(shape=(len(varylist),len(x))) |
| 4726 | |
| 4727 | dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x) |
| 4728 | |
| 4729 | if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints |
| 4730 | |
| 4731 | bBpos =varylist.index(hfx+'Back:0') |
| 4732 | |
| 4733 | dMdv[bBpos:bBpos+len(dMdb)] = dMdb |
| 4734 | |
| 4735 | names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU'] |
| 4736 | |
| 4737 | for name in varylist: |
| 4738 | |
| 4739 | if 'Debye' in name: |
| 4740 | |
| 4741 | id = int(name.split(':')[-1]) |
| 4742 | |
| 4743 | parm = name[:int(name.rindex(':'))] |
| 4744 | |
| 4745 | ip = names.index(parm) |
| 4746 | |
| 4747 | dMdv[varylist.index(name)] = dMddb[3*id+ip] |
| 4748 | |
| 4749 | if 'C' in calcControls[hfx+'histType']: |
| 4750 | |
| 4751 | dx = x[1]-x[0] |
| 4752 | |
| 4753 | shl = max(parmDict[hfx+'SH/L'],0.002) |
| 4754 | |
| 4755 | Ka2 = False |
| 4756 | |
| 4757 | if hfx+'Lam1' in parmDict.keys(): |
| 4758 | |
| 4759 | wave = parmDict[hfx+'Lam1'] |
| 4760 | |
| 4761 | Ka2 = True |
| 4762 | |
| 4763 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
| 4764 | |
| 4765 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
| 4766 | |
| 4767 | else: |
| 4768 | |
| 4769 | wave = parmDict[hfx+'Lam'] |
| 4770 | |
| 4771 | else: |
| 4772 | |
| 4773 | print 'TOF Undefined at present' |
| 4774 | |
| 4775 | raise ValueError |
| 4776 | |
| 4777 | for phase in Histogram['Reflection Lists']: |
| 4778 | |
| 4779 | refList = Histogram['Reflection Lists'][phase] |
| 4780 | |
| 4781 | Phase = Phases[phase] |
| 4782 | |
| 4783 | SGData = Phase['General']['SGData'] |
| 4784 | |
| 4785 | pId = Phase['pId'] |
| 4786 | |
| 4787 | pfx = '%d::'%(pId) |
| 4788 | |
| 4789 | phfx = '%d:%d:'%(pId,hId) |
| 4790 | |
| 4791 | A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] |
| 4792 | |
| 4793 | G,g = G2lat.A2Gmat(A) #recip & real metric tensors |
| 4794 | |
| 4795 | GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies |
| 4796 | |
| 4797 | if 'Pawley' not in Phase['General']['Type']: |
| 4798 | |
| 4799 | dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict) |
| 4800 | |
| 4801 | for iref,refl in enumerate(refList): |
| 4802 | |
| 4803 | if 'C' in calcControls[hfx+'histType']: #CW powder |
| 4804 | |
| 4805 | h,k,l = refl[:3] |
| 4806 | |
| 4807 | dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) |
| 4808 | |
| 4809 | if 'Pawley' in Phase['General']['Type']: |
| 4810 | |
| 4811 | try: |
| 4812 | |
| 4813 | refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) |
| 4814 | |
| 4815 | except KeyError: |
| 4816 | |
| 4817 | # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) |
| 4818 | |
| 4819 | continue |
| 4820 | |
| 4821 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
| 4822 | |
| 4823 | iBeg = np.searchsorted(x,refl[5]-fmin) |
| 4824 | |
| 4825 | iFin = np.searchsorted(x,refl[5]+fmax) |
| 4826 | |
| 4827 | if not iBeg+iFin: #peak below low limit - skip peak |
| 4828 | |
| 4829 | continue |
| 4830 | |
| 4831 | elif not iBeg-iFin: #peak above high limit - done |
| 4832 | |
| 4833 | break |
| 4834 | |
| 4835 | pos = refl[5] |
| 4836 | |
| 4837 | tanth = tand(pos/2.0) |
| 4838 | |
| 4839 | costh = cosd(pos/2.0) |
| 4840 | |
| 4841 | dMdpk = np.zeros(shape=(6,len(x))) |
| 4842 | |
| 4843 | dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) |
| 4844 | |
| 4845 | for i in range(1,5): |
| 4846 | |
| 4847 | dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i] |
| 4848 | |
| 4849 | dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0] |
| 4850 | |
| 4851 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} |
| 4852 | |
| 4853 | if Ka2: |
| 4854 | |
| 4855 | pos2 = refl[5]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) |
| 4856 | |
| 4857 | kdelt = int((pos2-refl[5])/dx) |
| 4858 | |
| 4859 | iBeg2 = min(lenX,iBeg+kdelt) |
| 4860 | |
| 4861 | iFin2 = min(lenX,iFin+kdelt) |
| 4862 | |
| 4863 | if iBeg2-iFin2: |
| 4864 | |
| 4865 | dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) |
| 4866 | |
| 4867 | for i in range(1,5): |
| 4868 | |
| 4869 | dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i] |
| 4870 | |
| 4871 | dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0] |
| 4872 | |
| 4873 | dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0] |
| 4874 | |
| 4875 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]} |
| 4876 | |
| 4877 | if 'Pawley' in Phase['General']['Type']: |
| 4878 | |
| 4879 | try: |
| 4880 | |
| 4881 | idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])) |
| 4882 | |
| 4883 | dMdv[idx] = dervDict['int']/refl[9] |
| 4884 | |
| 4885 | # Assuming Pawley variables not in constraints |
| 4886 | |
| 4887 | except ValueError: |
| 4888 | |
| 4889 | pass |
| 4890 | |
| 4891 | dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) |
| 4892 | |
| 4893 | names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], |
| 4894 | |
| 4895 | hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], |
| 4896 | |
| 4897 | hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], |
| 4898 | |
| 4899 | hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], |
| 4900 | |
| 4901 | hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], |
| 4902 | |
| 4903 | hfx+'DisplaceY':[dpdY,'pos'],} |
| 4904 | |
| 4905 | for name in names: |
| 4906 | |
| 4907 | item = names[name] |
| 4908 | |
| 4909 | if name in varylist: |
| 4910 | |
| 4911 | dMdv[varylist.index(name)] += item[0]*dervDict[item[1]] |
| 4912 | |
| 4913 | elif name in dependentVars: |
| 4914 | |
| 4915 | depDerivDict[name] += item[0]*dervDict[item[1]] |
| 4916 | |
| 4917 | |
| 4918 | |
| 4919 | for iPO in dIdPO: |
| 4920 | |
| 4921 | if iPO in varylist: |
| 4922 | |
| 4923 | dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int'] |
| 4924 | |
| 4925 | elif iPO in dependentVars: |
| 4926 | |
| 4927 | depDerivDict[iPO] = dIdPO[iPO]*dervDict['int'] |
| 4928 | |
| 4929 | |
| 4930 | |
| 4931 | for i,name in enumerate(['omega','chi','phi']): |
| 4932 | |
| 4933 | aname = pfx+'SH '+name |
| 4934 | |
| 4935 | if aname in varylist: |
| 4936 | |
| 4937 | dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int'] |
| 4938 | |
| 4939 | elif aname in dependentVars: |
| 4940 | |
| 4941 | depDerivDict[aname] += dFdSA[i]*dervDict['int'] |
| 4942 | |
| 4943 | for iSH in dFdODF: |
| 4944 | |
| 4945 | if iSH in varylist: |
| 4946 | |
| 4947 | dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int'] |
| 4948 | |
| 4949 | elif iSH in dependentVars: |
| 4950 | |
| 4951 | depDerivDict[iSH] += dFdODF[iSH]*dervDict['int'] |
| 4952 | |
| 4953 | cellDervNames = cellVaryDerv(pfx,SGData,dpdA) |
| 4954 | |
| 4955 | for name,dpdA in cellDervNames: |
| 4956 | |
| 4957 | if name in varylist: |
| 4958 | |
| 4959 | dMdv[varylist.index(name)] += dpdA*dervDict['pos'] |
| 4960 | |
| 4961 | elif name in dependentVars: |
| 4962 | |
| 4963 | depDerivDict[name] += dpdA*dervDict['pos'] |
| 4964 | |
| 4965 | dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) |
| 4966 | |
| 4967 | for name in dDijDict: |
| 4968 | |
| 4969 | if name in varylist: |
| 4970 | |
| 4971 | dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos'] |
| 4972 | |
| 4973 | elif name in dependentVars: |
| 4974 | |
| 4975 | depDerivDict[name] += dDijDict[name]*dervDict['pos'] |
| 4976 | |
| 4977 | gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) |
| 4978 | |
| 4979 | for name in gamDict: |
| 4980 | |
| 4981 | if name in varylist: |
| 4982 | |
| 4983 | dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam'] |
| 4984 | |
| 4985 | elif name in dependentVars: |
| 4986 | |
| 4987 | depDerivDict[name] += gamDict[name]*dervDict['gam'] |
| 4988 | |
| 4989 | |
| 4990 | |
| 4991 | elif 'T' in calcControls[hfx+'histType']: |
| 4992 | |
| 4993 | print 'TOF Undefined at present' |
| 4994 | |
| 4995 | raise Exception #no TOF yet |
| 4996 | |
| 4997 | #do atom derivatives - for F,X & U so far |
| 4998 | |
| 4999 | corr = dervDict['int']/refl[9] |
| 5000 | |
| 5001 | for name in varylist+dependentVars: |
| 5002 | |
| 5003 | try: |
| 5004 | |
| 5005 | aname = name.split(pfx)[1][:2] |
| 5006 | |
| 5007 | if aname not in ['Af','dA','AU']: continue # skip anything not an atom param |
| 5008 | |
| 5009 | except IndexError: |
| 5010 | |
| 5011 | continue |
| 5012 | |
| 5013 | if name in varylist: |
| 5014 | |
| 5015 | dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr |
| 5016 | |
| 5017 | elif name in dependentVars: |
| 5018 | |
| 5019 | depDerivDict[name] += dFdvDict[name][iref]*corr |
| 5020 | |
| 5021 | # now process derivatives in constraints |
| 5022 | |
| 5023 | G2mv.Dict2Deriv(varylist,depDerivDict,dMdv) |
| 5024 | |
| 5025 | return dMdv |
| 5026 | |
| 5027 | |
| 5028 | |
| 5029 | def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 5030 | |
| 5031 | parmdict.update(zip(varylist,values)) |
| 5032 | |
| 5033 | G2mv.Dict2Map(parmdict,varylist) |
| 5034 | |
| 5035 | Histograms,Phases = HistoPhases |
| 5036 | |
| 5037 | nvar = len(varylist) |
| 5038 | |
| 5039 | dMdv = np.empty(0) |
| 5040 | |
| 5041 | histoList = Histograms.keys() |
| 5042 | |
| 5043 | histoList.sort() |
| 5044 | |
| 5045 | for histogram in histoList: |
| 5046 | |
| 5047 | if 'PWDR' in histogram[:4]: |
| 5048 | |
| 5049 | Histogram = Histograms[histogram] |
| 5050 | |
| 5051 | hId = Histogram['hId'] |
| 5052 | |
| 5053 | hfx = ':%d:'%(hId) |
| 5054 | |
| 5055 | Limits = calcControls[hfx+'Limits'] |
| 5056 | |
| 5057 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 5058 | |
| 5059 | xB = np.searchsorted(x,Limits[0]) |
| 5060 | |
| 5061 | xF = np.searchsorted(x,Limits[1]) |
| 5062 | |
| 5063 | dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], |
| 5064 | |
| 5065 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
| 5066 | |
| 5067 | if len(dMdv): |
| 5068 | |
| 5069 | dMdv = np.concatenate((dMdv.T,dMdvh.T)).T |
| 5070 | |
| 5071 | else: |
| 5072 | |
| 5073 | dMdv = dMdvh |
| 5074 | |
| 5075 | return dMdv |
| 5076 | |
| 5077 | |
| 5078 | |
| 5079 | def ComputePowderHessian(args): |
| 5080 | |
| 5081 | Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args |
| 5082 | |
| 5083 | hId = Histogram['hId'] |
| 5084 | |
| 5085 | hfx = ':%d:'%(hId) |
| 5086 | |
| 5087 | Limits = calcControls[hfx+'Limits'] |
| 5088 | |
| 5089 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 5090 | |
| 5091 | dy = y-yc |
| 5092 | |
| 5093 | xB = np.searchsorted(x,Limits[0]) |
| 5094 | |
| 5095 | xF = np.searchsorted(x,Limits[1]) |
| 5096 | |
| 5097 | dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv( |
| 5098 | |
| 5099 | parmdict,x[xB:xF], |
| 5100 | |
| 5101 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
| 5102 | |
| 5103 | Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) |
| 5104 | |
| 5105 | Hess = np.inner(dMdvh,dMdvh) |
| 5106 | |
| 5107 | return Vec,Hess |
| 5108 | |
| 5109 | |
| 5110 | |
| 5111 | def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 5112 | |
| 5113 | parmdict.update(zip(varylist,values)) |
| 5114 | |
| 5115 | G2mv.Dict2Map(parmdict,varylist) |
| 5116 | |
| 5117 | Histograms,Phases = HistoPhases |
| 5118 | |
| 5119 | nvar = len(varylist) |
| 5120 | |
| 5121 | HessSum = np.zeros((nvar,nvar)) |
| 5122 | |
| 5123 | VecSum = np.zeros(nvar) |
| 5124 | |
| 5125 | histoList = Histograms.keys() |
| 5126 | |
| 5127 | histoList.sort() |
| 5128 | |
| 5129 | argList = [] |
| 5130 | |
| 5131 | for histogram in histoList: |
| 5132 | |
| 5133 | if 'PWDR' in histogram[:4]: |
| 5134 | |
| 5135 | Histogram = Histograms[histogram] |
| 5136 | |
| 5137 | argList.append([ |
| 5138 | |
| 5139 | Histogram,parmdict,varylist,Phases, |
| 5140 | |
| 5141 | calcControls,pawleyLookup]) |
| 5142 | |
| 5143 | if MaxProcess > 1: |
| 5144 | |
| 5145 | mpPool = mp.Pool(processes=MaxProcess) |
| 5146 | |
| 5147 | #results = mpPool.map(ComputePowderHessian,argList) |
| 5148 | |
| 5149 | #for Vec,Hess in results: |
| 5150 | |
| 5151 | for Vec,Hess in mpPool.imap_unordered(ComputePowderHessian,argList): |
| 5152 | |
| 5153 | VecSum += Vec |
| 5154 | |
| 5155 | HessSum += Hess |
| 5156 | |
| 5157 | else: |
| 5158 | |
| 5159 | for arg in argList: |
| 5160 | |
| 5161 | Vec,Hess = ComputePowderHessian(arg) |
| 5162 | |
| 5163 | VecSum += Vec |
| 5164 | |
| 5165 | HessSum += Hess |
| 5166 | |
| 5167 | return VecSum,HessSum |
| 5168 | |
| 5169 | |
| 5170 | |
| 5171 | def ComputePowderProfile(args): |
| 5172 | |
| 5173 | Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args |
| 5174 | |
| 5175 | hId = Histogram['hId'] |
| 5176 | |
| 5177 | print 'ComputePowderProfile',hId |
| 5178 | |
| 5179 | hfx = ':%d:'%(hId) |
| 5180 | |
| 5181 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 5182 | |
| 5183 | Limits = calcControls[hfx+'Limits'] |
| 5184 | |
| 5185 | xB = np.searchsorted(x,Limits[0]) |
| 5186 | |
| 5187 | xF = np.searchsorted(x,Limits[1]) |
| 5188 | |
| 5189 | yc,yb = getPowderProfile(parmdict,x[xB:xF], |
| 5190 | |
| 5191 | varylist,Histogram,Phases,calcControls, |
| 5192 | |
| 5193 | pawleyLookup) |
| 5194 | |
| 5195 | return xB,xF,yc,yb,Histogram['Reflection Lists'] |
| 5196 | |
| 5197 | |
| 5198 | |
| 5199 | def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
| 5200 | |
| 5201 | parmdict.update(zip(varylist,values)) |
| 5202 | |
| 5203 | Values2Dict(parmdict, varylist, values) |
| 5204 | |
| 5205 | G2mv.Dict2Map(parmdict,varylist) |
| 5206 | |
| 5207 | Histograms,Phases = HistoPhases |
| 5208 | |
| 5209 | M = np.empty(0) |
| 5210 | |
| 5211 | sumwYo = 0 |
| 5212 | |
| 5213 | Nobs = 0 |
| 5214 | |
| 5215 | histoList = Histograms.keys() |
| 5216 | |
| 5217 | histoList.sort() |
| 5218 | |
| 5219 | argList = [] |
| 5220 | |
| 5221 | for histogram in histoList: |
| 5222 | |
| 5223 | if 'PWDR' in histogram[:4]: |
| 5224 | |
| 5225 | Histogram = Histograms[histogram] |
| 5226 | |
| 5227 | argList.append( |
| 5228 | |
| 5229 | [Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup] |
| 5230 | |
| 5231 | ) |
| 5232 | |
| 5233 | if MaxProcess > 1: |
| 5234 | |
| 5235 | mpPool = mp.Pool(processes=MaxProcess) |
| 5236 | |
| 5237 | #results = mpPool.map(ComputePowderProfile,argList) |
| 5238 | |
| 5239 | #for arg,res in zip(argList,results): |
| 5240 | |
| 5241 | #results = mpPool.map(ComputePowderProfile,argList) |
| 5242 | |
| 5243 | #for i,res in enumerate(results): |
| 5244 | |
| 5245 | for i,res in enumerate( |
| 5246 | |
| 5247 | mpPool.imap_unordered(ComputePowderProfile,argList) |
| 5248 | |
| 5249 | ): |
| 5250 | |
| 5251 | print 'process',i |
| 5252 | |
| 5253 | xB,xF,ycSect,ybSect,RL = res |
| 5254 | |
| 5255 | Histogram = argList[i][0] |
| 5256 | |
| 5257 | Histogram['Reflection Lists'] = RL |
| 5258 | |
| 5259 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 5260 | |
| 5261 | yc *= 0.0 #zero full calcd profiles |
| 5262 | |
| 5263 | yb *= 0.0 |
| 5264 | |
| 5265 | yd *= 0.0 |
| 5266 | |
| 5267 | Histogram['Nobs'] = xF-xB |
| 5268 | |
| 5269 | Nobs += Histogram['Nobs'] |
| 5270 | |
| 5271 | Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) |
| 5272 | |
| 5273 | sumwYo += Histogram['sumwYo'] |
| 5274 | |
| 5275 | |
| 5276 | |
| 5277 | yc[xB:xF] = ycSect |
| 5278 | |
| 5279 | yb[xB:xF] = ybSect |
| 5280 | |
| 5281 | yc[xB:xF] += yb[xB:xF] |
| 5282 | |
| 5283 | yd[xB:xF] = y[xB:xF]-yc[xB:xF] |
| 5284 | |
| 5285 | Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) |
| 5286 | |
| 5287 | wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF]) |
| 5288 | |
| 5289 | Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) |
| 5290 | |
| 5291 | M = np.concatenate((M,wdy)) |
| 5292 | |
| 5293 | else: |
| 5294 | |
| 5295 | for arg in argList: |
| 5296 | |
| 5297 | xB,xF,ycSect,ybSect,RL = ComputePowderProfile(arg) |
| 5298 | |
| 5299 | Histogram = arg[0] |
| 5300 | |
| 5301 | x,y,w,yc,yb,yd = Histogram['Data'] |
| 5302 | |
| 5303 | yc *= 0.0 #zero full calcd profiles |
| 5304 | |
| 5305 | yb *= 0.0 |
| 5306 | |
| 5307 | yd *= 0.0 |
| 5308 | |
| 5309 | Histogram['Nobs'] = xF-xB |
| 5310 | |
| 5311 | Nobs += Histogram['Nobs'] |
| 5312 | |
| 5313 | Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) |
| 5314 | |
| 5315 | sumwYo += Histogram['sumwYo'] |
| 5316 | |
| 5317 | |
| 5318 | |
| 5319 | yc[xB:xF] = ycSect |
| 5320 | |
| 5321 | yb[xB:xF] = ybSect |
| 5322 | |
| 5323 | yc[xB:xF] += yb[xB:xF] |
| 5324 | |
| 5325 | yd[xB:xF] = y[xB:xF]-yc[xB:xF] |
| 5326 | |
| 5327 | Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) |
| 5328 | |
| 5329 | wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF]) |
| 5330 | |
| 5331 | Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) |
| 5332 | |
| 5333 | if dlg: |
| 5334 | |
| 5335 | dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0] |
| 5336 | |
| 5337 | M = np.concatenate((M,wdy)) |
| 5338 | |
| 5339 | |
| 5340 | |
| 5341 | Histograms['sumwYo'] = sumwYo |
| 5342 | |
| 5343 | Histograms['Nobs'] = Nobs |
| 5344 | |
| 5345 | Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.) |
| 5346 | |
| 5347 | if dlg: |
| 5348 | |
| 5349 | GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0] |
| 5350 | |
| 5351 | if not GoOn: |
| 5352 | |
| 5353 | parmDict['saved values'] = values |
| 5354 | |
| 5355 | raise Exception #Abort!! |
| 5356 | |
| 5357 | return M |
| 5358 | |
| 5359 | |
| 5360 | |
| 5361 | |
| 5362 | |
| 5363 | def Refine(GPXfile,dlg): |
| 5364 | |
| 5365 | import pytexture as ptx |
| 5366 | |
| 5367 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
| 5368 | |
| 5369 | |
| 5370 | |
| 5371 | ShowBanner() |
| 5372 | |
| 5373 | varyList = [] |
| 5374 | |
| 5375 | parmDict = {} |
| 5376 | |
| 5377 | calcControls = {} |
| 5378 | |
| 5379 | G2mv.InitVars() |
| 5380 | |
| 5381 | Controls = GetControls(GPXfile) |
| 5382 | |
| 5383 | ShowControls(Controls) |
| 5384 | |
| 5385 | constrDict,constrFlag,fixedList = GetConstraints(GPXfile) |
| 5386 | |
| 5387 | Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) |
| 5388 | |
| 5389 | if not Phases: |
| 5390 | |
| 5391 | print ' *** ERROR - you have no histograms to refine! ***' |
| 5392 | |
| 5393 | print ' *** Refine aborted ***' |
| 5394 | |
| 5395 | raise Exception |
| 5396 | |
| 5397 | if not Histograms: |
| 5398 | |
| 5399 | print ' *** ERROR - you have no data to refine with! ***' |
| 5400 | |
| 5401 | print ' *** Refine aborted ***' |
| 5402 | |
| 5403 | raise Exception |
| 5404 | |
| 5405 | Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases) |
| 5406 | |
| 5407 | calcControls['Natoms'] = Natoms |
| 5408 | |
| 5409 | calcControls['FFtables'] = FFtables |
| 5410 | |
| 5411 | calcControls['BLtables'] = BLtables |
| 5412 | |
| 5413 | hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms) |
| 5414 | |
| 5415 | calcControls.update(controlDict) |
| 5416 | |
| 5417 | histVary,histDict,controlDict = GetHistogramData(Histograms) |
| 5418 | |
| 5419 | calcControls.update(controlDict) |
| 5420 | |
| 5421 | varyList = phaseVary+hapVary+histVary |
| 5422 | |
| 5423 | parmDict.update(phaseDict) |
| 5424 | |
| 5425 | parmDict.update(hapDict) |
| 5426 | |
| 5427 | parmDict.update(histDict) |
| 5428 | |
| 5429 | GetFprime(calcControls,Histograms) |
| 5430 | |
| 5431 | # do constraint processing |
| 5432 | |
| 5433 | try: |
| 5434 | |
| 5435 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
| 5436 | |
| 5437 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList) |
| 5438 | |
| 5439 | except: |
| 5440 | |
| 5441 | print ' *** ERROR - your constraints are internally inconsistent ***' |
| 5442 | |
| 5443 | raise Exception(' *** Refine aborted ***') |
| 5444 | |
| 5445 | # check to see which generated parameters are fully varied |
| 5446 | |
| 5447 | msg = G2mv.SetVaryFlags(varyList) |
| 5448 | |
| 5449 | if msg: |
| 5450 | |
| 5451 | print ' *** ERROR - you have not set the refine flags for constraints consistently! ***' |
| 5452 | |
| 5453 | print msg |
| 5454 | |
| 5455 | raise Exception(' *** Refine aborted ***') |
| 5456 | |
| 5457 | G2mv.Map2Dict(parmDict,varyList) |
| 5458 | |
| 5459 | # print G2mv.VarRemapShow(varyList) |
| 5460 | |
| 5461 | |
| 5462 | |
390 | | def ShowBanner(): |
391 | | print 80*'*' |
392 | | print ' General Structure Analysis System-II Crystal Structure Refinement' |
393 | | print ' by Robert B. Von Dreele & Brian H. Toby' |
394 | | print ' Argonne National Laboratory(C), 2010' |
395 | | print ' This product includes software developed by the UChicago Argonne, LLC,' |
396 | | print ' as Operator of Argonne National Laboratory.' |
397 | | print 80*'*','\n' |
398 | | |
399 | | def ShowControls(Controls): |
400 | | print ' Least squares controls:' |
401 | | print ' Refinement type: ',Controls['deriv type'] |
402 | | if 'Hessian' in Controls['deriv type']: |
403 | | print ' Maximum number of cycles:',Controls['max cyc'] |
| 5944 | |
| 5945 | def ShowBanner(name): |
| 5946 | |
| 5947 | print 80*'*' |
| 5948 | |
| 5949 | print ' Interatomic Distances and Angles for phase '+name |
| 5950 | |
| 5951 | print 80*'*','\n' |
| 5952 | |
| 5953 | |
| 5954 | |
| 5955 | ShowBanner(DisAglCtls['Name']) |
| 5956 | |
| 5957 | SGData = DisAglData['SGData'] |
| 5958 | |
| 5959 | SGtext = G2spc.SGPrint(SGData) |
| 5960 | |
| 5961 | for line in SGtext: print line |
| 5962 | |
| 5963 | Cell = DisAglData['Cell'] |
| 5964 | |
| 5965 | |
| 5966 | |
| 5967 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
| 5968 | |
| 5969 | covData = {} |
| 5970 | |
| 5971 | if 'covData' in DisAglData: |
| 5972 | |
| 5973 | covData = DisAglData['covData'] |
| 5974 | |
| 5975 | covMatrix = covData['covMatrix'] |
| 5976 | |
| 5977 | varyList = covData['varyList'] |
| 5978 | |
| 5979 | pfx = str(DisAglData['pId'])+'::' |
| 5980 | |
| 5981 | A = G2lat.cell2A(Cell[:6]) |
| 5982 | |
| 5983 | cellSig = getCellEsd(pfx,SGData,A,covData) |
| 5984 | |
| 5985 | names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = '] |
| 5986 | |
| 5987 | valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)] |
| 5988 | |
| 5989 | line = '\n Unit cell:' |
| 5990 | |
| 5991 | for name,vals in zip(names,valEsd): |
| 5992 | |
| 5993 | line += name+vals |
| 5994 | |
| 5995 | print line |
| 5996 | |
| 5997 | else: |
| 5998 | |
| 5999 | print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \ |
| 6000 | |
| 6001 | ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \ |
| 6002 | |
| 6003 | '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6]) |
| 6004 | |
| 6005 | Factor = DisAglCtls['Factors'] |
| 6006 | |
| 6007 | Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii']))) |
| 6008 | |
| 6009 | Units = np.array([ #is there a nicer way to make this? |
| 6010 | |
| 6011 | [-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], |
| 6012 | |
| 6013 | [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], |
| 6014 | |
| 6015 | [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]]) |
| 6016 | |
| 6017 | origAtoms = DisAglData['OrigAtoms'] |
| 6018 | |
| 6019 | targAtoms = DisAglData['TargAtoms'] |
| 6020 | |
| 6021 | for Oatom in origAtoms: |
| 6022 | |
| 6023 | OxyzNames = '' |
| 6024 | |
| 6025 | IndBlist = [] |
| 6026 | |
| 6027 | Dist = [] |
| 6028 | |
| 6029 | Vect = [] |
| 6030 | |
| 6031 | VectA = [] |
| 6032 | |
| 6033 | angles = [] |
| 6034 | |
| 6035 | for Tatom in targAtoms: |
| 6036 | |
| 6037 | Xvcov = [] |
| 6038 | |
| 6039 | TxyzNames = '' |
| 6040 | |
| 6041 | if 'covData' in DisAglData: |
| 6042 | |
| 6043 | OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])] |
| 6044 | |
| 6045 | TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])] |
| 6046 | |
| 6047 | Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix) |
| 6048 | |
| 6049 | result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False) |
| 6050 | |
| 6051 | BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0] |
| 6052 | |
| 6053 | AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1] |
| 6054 | |
| 6055 | for Txyz,Top,Tunit in result: |
| 6056 | |
| 6057 | Dx = (Txyz-np.array(Oatom[3:6]))+Units |
| 6058 | |
| 6059 | dx = np.inner(Amat,Dx) |
| 6060 | |
| 6061 | dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5) |
| 6062 | |
| 6063 | IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.)) |
| 6064 | |
| 6065 | if np.any(IndB): |
| 6066 | |
| 6067 | for indb in IndB: |
| 6068 | |
| 6069 | for i in range(len(indb)): |
| 6070 | |
| 6071 | if str(dx.T[indb][i]) not in IndBlist: |
| 6072 | |
| 6073 | IndBlist.append(str(dx.T[indb][i])) |
| 6074 | |
| 6075 | unit = Units[indb][i] |
| 6076 | |
| 6077 | tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2]) |
| 6078 | |
| 6079 | pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData) |
| 6080 | |
| 6081 | sig = 0.0 |
| 6082 | |
| 6083 | if len(Xvcov): |
| 6084 | |
| 6085 | sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx))) |
| 6086 | |
| 6087 | Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig]) |
| 6088 | |
| 6089 | if (Dist[-1][-1]-AsumR) <= 0.: |
| 6090 | |
| 6091 | Vect.append(dx.T[indb][i]/Dist[-1][-2]) |
| 6092 | |
| 6093 | VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top]) |
| 6094 | |
| 6095 | else: |
| 6096 | |
| 6097 | Vect.append([0.,0.,0.]) |
| 6098 | |
| 6099 | VectA.append([]) |
| 6100 | |
| 6101 | Vect = np.array(Vect) |
| 6102 | |
| 6103 | angles = np.zeros((len(Vect),len(Vect))) |
| 6104 | |
| 6105 | angsig = np.zeros((len(Vect),len(Vect))) |
| 6106 | |
| 6107 | for i,veca in enumerate(Vect): |
| 6108 | |
| 6109 | if np.any(veca): |
| 6110 | |
| 6111 | for j,vecb in enumerate(Vect): |
| 6112 | |
| 6113 | if np.any(vecb): |
| 6114 | |
| 6115 | angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData) |
| 6116 | |
| 6117 | line = '' |
| 6118 | |
| 6119 | for i,x in enumerate(Oatom[3:6]): |
| 6120 | |
| 6121 | if len(Xvcov): |
| 6122 | |
| 6123 | line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True)) |
| 6124 | |
| 6125 | else: |
| 6126 | |
| 6127 | line += '%12.5f'%(x) |
| 6128 | |
| 6129 | print '\n Distances & angles for ',Oatom[1],' at ',line |
| 6130 | |
| 6131 | print 80*'*' |
| 6132 | |
| 6133 | line = '' |
| 6134 | |
| 6135 | for dist in Dist[:-1]: |
| 6136 | |
| 6137 | line += '%12s'%(dist[1].center(12)) |
| 6138 | |
| 6139 | print ' To cell +(sym. op.) dist. ',line |
| 6140 | |
| 6141 | for i,dist in enumerate(Dist): |
| 6142 | |
| 6143 | line = '' |
| 6144 | |
| 6145 | for j,angle in enumerate(angles[i][0:i]): |
| 6146 | |
| 6147 | sig = angsig[i][j] |
| 6148 | |
| 6149 | if angle: |
| 6150 | |
| 6151 | if sig: |
| 6152 | |
| 6153 | line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12)) |
| 6154 | |
| 6155 | else: |
| 6156 | |
| 6157 | val = '%.3f'%(angle) |
| 6158 | |
| 6159 | line += '%12s'%(val.center(12)) |
| 6160 | |
| 6161 | else: |
| 6162 | |
| 6163 | line += 12*' ' |
| 6164 | |
| 6165 | if dist[5]: #sig exists! |
| 6166 | |
| 6167 | val = G2mth.ValEsd(dist[4],dist[5]) |
| 6168 | |
| 6169 | else: |
| 6170 | |
| 6171 | val = '%8.4f'%(dist[4]) |
| 6172 | |
| 6173 | print ' %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line |
| 6174 | |
| 6175 | |
| 6176 | |
| 6177 | def Torsion(TorsionData): |
| 6178 | |
| 6179 | |
| 6180 | |
| 6181 | def ShowBanner(name): |
| 6182 | |
| 6183 | print 80*'*' |
| 6184 | |
| 6185 | print ' Torsion angle for phase '+name |
| 6186 | |
| 6187 | print 80*'*' |
| 6188 | |
| 6189 | |
| 6190 | |
| 6191 | ShowBanner(TorsionData['Name']) |
| 6192 | |
| 6193 | SGData = TorsionData['SGData'] |
| 6194 | |
| 6195 | Cell = TorsionData['Cell'] |
| 6196 | |
| 6197 | |
| 6198 | |
| 6199 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
| 6200 | |
| 6201 | covData = {} |
| 6202 | |
| 6203 | pfx = '' |
| 6204 | |
| 6205 | if 'covData' in TorsionData: |
| 6206 | |
| 6207 | covData = TorsionData['covData'] |
| 6208 | |
| 6209 | covMatrix = covData['covMatrix'] |
| 6210 | |
| 6211 | varyList = covData['varyList'] |
| 6212 | |
| 6213 | pfx = str(TorsionData['pId'])+'::' |
| 6214 | |
| 6215 | #find one end of 4 atom string - involved in longest distance |
| 6216 | |
| 6217 | dist = {} |
| 6218 | |
| 6219 | for i,X1 in enumerate(TorsionData['Datoms']): |
| 6220 | |
| 6221 | for j,X2 in enumerate(TorsionData['Datoms'][:i]): |
| 6222 | |
| 6223 | dist[np.sqrt(np.sum(np.inner(Amat,np.array(X2[3:6])-np.array(X1[3:6]))**2))] = [i,j] |
| 6224 | |
| 6225 | sortdist = dist.keys() |
| 6226 | |
| 6227 | sortdist.sort() |
| 6228 | |
| 6229 | end = dist[sortdist[-1]][0] |
| 6230 | |
| 6231 | #order atoms in distance from end - defines sequence of atoms for the torsion angle |
| 6232 | |
| 6233 | dist = {} |
| 6234 | |
| 6235 | X1 = TorsionData['Datoms'][end] |
| 6236 | |
| 6237 | for i,X2 in enumerate(TorsionData['Datoms']): |
| 6238 | |
| 6239 | dist[np.sqrt(np.sum(np.inner(Amat,np.array(X2[3:6])-np.array(X1[3:6]))**2))] = i |
| 6240 | |
| 6241 | sortdist = dist.keys() |
| 6242 | |
| 6243 | sortdist.sort() |
| 6244 | |
| 6245 | Datoms = [] |
| 6246 | |
| 6247 | Oatoms = [] |
| 6248 | |
| 6249 | for d in sortdist: |
| 6250 | |
| 6251 | atom = TorsionData['Datoms'][dist[d]] |
| 6252 | |
| 6253 | symop = atom[-1].split('+') |
| 6254 | |
| 6255 | if len(symop) == 1: |
| 6256 | |
| 6257 | symop.append('0,0,0') |
| 6258 | |
| 6259 | symop[0] = int(symop[0]) |
| 6260 | |
| 6261 | symop[1] = eval(symop[1]) |
| 6262 | |
| 6263 | atom[-1] = symop |
| 6264 | |
| 6265 | Datoms.append(atom) |
| 6266 | |
| 6267 | oatom = TorsionData['Oatoms'][dist[d]] |
| 6268 | |
| 6269 | names = ['','',''] |
| 6270 | |
| 6271 | if pfx: |
| 6272 | |
| 6273 | names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])] |
| 6274 | |
| 6275 | oatom += [names,] |
| 6276 | |
| 6277 | Oatoms.append(oatom) |
| 6278 | |
| 6279 | Tors,sig = G2mth.GetTorsionSig(Oatoms,Datoms,Amat,SGData,covData) |
| 6280 | |
| 6281 | print ' Torsion angle for atom sequence: ',[Datoms[i][1] for i in range(4)],'=',G2mth.ValEsd(Tors,sig) |
| 6282 | |
| 6283 | print ' NB: Atom sequence determined by interatomic distances' |
| 6284 | |
| 6285 | |
| 6286 | |
| 6287 | def BestPlane(PlaneData): |
| 6288 | |
| 6289 | |
| 6290 | |
| 6291 | def ShowBanner(name): |
| 6292 | |
| 6293 | print 80*'*' |
| 6294 | |
| 6295 | print ' Best plane result for phase '+name |
| 6296 | |
| 6297 | print 80*'*','\n' |
| 6298 | |
| 6299 | |
| 6300 | |
| 6301 | ShowBanner(PlaneData['Name']) |
| 6302 | |
| 6303 | |
| 6304 | |
| 6305 | Cell = PlaneData['Cell'] |
| 6306 | |
| 6307 | Amat,Bmat = G2lat.cell2AB(Cell[:6]) |
| 6308 | |
| 6309 | Atoms = PlaneData['Atoms'] |
| 6310 | |
| 6311 | sumXYZ = np.zeros(3) |
| 6312 | |
| 6313 | XYZ = [] |
| 6314 | |
| 6315 | Natoms = len(Atoms) |
| 6316 | |
| 6317 | for atom in Atoms: |
| 6318 | |
| 6319 | xyz = np.array(atom[3:6]) |
| 6320 | |
| 6321 | XYZ.append(xyz) |
| 6322 | |
| 6323 | sumXYZ += xyz |
| 6324 | |
| 6325 | sumXYZ /= Natoms |
| 6326 | |
| 6327 | XYZ = np.array(XYZ)-sumXYZ |
| 6328 | |
| 6329 | XYZ = np.inner(Amat,XYZ).T |
| 6330 | |
| 6331 | Zmat = np.zeros((3,3)) |
| 6332 | |
| 6333 | for i,xyz in enumerate(XYZ): |
| 6334 | |
| 6335 | Zmat += np.outer(xyz.T,xyz) |
| 6336 | |
| 6337 | print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]) |
| 6338 | |
| 6339 | Evec,Emat = nl.eig(Zmat) |
| 6340 | |
| 6341 | Evec = np.sqrt(Evec)/(Natoms-3) |
| 6342 | |
| 6343 | Order = np.argsort(Evec) |
| 6344 | |
| 6345 | XYZ = np.inner(XYZ,Emat.T).T |
| 6346 | |
| 6347 | XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T |
| 6348 | |
| 6349 | print ' Atoms in Cartesian best plane coordinates:' |
| 6350 | |
| 6351 | print ' Name X Y Z' |
| 6352 | |
| 6353 | for i,xyz in enumerate(XYZ): |
| 6354 | |
| 6355 | print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]) |
| 6356 | |
| 6357 | print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]) |
| 6358 | |
| 6359 | |
| 6360 | |
| 6361 | |
| 6362 | |
| 6363 | def main(): |
| 6364 | |
| 6365 | arg = sys.argv |
| 6366 | |
| 6367 | try: |
| 6368 | |
| 6369 | mkl.set_num_threads(MaxThreads) |
| 6370 | |
| 6371 | print "Max threads ",mkl.get_max_threads() |
| 6372 | |
| 6373 | ActualThreads = mkl.get_max_threads() |
| 6374 | |
| 6375 | except: |
| 6376 | |
| 6377 | print "MKL module not present" |
| 6378 | |
| 6379 | ActualThreads = 1 |
| 6380 | |
| 6381 | if len(arg) > 1: |
| 6382 | |
| 6383 | GPXfile = arg[1] |
| 6384 | |
| 6385 | if not ospath.exists(GPXfile): |
| 6386 | |
| 6387 | print 'ERROR - ',GPXfile," doesn't exist!" |
| 6388 | |
| 6389 | exit() |
| 6390 | |
| 6391 | GPXpath = ospath.dirname(arg[1]) |
| 6392 | |
| 6393 | Refine(GPXfile,None) |
| 6394 | |