Changeset 3297


Ignore:
Timestamp:
Feb 27, 2018 4:25:17 PM (4 years ago)
Author:
vondreele
Message:

progress on magnetism - BNS system ok for monoclinic, orthorhombic, tetragonal, hexagonal & trigonal
a few exceptions remain (cubic, rhombohedral, etc.)

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIphsGUI.py

    r3296 r3297  
    118118                        tableSizer.Add(text,0,WACV)
    119119                text = wx.StaticText(self.panel,label=' (%s) '%(self.names[j]))
    120                 if self.spins[j+ic*lentable] < 0 or Red:
    121                     text.SetForegroundColour('Red')
     120                try:
     121                    if self.spins[j+ic*lentable] < 0 or Red:
     122                        text.SetForegroundColour('Red')
     123                except IndexError:
     124                    print(self.spins,j,ic,lentable,self.names[j])
    122125                tableSizer.Add(text,0,WACV)
    123126                if not j%2:
     
    440443        def OnMag(event):
    441444            self.ifMag = mag.GetValue()
     445            wx.CallAfter(self.Draw)
    442446           
    443447        def OnConstr(event):
    444448            self.ifConstr = constr.GetValue()
    445 
     449           
     450        def OnBNSlatt(event):
     451            Obj = event.GetEventObject()
     452            BNSlatt = Obj.GetValue()
     453            SGData['BNSlattsym'] = [BNSlatt,BNSsym[BNSlatt]]
     454            self.Trans = G2spc.ApplyBNSlatt(SGData,SGData['BNSlattsym'])
     455            wx.CallAfter(self.Draw)
     456
     457        SGData = self.Phase['General']['SGData']
    446458        self.panel.Destroy()
    447459        self.panel = wx.Panel(self)
     
    509521            constr.Bind(wx.EVT_CHECKBOX,OnConstr)
    510522            mainSizer.Add(constr,0,WACV)
    511 
     523        if self.ifMag:
     524            GenSym,GenFlg,BNSsym = G2spc.GetGenSym(SGData)
     525            BNSizer = wx.BoxSizer(wx.HORIZONTAL)
     526            BNSizer.Add(wx.StaticText(self.panel,label=' BNS lattice:'),0,WACV)
     527            BNS = wx.ComboBox(self.panel,value=SGData['BNSlattsym'][0],choices=list(BNSsym.keys()),style=wx.CB_READONLY|wx.CB_DROPDOWN)
     528            BNS.Bind(wx.EVT_COMBOBOX,OnBNSlatt)
     529            BNSizer.Add(BNS,0,WACV)
     530            mainSizer.Add(BNSizer,0,WACV)
    512531        TestBtn = wx.Button(self.panel,-1,"Test")
    513532        TestBtn.Bind(wx.EVT_BUTTON, OnTest)
     
    18201839                OprNames = G2spc.GenMagOps(SGData)[0]
    18211840            else:
    1822                 if not len(GenSym) or SGData['SGGray']:
     1841                if not len(GenSym) or SGData['SGGray'] or '(' in MagSym:    #not if BNS centered either
    18231842                    magSizer.Add(wx.StaticText(General,label=' No spin inversion allowed'),0,WACV)
    18241843                    OprNames,SpnFlp = G2spc.GenMagOps(SGData)                   
     
    22572276        if generalData['Type'] == 'magnetic':
    22582277            if not generalData['SGData']['SGFixed']:
    2259                 GenSym,GenFlg = G2spc.GetGenSym(generalData['SGData'])
     2278                GenSym,GenFlg,BNSsym = G2spc.GetGenSym(generalData['SGData'])
    22602279                generalData['SGData']['GenSym'] = GenSym
    22612280                generalData['SGData']['GenFlg'] = GenFlg
     
    23262345                finally:
    23272346                    dlg.Destroy()
    2328                 SGData['GenSym'],SGData['GenFlg'] = G2spc.GetGenSym(SGData)
     2347                SGData['GenSym'],SGData['GenFlg'],BNSsym = G2spc.GetGenSym(SGData)
    23292348                SGData['MagSpGrp'] = G2spc.MagSGSym(SGData)
    23302349                SGData['OprNames'],SGData['SpnFlp'] = G2spc.GenMagOps(SGData)
     
    27322751                    if not SGData['SGGray']:
    27332752                        CSI = G2spc.GetCSpqinel(SytSym,SpnFlp,dupDir)
    2734                     print (SytSym,Nop,SpnFlp[Nop],CSI,dupDir)
     2753#                    print (SytSym,Nop,SpnFlp[Nop],CSI,dupDir)
    27352754                    for i in range(3):
    27362755                        ci = i+colM
  • trunk/GSASIIspc.py

    r3295 r3297  
    6262             * 'SSGKl': default internal (Kl) part of supersymmetry point group; modified
    6363                in supersymmetry stuff depending on chosen modulation vector for Mono & Ortho
     64             * 'BNSlattsym': BNS lattice symbol & cenering op - used for magnetic structures
    6465
    6566    """
     
    8485    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
    8586    SGData['SGFixed'] = False
    86     if SGData['SGLatt'] == 'P':
    87         SGData['SGCen'] = np.array(([0,0,0],))
    88     elif SGData['SGLatt'] == 'A':
    89         SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
    90     elif SGData['SGLatt'] == 'B':
    91         SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
    92     elif SGData['SGLatt'] == 'C':
    93         SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
    94     elif SGData['SGLatt'] == 'I':
    95         SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
    96     elif SGData['SGLatt'] == 'F':
    97         SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
    98     elif SGData['SGLatt'] == 'R':
    99         SGData['SGCen'] = np.array(([0,0,0],[1./3.,2./3.,2./3.],[2./3.,1./3.,1./3.]))
    10087    SGData['SGOps'] = []
    10188    SGData['SGGen'] = []
    102     SGData['SGSpin'] = []
    10389    for i in range(SGInfo[5]):
    10490        Mat = np.array(SGInfo[6][i])
     
    10793        if 'array' in str(type(SGInfo[8])):        #patch for old fortran bin?
    10894            SGData['SGGen'].append(int(SGInfo[8][i]))
    109         SGData['SGSpin'].append(1)
    110     if SGData['SGLaue'] == '2/m' and SGData['SGLatt'] != 'P' and '/' in SGData['SpGrp']:
    111         SGData['SGSpin'].append(1)  #fix bug in fortran
    112     if 'F' in SGData['SpGrp']:
    113         SGData['SGSpin'] += [1,1,1,1]
    114     elif 'R' in SGData['SpGrp']:
    115         SGData['SGSpin'] += [1,1,1]
    116     elif SGData['SpGrp'][0] in ['A','B','C','I']:
    117         SGData['SGSpin'] += [1,]
     95    SGData['BNSlattsym'] = [LattSym[SGInfo[2]-1],[0,0,0]]
     96    lattSpin = []
     97    if SGData['SGLatt'] == 'P':
     98        SGData['SGCen'] = np.array(([0,0,0],))
     99    elif SGData['SGLatt'] == 'A':
     100        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
     101        lattSpin += [1,]
     102    elif SGData['SGLatt'] == 'B':
     103        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
     104        lattSpin += [1,]
     105    elif SGData['SGLatt'] == 'C':
     106        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
     107        lattSpin += [1,]
     108    elif SGData['SGLatt'] == 'I':
     109        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
     110        lattSpin += [1,]
     111    elif SGData['SGLatt'] == 'F':
     112        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
     113        lattSpin += [1,1,1,1]
     114    elif SGData['SGLatt'] == 'R':
     115        SGData['SGCen'] = np.array(([0,0,0],[1./3.,2./3.,2./3.],[2./3.,1./3.,1./3.]))
     116        lattSpin += [1,1,]
     117
     118#    if SGData['SGLaue'] == '2/m' and SGData['SGLatt'] != 'P' and '/' in SGData['SpGrp']:
     119#        SGData['SGSpin'].append(1)  #fix bug in fortran
     120#    if 'F' in SGData['SpGrp']:
     121#        SGData['SGSpin'] += [1,1,1,1]
     122#    elif 'R' in SGData['SpGrp']:
     123#        SGData['SGSpin'] += [1,1,1]
     124#    elif SGData['SpGrp'][0] in ['A','B','C','I']:
     125#        SGData['SGSpin'] += [1,]
     126
    118127    if SGData['SGInv']:
    119128        if SGData['SGLaue'] in ['-1','2/m','mmm']:
     
    196205    SGData['SGPolax'] = SGpolar(SGData)
    197206    SGData['SGPtGrp'],SGData['SSGKl'] = SGPtGroup(SGData)
     207
     208    if SGData['SGLatt'] == 'R':
     209        if SGData['SGPtGrp'] in ['3']:
     210            SGData['SGSpin'] = lattSpin+[1,]
     211        else:
     212            SGData['SGSpin'] = lattSpin+[1,1,1,]
     213           
     214         
     215    else:
     216        if SGData['SGPtGrp'] in ['1','3']:
     217            SGData['SGSpin'] = lattSpin
     218        elif SGData['SGPtGrp'] in ['-1','2','m','4','-4','-3','312','321','3m1','31m','6','-6','432','-43m',]:
     219            SGData['SGSpin'] = lattSpin+[1,]
     220        elif SGData['SGPtGrp'] in ['2/m','4/m','422','4mm','-42m','-4m2','-3m1','-31m',
     221            '6/m','622','6mm','-6m2','-62m',]:
     222            SGData['SGSpin'] = lattSpin+[1,1,]
     223        elif SGData['SGPtGrp'] in ['3',]:
     224            SGData['SGSpin'] = lattSpin+[1,1,1,1,]
     225        else: #'222'-'mmm','4/mmm','6/mmm'
     226            SGData['SGSpin'] = lattSpin+[1,1,1,]
     227   
    198228    return SGInfo[-1],SGData
    199229
     
    382412    SGText = []
    383413    SGText.append(' Space Group: '+SGData['SpGrp'])
    384     if SGData['SGGray']: SGText[-1] += " 1'"
     414    if SGData.get('SGGray',False): SGText[-1] += " 1'"
    385415    CentStr = 'centrosymmetric'
    386416    if not SGData['SGInv']:
     
    668698            UsymOp.append(' m110 ')
    669699            OprFlg.append(24)
    670     ncv = len(SGData['SGCen'])
    671     if ncv > 1:
    672         for icv in range(ncv):
    673             if SGData['SpGrp'] in ['F d d 2','F d 2 d','F 2 d d','F d d d']:
    674                 break
    675             if 'F' in SGData['SpGrp'] and SGData['SGSys'] == 'cubic':
    676                 break
    677             if icv:
    678                 if SGData['SGCen'][icv][0] == 0.5:
    679                     if SGData['SGCen'][icv][1] == 0.5:
    680                         if SGData['SGCen'][icv][2] == 0.5:
    681                             if not SGData['SpGrp'] in ['I 41/a','I 41 m d',
    682                                 'I 41 c d','I -4 2 d','I -4 3 d',
    683                                 'I a 3 d','I a -3 d','I b 3 d','I b -3 d']:
    684                                 UsymOp.append(' Icen ')
    685                         else:
    686                             UsymOp.append(' Ccen ')
    687                     else:
    688                         UsymOp.append(' Bcen ')
    689                 elif SGData['SGCen'][icv][1] == 0.5:
    690                     UsymOp.append(' Acen ')
    691     return UsymOp,OprFlg
     700           
     701    if 'P' in SGData['SGLatt']:
     702        if SGData['SGSys'] == 'triclinic':
     703            BNSsym = {'P':[0,0,0],'P(a)':[.5,0,0],'P(b)':[0,.5,0],'P(c)':[0,0,.5]}           
     704        elif SGData['SGSys'] == 'monoclinic':
     705            BNSsym = {'P':[0,0,0],'P(a)':[.5,0,0],'P(b)':[0,.5,0],'P(c)':[0,0,.5]}
     706            if SGData['SGUniq'] == 'a':
     707                BNSsym.update({'P(B)':[.5,0,.5],'P(C)':[.5,.5,0]})
     708            elif SGData['SGUniq'] == 'b':
     709                BNSsym.update({'P(A)':[.5,.5,0],'P(C)':[0,.5,.5]})
     710            elif SGData['SGUniq'] == 'c':
     711                BNSsym.update({'P(A)':[0,.5,.5],'P(B)':[.5,0,.5]})
     712        elif SGData['SGSys'] == 'orthorhombic':
     713            BNSsym = {'P':[0,0,0],'P(a)':[.5,0,0],'P(b)':[0,.5,0],'P(c)':[0,0,.5],
     714                'P(A)':[0,.5,.5],'P(B)':[.5,0,.5],'P(C)':[.5,.5,0],'P(I)':[.5,.5,.5]}
     715        elif SGData['SGSys'] == 'tetragonal':
     716            BNSsym = {'P':[0,0,0],'P(c)':[0,0,.5],'P(C)':[.5,.5,0],'P(I)':[.5,.5,.5]}           
     717        elif SGData['SGSys'] in ['trigonal','hexagonal']:
     718            BNSsym = {'P':[0,0,0],'P(c)':[0,0,.5]}           
     719        elif SGData['SGSys'] == 'cubic':
     720            BNSsym = {'P':[0,0,0],'P(I)':[.5,.5,.5]}           
     721           
     722    elif 'A' in SGData['SGLatt']:
     723        if SGData['SGSys'] == 'monoclinic':
     724            BNSsym = {'A':[0,0,0],}
     725            if SGData['SGUniq'] == 'b':
     726                BNSsym.update({'A(a)':[.5,0,0],'A(c)':[0,0,.5]})
     727            elif SGData['SGUniq'] == 'c':
     728                BNSsym.update({'A(a)':[.5,0,0],'A(b)':[0,.5,0]})
     729        elif SGData['SGSys'] == 'orthorhombic':
     730            BNSsym = {'A':[0,0,0],'A(a)':[.5,0,0],'A(b)':[0,.5,0],'A(c)':[0,0,.5],
     731               'A(B)':[.5,0,.5],'A(C)':[.5,.5,0]}   
     732           
     733    elif 'B' in SGData['SGLatt']:
     734        if SGData['SGSys'] == 'monoclinic':
     735            BNSsym = {'B':[0,0,0],}
     736            if SGData['SGUniq'] == 'a':
     737                BNSsym.update({'B(b)':[0,.5,0],'B(c)':[0,0,.5]})
     738            elif SGData['SGUniq'] == 'c':
     739                BNSsym.update({'B(a)':[.5,0,0],'B(b)':[0,.5,0]})
     740        elif SGData['SGSys'] == 'orthorhombic':
     741            BNSsym = {'B':[0,0,0],'B(a)':[.5,0,0],'B(b)':[0,.5,0],'B(c)':[0,0,.5],
     742                'B(A)':[0,.5,.5],'B(C)':[.5,.5,0]}     
     743           
     744    elif 'C' in SGData['SGLatt']:
     745        if SGData['SGSys'] == 'monoclinic':
     746            BNSsym = {'C':[0,0,0],}
     747            if SGData['SGUniq'] == 'a':
     748                BNSsym.update({'C(b)':[0,.5,.0],'C(c)':[0,0,.5]})
     749            elif SGData['SGUniq'] == 'b':
     750                BNSsym.update({'C(a)':[.5,0,0],'C(c)':[0,0,.5]})
     751        elif SGData['SGSys'] == 'orthorhombic':
     752            BNSsym = {'C':[0,0,0],'C(a)':[.5,0,0],'C(b)':[0,.5,0],'C(c)':[0,0,.5],
     753                'C(A)':[0,.5,.5],'C(B)':[.5,0,.5]}     
     754           
     755    elif 'I' in SGData['SGLatt']:
     756        if SGData['SGSys'] in ['monoclinic','orthorhombic']:
     757            BNSsym = {'I':[0,0,0],'I(a)':[.5,0,0],'I(b)':[0,.5,0],'I(c)':[0,0,.5]}
     758        elif SGData['SGSys'] == 'tetragonal':
     759            BNSsym = {'I':[0,0,0],'I(c)':[0,0,.5]}
     760        elif SGData['SGSys'] == 'cubic':
     761            BNSsym = {'I':[0,0,0]}
     762           
     763    elif 'F' in SGData['SGLatt']:
     764        if SGData['SGSys'] in ['monoclinic','orthorhombic','cubic']:
     765            BNSsym = {'F':[0,0,0],'F(S)':[.5,.5,.5]}
     766           
     767    elif 'R' in SGData['SGLatt']:
     768        BNSsym = {'R':[0,0,0],'R(I)':[0,0,.5]}
     769    return UsymOp,OprFlg,BNSsym
     770
     771def ApplyBNSlatt(SGData,BNSlatt):
     772    Tmat = np.eye(3)
     773    BNS = BNSlatt[0]
     774    A = np.array(BNSlatt[1])
     775    SGCen = SGData['SGCen']
     776    if '(a)' in BNS:
     777        Tmat[0,0] = 2.0
     778    elif '(b)' in BNS:
     779        Tmat[1,1] = 2.0
     780    elif '(c)' in BNS:
     781        Tmat[2,2] = 2.0
     782    elif '(A)' in BNS:
     783        Tmat[0,0] = 2.0
     784    elif '(B)' in BNS:
     785        Tmat[1,1] = 2.0
     786    elif '(C)' in BNS:
     787        Tmat[2,2] = 2.0
     788    elif '(I)' in BNS:
     789        Tmat *= 2.0
     790    elif '(S)' in BNS:
     791        SGData['SGSpin'][-1] = -1
     792        SGData['SGSpin'] += [-1,-1,-1,]
     793        Tmat *= 2.0
     794    SGData['SGSpin'].append(-1)
     795    if 'P' not in BNS:
     796        SGData['SGSpin'].append(-1)
     797    C = SGCen+A
     798    SGData['SGCen'] = np.vstack((SGCen,C))%1.
     799    return Tmat
    692800   
    693801def CheckSpin(isym,SGData):
     
    734842    SGLaue = SGData['SGLaue']
    735843    SpnFlp = SGData['SGSpin']
     844#    print('SpnFlp',SpnFlp)
    736845    GenSym = SGData['GenSym']
    737846    SGPtGrp = SGData['SGPtGrp']
     
    745854            magSym[1] += "'"
    746855            SGData['MagPtGp'] += "'"
    747         if magSym[0] in ['A','B','C','I'] and SGData['SpGrp'] != 'I 41/a':
    748             if SpnFlp[1] < 0:
    749                 magSym[0] += '(P)'
     856    elif SGLaue in ['2/m','4/m','6/m']: #all ok
     857        Uniq = {'a':1,'b':2,'c':3,'':1}
     858        id = [0,1]
     859        if len(magSym) > 2:
     860            id = [0,Uniq[SGData['SGUniq']]]
     861        sym = magSym[id[1]].split('/')
     862        Ptsym = SGLaue.split('/')
     863        if len(GenSym) == 3:
     864            for i in [0,1,2]:
     865                if SpnFlp[i] < 0:
     866                    sym[i] += "'"
     867                    Ptsym[i] += "'"
     868        else:
     869            for i in range(len(GenSym)):
     870                if SpnFlp[i] < 0:                     
     871                    sym[i] += "'"
     872                    Ptsym[i] += "'"
     873        SGData['MagPtGp'] = '/'.join(Ptsym)
     874        magSym[id[1]] = '/'.join(sym)
    750875    elif SGPtGrp in ['mmm','mm2','m2m','2mm','222']:
    751876        SGData['MagPtGp'] = ''
     
    755880                magSym[i+1] += "'"
    756881                SGData['MagPtGp'] += "'"
    757         if len(GenSym) > 3:
    758             if magSym[0] == 'F':
    759                 if SpnFlp[3]+SpnFlp[4]+SpnFlp[5] < 0:
    760                     if SpnFlp[3] > 0:
    761                         magSym[0] += '(A)'
    762                     elif SpnFlp[4] > 0:
    763                         magSym[0] += '(B)'
    764                     elif SpnFlp[5] > 0:
    765                         magSym[0] += '(C)'
    766             else:
    767                 if SpnFlp[3] < 0:
    768                     magSym[0] += '(P)'
    769882    elif SGLaue == '6/mmm': #ok
     883        magPtGp = list(SGPtGrp)
    770884        if len(GenSym) == 2:
    771             magPtGp = ['6','m','m']
    772885            for i in [0,1]:
    773886                if SpnFlp[i] < 0:
     
    796909        SGData['MagPtGp'] = ''.join(magPtGp)
    797910    elif SGLaue == '4/mmm':
     911        magPtGp = list(SGPtGrp)
    798912        if len(GenSym) == 2:
    799             magPtGp = ['4','m','m']
    800913            for i in [0,1]:
    801914                if SpnFlp[i] < 0:
     
    823936                magSym[1] = '/'.join(sym)
    824937                magPtGp[0] = '/'.join(Ptsym)
    825                 if SpnFlp[3] < 0:
    826                     magSym[0] += '(P)'
    827938            else:
    828939                for i in [0,1]:
     
    831942                if SpnFlp[0]*SpnFlp[1] < 0:
    832943                    magSym[1] += "'"
    833                 if SpnFlp[2] < 0:
    834                     magSym[0] += '(P)'
    835944        SGData['MagPtGp'] = ''.join(magPtGp)
    836     elif SGLaue in ['2/m','4/m','6/m']: #all ok
    837         Uniq = {'a':1,'b':2,'c':3,'':1}
    838         id = [0,1]
    839         if len(magSym) > 2:
    840             id = [0,Uniq[SGData['SGUniq']]]
    841         sym = magSym[id[1]].split('/')
    842         Ptsym = SGLaue.split('/')
    843         if len(GenSym) == 3:
    844             for i in [0,1,2]:
    845                 if SpnFlp[i] < 0:
    846                     if i == 2:
    847                         magSym[0] += '(P)'
    848                     else:
    849                         sym[i] += "'"
    850                         Ptsym[i] += "'"
    851         else:
    852             for i in range(len(GenSym)):
    853                 if SpnFlp[i] < 0:                     
    854                     if i and magSym[0] in ['A','B','C','I'] and SGData['SpGrp'] != 'I 41/a':
    855                         magSym[0] += '(P)'
    856                     else:
    857                         sym[i] += "'"
    858                         Ptsym[i] += "'"
    859         SGData['MagPtGp'] = '/'.join(Ptsym)
    860         magSym[id[1]] = '/'.join(sym)
    861945    elif SGLaue in ['3','3m1','31m']:   #ok
    862 #        GSASIIpath.IPyBreak()
    863946        Ptsym = list(SGLaue)
    864947        if len(GenSym) == 1:    #all ok
     
    9181001    elif SGData['SGPtGrp'] == '23' and len(magSym):
    9191002        SGData['MagPtGp'] = '23'
    920         if SpnFlp[0] < 0:
    921             magSym[0] += '(P)'
    9221003    elif SGData['SGPtGrp'] == 'm3':
    9231004        SGData['MagPtGp'] = "m3"
     
    9271008            SGData['MagPtGp'] = "m'3'"
    9281009        if SpnFlp[1] < 0:
    929             magSym[0] += '(P)'
    9301010            if not 'm' in magSym[1]:    #only Ia3
    9311011                magSym[1].strip("'")
     
    9381018            magSym[3] += "'"
    9391019            Ptsym[1] += "'"
    940         if SpnFlp[1] < 0:
    941             magSym[0] += '(P)'
    9421020        SGData['MagPtGp'] = '3'.join(Ptsym)
    9431021    elif SGData['SGPtGrp'] == 'm-3m':
     
    9581036            magSym[3] += "'"
    9591037            Ptsym[2] += "'"
    960         if SpnFlp[2] < 0:
    961             magSym[0] += '(P)'
    9621038        SGData['MagPtGp'] = ''.join(Ptsym)
    9631039#    print SpnFlp
     1040    magSym[0] = SGData.get('BNSlattsym',[SGData['SGLatt'],[0,0,0]])[0]
    9641041    return ' '.join(magSym)
    9651042   
     
    21022179        else:
    21032180            csi = CSxinel[indx[3]]  #Q
    2104         print(opr,SpnFlp[dupDir[opr]],indx,csi,CSI)
     2181#        print(opr,SpnFlp[dupDir[opr]],indx,csi,CSI)
    21052182        if not len(csi):
    21062183            return [[0,0,0],[0.,0.,0.]]
  • trunk/GSASIIstrMath.py

    r3268 r3297  
    689689    if parmDict[pfx+'isMag']:       #TODO: fix the math - mag moments now along crystal axes
    690690        Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
    691 #        Mag = np.sqrt(np.inner(np.inner(Gdata.T,G),Gdata.T))
    692691        Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
    693 #        Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
    694692        Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
    695693        if SGData['SGInv'] and not SGData['SGFixed']:
     
    697695        Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
    698696        Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
    699 #        Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
    700 #        Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
    701697        Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
    702698        if SGData['SGInv'] and not SGData['SGFixed']:
     
    13081304    SGInv = SGData['SGInv']
    13091305    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1306    Ncen = len(SGData['SGCen'])
     1307    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
    13101308    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    13111309    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     
    13131311    BLtables = calcControls['BLtables']
    13141312    MFtables = calcControls['MFtables']
     1313    Amat,Bmat = G2lat.Gmat2AB(G)
    13151314    Flack = 1.0
    13161315    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     
    13201319    if not Xdata.size:          #no atoms in phase!
    13211320        return
     1321
     1322    if parmDict[pfx+'isMag']:       #TODO: fix the math - mag moments now along crystal axes
     1323        Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
     1324        Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
     1325        Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     1326        if SGData['SGInv'] and not SGData['SGFixed']:
     1327            Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     1328        Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
     1329        Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
     1330        Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
     1331        if SGData['SGInv'] and not SGData['SGFixed']:
     1332            Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
     1333
     1334
    13221335    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    13231336    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
     
    13831396        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    13841397        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
    1385         if 'T' in calcControls[hfx+'histType']:
    1386             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    1387             fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1398
     1399        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: math here??
     1400            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
     1401            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
     1402            if SGData['SGInv'] and not SGData['SGFixed']:
     1403                mphase = np.hstack((phase,-phase))
     1404            else:
     1405                mphase = phase
     1406            mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
     1407            mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
     1408            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
     1409            cosm = np.cos(mphase)                               #ditto
     1410            HM = np.inner(Bmat.T,H)                             #put into cartesian space
     1411            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
     1412            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
     1413            Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     1414            fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     1415            fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     1416            fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
     1417            fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
     1418            refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     1419            refl.T[7] = np.copy(refl.T[9])               
     1420            refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     1421
     1422
    13881423        else:
    1389             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    1390             fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1391         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
    1392         fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
    1393         fbg = fb*GfpuA[0]+fa*GfpuA[1]
    1394         fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
    1395         fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
    1396 #        GSASIIpath.IPyBreak()
    1397         if 'P' in calcControls[hfx+'histType']:
    1398             refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1399 #            refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
    1400             refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    1401         else:
    1402             refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
    1403             refl.T[8] = np.copy(refl.T[10])               
    1404             refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1424            if 'T' in calcControls[hfx+'histType']:
     1425                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     1426                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1427            else:
     1428                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     1429                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1430            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
     1431            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
     1432            fbg = fb*GfpuA[0]+fa*GfpuA[1]
     1433            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
     1434            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
     1435            if 'P' in calcControls[hfx+'histType']:
     1436                refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
     1437                refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1438            else:
     1439                refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
     1440                refl.T[8] = np.copy(refl.T[10])               
     1441                refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    14051442        iBeg += blkSize
    14061443    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
Note: See TracChangeset for help on using the changeset viewer.