Changeset 5663


Ignore:
Timestamp:
Sep 29, 2023 3:47:55 PM (21 months ago)
Author:
vondreele
Message:

Implement 1st draft of atom deformation refinement - functions from Phil Coppens. Functions apparently correct; derivatives need work. GUI interface done, constraints available for all deformation parameters.
New ORBtables available for orbital form factor coefficients
Add factor sizer to bond (only) control - was only for bond & angle control
new routine GetHistogramTypes? - from tree in G2dataGUI

Add optional comment line to G2ScrolledGrid - shows below table
Move Parameter Impact menu item to be under Compute partials

Location:
trunk
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified trunk/GSASIIElem.py

    r5572 r5663  
    1717import GSASIIpath
    1818GSASIIpath.SetVersionNumber("$Revision$")
     19import copy
    1920import numpy as np
    2021import atmdata
    … …  
    100101                FFtable[El] = item
    101102    return FFtable
     103   
     104def GetORBtable(atomTypes):
     105    ''' returns a dictionary of orbital form factor data for atom types found in atomTypes
     106
     107    :param list atomTypes: list of atom types
     108    :return: ORBtable, dictionary of orbital form factor data; key is atom type
     109
     110    '''
     111    ORBtable = {}
     112    for El in atomTypes:
     113        ORBtable[El] = copy.deepcopy(atmdata.OrbFF[El])
     114    return ORBtable
    102115   
    103116def GetMFtable(atomTypes,Landeg):
    … …  
    422435    t = -fb[:,np.newaxis]*SQ
    423436    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El.get('fc',0.0)
     437
     438def ScatFacDer(El, SQ):
     439    """compute derivative of form factor wrt SQ
     440
     441    :param El: element dictionary defined in GetFormFactorCoeff
     442    :param SQ: (sin-theta/lambda)**2
     443    :return: real part of form factor
     444    """
     445    fa = np.array(El['fa'])
     446    fb = np.array(El['fb'])
     447    t = -fb[:,np.newaxis]*SQ
     448    return -np.sum(fa[:,np.newaxis]*fb[:,np.newaxis]*np.exp(t)[:],axis=0)
    424449       
    425450def MagScatFac(El, SQ):
    … …  
    443468    MF0 = MMF0+(2.0/El['gfac']-1.0)*NMF0
    444469    return (MMF+(2.0/El['gfac']-1.0)*NMF)/MF0
     470
     471def ClosedFormFF(El,SQ,k,N):
     472    """Closed form expressions for FT Slater fxns. IT B Table 1.2.7.4
     473    (not used at present - doesn't make sense yet)
     474   
     475    :param El: element dictionary defined in GetFormFactorCoeff
     476    :param SQ: (sin-theta/lambda)**2
     477    :param k: int principal Bessel fxn order as in <jk>
     478    :param N: int power
     479   
     480    return: form factor
     481    """
     482    Z = El['Z']
     483    Z2 = Z**2
     484    K2 = 16.0*SQ*np.pi**2
     485    K2pZ2 = K2+Z2
     486    K = np.sqrt(K2)
     487    if not k: #==0
     488        if N == 1:
     489            return 1.0/K2pZ2
     490        elif N == 2:
     491            return 2.0*Z/K2pZ2**2
     492        elif N == 3:
     493            return 2.0*(3.0*Z2-K2)/K2pZ2**3
     494        elif N == 4:
     495            return 24.0*Z*(Z2-K2)/K2pZ2**4
     496        elif N == 5:
     497            return 24.0*(5.0*Z2-10.0*K2*Z2+K2**2)/K2pZ2**5
     498        elif N == 6:
     499            return 240.0*(K2-3.0*Z2)*(3.0*K2-Z2)/K2pZ2**6
     500        elif N == 7:
     501            return 720.0*(7.0*Z2**3-35.0*K2*Z2**2+21.0*Z2*K2**2-K2**3)/K2pZ2**7
     502        elif N == 8:
     503            return 40320.0*(Z*Z2**3-7.0*K2*Z*Z2**2+7.0*K2**2*Z*Z2-Z*K2**3)/K2pZ2**8
     504    elif k == 1:
     505        if N == 2:
     506            return 2.0*K/K2pZ2**2
     507        elif N == 3:
     508            return 8.0*K*Z/K2pZ2**3
     509        elif N == 4:
     510            return 8.0*K*(5.0*Z2-K2)/K2pZ2**4
     511        elif N == 5:
     512            return 48.0*K*Z*(5.0*Z2-3.0*K2)/K2pZ2**5
     513        elif N == 6:
     514            return 48.0*K*(35.0*Z2**2-42.0*K2*Z2+3.0*K2**2)/K2pZ2**6
     515        elif N == 7:
     516            return 1920.0*K*Z*(7.0*Z2**2-14.0*K2*Z2+3.0*K2**2)/K2pZ2**7
     517        elif N == 8:
     518            return 5760.0*K*(21.0*Z2**3-63.0*K2*Z2**2+27.0*K2**2*Z2-K2**3)/K2pZ2**8
     519    elif k == 2:
     520        if N == 3:
     521            return 8.0*K2/K2pZ2**3
     522        elif N == 4:
     523            return 48.0*K2*Z/K2pZ2**4
     524        elif N == 5:
     525            return 48.0*K2*(7.0*Z2-K2)/K2pZ2**5
     526        elif N == 6:
     527            return 384.0*K2*Z*(7.0*Z2-3.0*K2)/K2pZ2**6
     528        elif N == 7:
     529            return 1152.0*K2*(21.0*Z2**2-18.0*K2*Z2+K2**2)/K2pZ2**7
     530        elif N == 8:
     531            return 11520.0*K2*Z*(21.0*Z2**2-30.0*K2*Z2+5.0*K2**2)/K2pZ2**8
     532    elif k == 3:
     533        if N == 4:
     534            return 48.0*K**3/K2pZ2**4
     535        elif N == 5:
     536            return 384.0*K**3*Z/K2pZ2**5
     537        elif N == 6:
     538            return 384.0*K**3*(9.0*Z2-K2)/K2pZ2**6
     539        elif N == 7:
     540            return 11520.0*K**3*Z*(3.0*Z2-K2)/K2pZ2**7
     541        elif N == 8:
     542            return 11520.0*K**3*(33.0*Z2**2-22.0*K2*Z2+K2**2)/K2pZ2**8
     543    elif k == 4:
     544        if N == 5:
     545            return 384.0*K2**2/K2pZ2**5
     546        elif N == 6:
     547            return 3840.0*K2**2*Z/K2pZ2**6
     548        elif N == 7:
     549            return 3840.0*K2**2*(11.0*Z2-K2)/K2pZ2**7
     550        elif N == 8:
     551            return 46080.0*K**5*(13.0*Z2-K2)/K2pZ2**8
     552    elif k == 5:
     553        if N == 6:
     554            return 3840.0*K**5/K2pZ2**6
     555        elif N == 7:
     556            return 46080.0*Z*K**5/K2pZ2**7
     557        elif N == 8:
     558            return 46080.0*K**5*(13.0*Z2-K2)/K2pZ2**8
     559    elif k == 6:
     560        if N == 7:
     561            return 46080.0*K**6/K2pZ2**7
     562        elif N == 8:
     563            return 645120.0*Z*K2**3/K2pZ2**8
     564    elif k == 7:
     565        if N == 8:
     566            return 645120.0*K**7/K2pZ2**8
    445567       
     568       
     569           
     570   
     571   
     572   
    446573def BlenResCW(Els,BLtables,wave):
    447574    FP = np.zeros(len(Els))

  • TabularUnified trunk/GSASIIIO.py

    r5634 r5663  
    14731473        self.parmDict.update(rbDict)
    14741474        rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    1475         Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave =  \
     1475        Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,ORBtables,BLtables,MFtables,maxSSwave =  \
    14761476            G2stIO.GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
    14771477        self.parmDict.update(phaseDict)
    … …  
    15101510        rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    15111511        rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)  # done twice, needed?
    1512         Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave = \
     1512        Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,ORBtables,BLtables,MFtables,maxSSwave = \
    15131513            G2stIO.GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
    15141514        msg = G2mv.EvaluateMultipliers(constrDict,phaseDict)

  • TabularUnified trunk/GSASIIconstrGUI.py

    r5586 r5663  
    385385    rbVary, rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict, Print=False)
    386386    parmDict.update(rbDict)
    387     (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave) = \
     387    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,ORBtables,BLtables,MFtables,maxSSwave) = \
    388388        G2stIO.GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
    389389    parmDict.update(phaseDict)
    … …  
    446446        elif 'AU' in name:
    447447            namelist = ['AUiso','AU11','AU22','AU33','AU12','AU13','AU23']
     448        elif 'Akappa' in name:
     449            namelist = ['Akappa%d'%i for i in range(6)]
     450        elif 'ANe' in name:
     451            namelist = ['ANe0','ANe1']
     452        elif 'AD(1' in name:
     453            orb = name.split(':')[2][-1]
     454            namelist = ['AD(1,0)'+orb,'AD(1,1)'+orb,'AD(1,-1)'+orb]
     455        elif 'AD(2' in name:
     456            orb = name.split(':')[2][-1]
     457            namelist = ['AD(2,0)'+orb,'AD(2,1)'+orb,'AD(2,-1)'+orb,'AD(2,2)'+orb,'AD(2,-2)'+orb]
     458        elif 'AD(4' in name:
     459            orb = name.split(':')[2][-1]
     460            namelist = ['AD(4,0)'+orb,'AD(4,1)'+orb,'AD(4,-1)'+orb,'AD(4,2)'+orb,'AD(4,-2)'+orb,
     461                'AD(4,3)'+orb,'AD(4,-3)'+orb,'AD(4,4)'+orb,'AD(4,-4)'+orb]
    448462        elif 'AM' in name:
    449463            namelist = ['AMx','AMy','AMz']
    … …  
    15821596    # create a list of the phase variables
    15831597    symHolds = []
    1584     (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,EFtable,BLtable,MFtable,maxSSwave) = \
     1598    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,EFtable,ORBtables,BLtable,MFtable,maxSSwave) = \
    15851599        G2stIO.GetPhaseData(Phases,rbIds=rbIds,Print=False,symHold=symHolds)
    15861600    phaseList = []
    … …  
    16811695    if errmsg:
    16821696        G2frame.ErrorDialog('Constraint Error',
    1683                             'Error in constraints.\nCheck console output for more information'+
    1684                             ' or press "Show Errors" & "Show Warnings" buttons',
    1685                             parent=G2frame)
     1697            'Error in constraints.\nCheck console output for more information'+
     1698            ' or press "Show Errors" & "Show Warnings" buttons',
     1699            parent=G2frame)
    16861700        if seqhistnum is None:
    16871701            print ('\nError message(s):\n',errmsg)

  • TabularUnified trunk/GSASIIctrlGUI.py

    r5656 r5663  
    26552655    dlg.Destroy()
    26562656   
    2657 def G2ScrolledGrid(G2frame,lbl,title,tbl,colLbls,colTypes,maxSize=(600,300)):
     2657def G2ScrolledGrid(G2frame,lbl,title,tbl,colLbls,colTypes,maxSize=(600,300),comment=''):
    26582658    '''Display a scrolled table of information in a dialog window
    26592659
    … …  
    26672667    :param list maxSize: Maximum size for the table in points. Defaults to
    26682668      (600,300)
     2669      :param str comment: optional text that appears below table
    26692670
    26702671    Example::
    … …  
    26762677
    26772678    '''
     2679       
    26782680    dlg = wx.Dialog(G2frame,title=title,style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
    26792681    sizer = wx.BoxSizer(wx.VERTICAL)
    … …  
    27012703    scrGrid.Scroll(0,0)
    27022704    sizer.Add(scrGrid,1,wx.EXPAND,1)
     2705    if len(comment):
     2706        sizer.Add(wx.StaticText(dlg,label=comment))
    27032707       
    27042708    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    41604164                radiiSizer.Add(aRadii,0,WACV)
    41614165        mainSizer.Add(radiiSizer,0,wx.EXPAND)
     4166        Names = ['Bond']
     4167        factorSizer = wx.FlexGridSizer(0,2,5,5)
    41624168        if self.Angle:
    4163             factorSizer = wx.FlexGridSizer(0,2,5,5)
    41644169            Names = ['Bond','Angle']
    4165             for i,name in enumerate(Names):
    4166                 factorSizer.Add(wx.StaticText(self.panel,-1,name+' search factor'),0,WACV)
    4167                 bondFact = ValidatedTxtCtrl(self.panel,data['Factors'],i,nDig=(10,3))
    4168                 factorSizer.Add(bondFact)
    4169             mainSizer.Add(factorSizer,0,wx.EXPAND)
     4170        for i,name in enumerate(Names):
     4171            factorSizer.Add(wx.StaticText(self.panel,-1,name+' search factor'),0,WACV)
     4172            bondFact = ValidatedTxtCtrl(self.panel,data['Factors'],i,nDig=(10,3))
     4173            factorSizer.Add(bondFact)
     4174        mainSizer.Add(factorSizer,0,wx.EXPAND)
    41704175       
    41714176        OkBtn = wx.Button(self.panel,-1,"Ok")
    … …  
    66926697        if style is None: style = wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER
    66936698
    6694         super(self.__class__,self).__init__(parent, wx.ID_ANY, title,
    6695                                             style=style, size=(-1,-1))
     6699        super(self.__class__,self).__init__(parent, wx.ID_ANY, title,style=style, size=(-1,-1))
    66966700        self.Bind(wx.EVT_CLOSE, self._onClose)
    66976701        mainSizer = wx.BoxSizer(wx.VERTICAL)
    … …  
    67886792        lblSizer.Add(self.Labels)
    67896793        self.RwPanel = wx.lib.scrolledpanel.ScrolledPanel(self, wx.ID_ANY, size=(180, 130),
    6790                                  style = wx.SUNKEN_BORDER)
     6794            style = wx.SUNKEN_BORDER)
    67916795        lblSizer.Add(self.RwPanel,1,wx.EXPAND,1)
    67926796
    67936797        self.Labels.label = {}
    67946798        for i in range(self.rows):
    6795             self.Labels.label[i] = wx.StaticText(self,wx.ID_ANY,'',
    6796                                                 style=wx.ALIGN_CENTER_VERTICAL)
     6799            self.Labels.label[i] = wx.StaticText(self,wx.ID_ANY,'',style=wx.ALIGN_CENTER_VERTICAL)
    67976800            self.Labels.Add(self.Labels.label[i],(i,0))
    67986801        mainsizer = wx.BoxSizer(wx.VERTICAL)
    … …  
    68286831        lbls = []
    68296832        for i in range(self.rows-1):
    6830             lbls.append(wx.StaticText(self.RwPanel,wx.ID_ANY,'',
    6831                                           style=wx.ALIGN_CENTER))
     6833            lbls.append(wx.StaticText(self.RwPanel,wx.ID_ANY,'',style=wx.ALIGN_CENTER))
    68326834            self.gridSiz.Add(lbls[-1],(1+i,col))
    68336835        return col,lbls
    … …  
    68946896            self.seqHist = nextHist
    68956897            self.SeqCount += 1
    6896             self.gauge.SetValue(int(min(self.gaugemaximum,
    6897                                 100.*self.SeqCount/self.SeqLen)))
     6898            self.gauge.SetValue(int(min(self.gaugemaximum,100.*self.SeqCount/self.SeqLen)))
    68986899
    68996900    def _plotBar(self,h):
    … …  
    69106911            self.histOff[h] = len([i for i in self.plotted if i >= 0])/wid
    69116912        self.plotaxis.bar(np.array(range(l))+self.histOff[h],self.fitVals[h],
    6912                               width=1/wid,label=lbl,color=sym)
     6913            width=1/wid,label=lbl,color=sym)
    69136914        self.plotted.append(h)
    69146915       

  • TabularUnified trunk/GSASIIdataGUI.py

    r5661 r5663  
    823823        item.Enable(state) # disabled during sequential fits
    824824        self.Bind(wx.EVT_MENU, self.OnRefinePartials, id=item.GetId())
     825        item = parent.Append(wx.ID_ANY,'&Parameter Impact\tCTRL+I','Perform a derivative calculation')
     826        self.Bind(wx.EVT_MENU, self.OnDerivCalc, id=item.GetId())
    825827       
    826828        item = parent.Append(wx.ID_ANY,'Save partials as csv','Save the computed partials as a csv file')
    … …  
    838840        item = parent.Append(wx.ID_ANY,'&Run PlotXNFF','Plot X-ray, neutron & magnetic form factors')
    839841        self.Bind(wx.EVT_MENU, self.OnRunPlotXNFF, id=item.GetId())
    840         item = parent.Append(wx.ID_ANY,'&Parameter Impact\tCTRL+I','Perform a derivative calculation')
    841         self.Bind(wx.EVT_MENU, self.OnDerivCalc, id=item.GetId())
    842842
    843843#        if GSASIIpath.GetConfigValue('debug'): # allow exceptions for debugging
    … …  
    52705270        return phaseNames
    52715271   
     5272    def GetHistogramTypes(self):
     5273        """ Returns a list of histogram types found in the GSASII data tree
     5274       
     5275        :return: list of histogram types
     5276       
     5277        """
     5278        HistogramTypes = []
     5279        if self.GPXtree.GetCount():
     5280            item, cookie = self.GPXtree.GetFirstChild(self.root)
     5281            while item:
     5282                name = self.GPXtree.GetItemText(item)
     5283                if name[:4] in ['PWDR','HKLF']:
     5284                    Inst = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,item,'Instrument Parameters'))
     5285                    HistogramTypes.append(Inst[0]['Type'][0])
     5286                item, cookie = self.GPXtree.GetNextChild(self.root, cookie)
     5287        return HistogramTypes
     5288   
    52725289    def GetHistogramNames(self,hType):
    52735290        """ Returns a list of histogram names found in the GSASII data tree
    … …  
    53995416        rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
    54005417        rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[],'Spin':[]})
    5401         Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,EFtable,BLtable,MFtable,maxSSwave = \
     5418        Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,EFtable,ORBtables,BLtable,MFtable,maxSSwave = \
    54025419            G2stIO.GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)       
    54035420        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,histDict,Print=False,resetRefList=False)
    … …  
    54965513        Controls = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,self.root, 'Controls'))
    54975514        wx.BeginBusyCursor()
    5498         pdlg = wx.ProgressDialog('Computing derivatives',
    5499                     'Impact computation in progress',100,
    5500                     parent=self,
    5501 #                    style = wx.PD_ELAPSED_TIME)
    5502                     style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
     5515        pdlg = wx.ProgressDialog('Computing derivatives','Impact computation in progress',100,
     5516            parent=self,style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
    55035517        pdlg.CenterOnParent()
    55045518        derivCalcs,varyList = G2stMn.Refine(self.GSASprojectfile,pdlg,allDerivs=True)
    … …  
    55185532            txt = txt.replace('Ph=','Phase: ')
    55195533            txt = txt.replace('Pwd=','Histogram: ')
    5520             #print(f'{x}\t{derivCalcs[x][1]}\t{txt}')
    55215534            tbl.append([x,x in varyList,derivCalcs[x][1],txt])
    5522         G2G.G2ScrolledGrid(self,'Parameter Impact Results',
    5523                            'Impact Results',tbl,colLbls,
    5524                            colTypes,maxSize=(700,400))
     5535        G2G.G2ScrolledGrid(self,'Parameter Impact Results','Impact Results',tbl,colLbls,colTypes,
     5536            maxSize=(700,400),comment=' Cite: B.H. Toby, IUCrJ, to be published')
    55255537
    55265538    def OnRefine(self,event):
    … …  
    71027114        self.PostfillDataMenu()
    71037115       
     7116        # Phase / Deformation tab
     7117        G2G.Define_wxId('wxID_DEFORMSETSEL','wxID_DEFORMDISTSET')
     7118        self.DeformationMenu = wx.MenuBar()
     7119        self.PrefillDataMenu(self.DeformationMenu)
     7120        self.DeformationMenu.Append(menu=wx.Menu(title=''),title='Select tab')
     7121        self.DeformationEdit = wx.Menu(title='')
     7122        self.DeformationMenu.Append(menu=self.DeformationEdit, title='Edit Deformations')
     7123        self.DeformationEdit.Append(G2G.wxID_DEFORMSETSEL,'Select atom','Select atom for deformation study')
     7124        self.DeformationEdit.Append(G2G.wxID_DEFORMDISTSET,'Set bond parms','Set bond selection parameters')
     7125        self.PostfillDataMenu()
     7126       
    71047127        # Phase / Draw Atoms tab
    71057128        G2G.Define_wxId('wxID_DRAWATOMSTYLE', 'wxID_DRAWATOMLABEL', 'wxID_DRAWATOMCOLOR', 'wxID_DRAWATOMRESETCOLOR',
    … …  
    72437266    # end of GSAS-II menu definitions
    72447267   
    7245 #####  Notebook Tree Item editor ##############################################
     7268####  Notebook Tree Item editor ##############################################
    72467269def UpdateNotebook(G2frame,data):
    72477270    '''Called when the data tree notebook entry is selected. Allows for
    … …  
    72767299    G2frame.SendSizeEvent()
    72777300
    7278 #####  Comments ###############################################################
     7301####  Comments ###############################################################
    72797302def UpdateComments(G2frame,data):
    72807303    '''Place comments into the data window
    … …  
    72947317
    72957318           
    7296 #####  Controls Tree Item editor ##############################################
     7319####  Controls Tree Item editor ##############################################
    72977320def UpdateControls(G2frame,data):
    72987321    '''Edit overall GSAS-II controls in main Controls data tree entry
    … …  
    76007623    G2frame.SendSizeEvent()
    76017624   
    7602 #####  Main PWDR panel ########################################################   
     7625####  Main PWDR panel ########################################################   
    76037626def UpdatePWHKPlot(G2frame,kind,item):
    76047627    '''Called when the histogram main tree entry is called. Displays the

  • TabularUnified trunk/GSASIIexprGUI.py

    r5574 r5663  
    920920        bNames = []
    921921        if self.Oatom:
    922             neigh = G2mth.FindAllNeighbors(Phase,self.Oatom,aNames,
    923                                            searchType='Angle')[0]
     922            neigh = G2mth.FindAllNeighbors(Phase,self.Oatom,aNames,searchType='Angle')[0]
    924923        if neigh:
    925924            for iA,aName in enumerate(neigh):

  • TabularUnified trunk/GSASIImath.py

    r5639 r5663  
    932932    return drawAtoms
    933933   
     934def FindOctahedron(results):
     935    Octahedron = np.array([[1.,0,0],[0,1.,0],[0,0,1.],[-1.,0,0],[0,-1.,0],[0,0,-1.]])
     936    Polygon = np.array([result[3] for result in results])
     937    Dists = np.array([np.sqrt(np.sum(axis**2)) for axis in Polygon])
     938    bond = np.mean(Dists)
     939    std = np.std(Dists)
     940    Norms = Polygon/Dists[:,nxs]
     941    Tilts = acosd(np.dot(Norms,Octahedron[0]))
     942    iAng = np.argmin(Tilts)
     943    Qavec = np.cross(Norms[iAng],Octahedron[0])
     944    QA = AVdeg2Q(Tilts[iAng],Qavec)
     945    newNorms = prodQVQ(QA,Norms)
     946    Rots = acosd(np.dot(newNorms,Octahedron[1]))
     947    jAng = np.argmin(Rots)
     948    Qbvec = np.cross(Norms[jAng],Octahedron[1])
     949    QB = AVdeg2Q(Rots[jAng],Qbvec)
     950    QQ = prodQQ(QA,QB)   
     951    newNorms = prodQVQ(QQ,Norms)
     952    dispVecs = np.array([norm[:,nxs]-Octahedron.T for norm in newNorms])
     953    disp = np.sqrt(np.sum(dispVecs**2,axis=1))
     954    dispids = np.argmin(disp,axis=1)
     955    vecDisp = np.array([dispVecs[i,:,dispid] for i,dispid in enumerate(dispids)])
     956    Disps = np.array([disp[i,dispid] for i,dispid in enumerate(dispids)])
     957    meanDisp = np.mean(Disps)
     958    stdDisp = np.std(Disps)
     959    A,V = Q2AVdeg(QQ)
     960    return bond,std,meanDisp,stdDisp,A,V,vecDisp
     961   
     962def FindTetrahedron(results):
     963    Tetrahedron = np.array([[1.,1,1],[1,-1,-1],[-1,1,-1],[-1,-1,1]])/np.sqrt(3.)
     964    Polygon = np.array([result[3] for result in results])
     965    Dists = np.array([np.sqrt(np.sum(axis**2)) for axis in Polygon])
     966    bond = np.mean(Dists)
     967    std = np.std(Dists)
     968    Norms = Polygon/Dists[:,nxs]
     969    Tilts = acosd(np.dot(Norms,Tetrahedron[0]))
     970    iAng = np.argmin(Tilts)
     971    Qavec = np.cross(Norms[iAng],Tetrahedron[0])
     972    QA = AVdeg2Q(Tilts[iAng],Qavec)
     973    newNorms = prodQVQ(QA,Norms)
     974    Rots = acosd(np.dot(newNorms,Tetrahedron[1]))
     975    jAng = np.argmin(Rots)
     976    Qbvec = np.cross(Norms[jAng],Tetrahedron[1])
     977    QB = AVdeg2Q(Rots[jAng],Qbvec)
     978    QQ = prodQQ(QA,QB)   
     979    newNorms = prodQVQ(QQ,Norms)
     980    dispVecs = np.array([norm[:,nxs]-Tetrahedron.T for norm in newNorms])
     981    disp = np.sqrt(np.sum(dispVecs**2,axis=1))
     982    dispids = np.argmin(disp,axis=1)
     983    vecDisp = np.array([dispVecs[i,:,dispid] for i,dispid in enumerate(dispids)])
     984    Disps = np.array([disp[i,dispid] for i,dispid in enumerate(dispids)])
     985    meanDisp = np.mean(Disps)
     986    stdDisp = np.std(Disps)
     987    A,V = Q2AVdeg(QQ)
     988    return bond,std,meanDisp,stdDisp,A,V,vecDisp
     989   
    934990def FindNeighbors(phase,FrstName,AtNames,notName=''):
    935991    General = phase['General']
    … …  
    9641020    return Neigh,[OId,Ids]
    9651021
    966 def FindOctahedron(results):
    967     Octahedron = np.array([[1.,0,0],[0,1.,0],[0,0,1.],[-1.,0,0],[0,-1.,0],[0,0,-1.]])
    968     Polygon = np.array([result[3] for result in results])
    969     Dists = np.array([np.sqrt(np.sum(axis**2)) for axis in Polygon])
    970     bond = np.mean(Dists)
    971     std = np.std(Dists)
    972     Norms = Polygon/Dists[:,nxs]
    973     Tilts = acosd(np.dot(Norms,Octahedron[0]))
    974     iAng = np.argmin(Tilts)
    975     Qavec = np.cross(Norms[iAng],Octahedron[0])
    976     QA = AVdeg2Q(Tilts[iAng],Qavec)
    977     newNorms = prodQVQ(QA,Norms)
    978     Rots = acosd(np.dot(newNorms,Octahedron[1]))
    979     jAng = np.argmin(Rots)
    980     Qbvec = np.cross(Norms[jAng],Octahedron[1])
    981     QB = AVdeg2Q(Rots[jAng],Qbvec)
    982     QQ = prodQQ(QA,QB)   
    983     newNorms = prodQVQ(QQ,Norms)
    984     dispVecs = np.array([norm[:,nxs]-Octahedron.T for norm in newNorms])
    985     disp = np.sqrt(np.sum(dispVecs**2,axis=1))
    986     dispids = np.argmin(disp,axis=1)
    987     vecDisp = np.array([dispVecs[i,:,dispid] for i,dispid in enumerate(dispids)])
    988     Disps = np.array([disp[i,dispid] for i,dispid in enumerate(dispids)])
    989     meanDisp = np.mean(Disps)
    990     stdDisp = np.std(Disps)
    991     A,V = Q2AVdeg(QQ)
    992     return bond,std,meanDisp,stdDisp,A,V,vecDisp
    993    
    994 def FindTetrahedron(results):
    995     Tetrahedron = np.array([[1.,1,1],[1,-1,-1],[-1,1,-1],[-1,-1,1]])/np.sqrt(3.)
    996     Polygon = np.array([result[3] for result in results])
    997     Dists = np.array([np.sqrt(np.sum(axis**2)) for axis in Polygon])
    998     bond = np.mean(Dists)
    999     std = np.std(Dists)
    1000     Norms = Polygon/Dists[:,nxs]
    1001     Tilts = acosd(np.dot(Norms,Tetrahedron[0]))
    1002     iAng = np.argmin(Tilts)
    1003     Qavec = np.cross(Norms[iAng],Tetrahedron[0])
    1004     QA = AVdeg2Q(Tilts[iAng],Qavec)
    1005     newNorms = prodQVQ(QA,Norms)
    1006     Rots = acosd(np.dot(newNorms,Tetrahedron[1]))
    1007     jAng = np.argmin(Rots)
    1008     Qbvec = np.cross(Norms[jAng],Tetrahedron[1])
    1009     QB = AVdeg2Q(Rots[jAng],Qbvec)
    1010     QQ = prodQQ(QA,QB)   
    1011     newNorms = prodQVQ(QQ,Norms)
    1012     dispVecs = np.array([norm[:,nxs]-Tetrahedron.T for norm in newNorms])
    1013     disp = np.sqrt(np.sum(dispVecs**2,axis=1))
    1014     dispids = np.argmin(disp,axis=1)
    1015     vecDisp = np.array([dispVecs[i,:,dispid] for i,dispid in enumerate(dispids)])
    1016     Disps = np.array([disp[i,dispid] for i,dispid in enumerate(dispids)])
    1017     meanDisp = np.mean(Disps)
    1018     stdDisp = np.std(Disps)
    1019     A,V = Q2AVdeg(QQ)
    1020     return bond,std,meanDisp,stdDisp,A,V,vecDisp
    1021    
    1022 def FindAllNeighbors(phase,FrstName,AtNames,notName='',Orig=None,Short=False,
    1023                      searchType='Bond'):
     1022def FindAllNeighbors(phase,FrstName,AtNames,notName='',Orig=None,Short=False,searchType='Bond'):
    10241023    '''Find neighboring atoms
    10251024    Uses Bond search criteria unless searchType is set to non-default
    … …  
    10591058    Ids = []
    10601059    sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
    1061     sumR = np.reshape(np.tile(sumR,27),(27,-1))
     1060    sumR = np.reshape(np.tile(sumR,27),(27,-1))     #27 = 3x3x3 unit cell block
    10621061    results = []
    10631062    for xyz in XYZ:
    … …  
    10651064    for iA,result in enumerate(results):
    10661065        for [Txyz,Top,Tunit,Spn] in result:
    1067             Dx = np.array([Txyz-Oxyz+unit for unit in Units])
     1066            Dx = (Txyz-np.array(Oxyz))+Units
    10681067            dx = np.inner(Dx,Amat)
    10691068            dist = np.sqrt(np.sum(dx**2,axis=1))

  • TabularUnified trunk/GSASIIobj.py

    r5578 r5663  
    625625        'Amul': 'Atomic site multiplicity value',
    626626        'AM([xyz])$' : 'Atomic magnetic moment parameter, \\1',
     627        # Atom deformation parameters
     628        'Akappa([0-6])'  : ' Atomic orbital softness for orbital, \\1',
     629        'ANe([01])' : ' Atomic <j0> orbital population for orbital, \\1',
     630        'AD\\([0-6],[0-6]\\)([0-6])' : ' Atomic sp. harm. coeff for orbital, \\1',
     631        'AD\\([0-6],-[0-6]\\)([0-6])' : ' Atomic sp. harm. coeff for orbital, \\1',     #need both!
    627632        # Hist (:h:<var>) & Phase (HAP) vars (p:h:<var>)
    628633        'Back(.*)': 'Background term #\\1',
    … …  
    786791        if m:
    787792            reVarDesc[key]
    788             return m.expand(reVarDesc[key])
     793            try:
     794                return m.expand(reVarDesc[key])
     795            except:
     796                print('Error in key: %s'%key)
    789797    return None
    790798

  • TabularUnified trunk/GSASIIphsGUI.py

    r5649 r5663  
    5959import numpy as np
    6060import numpy.linalg as nl
     61import numpy.ma as ma
    6162import atmdata
    6263import ISODISTORT as ISO
    … …  
    1004610047            choices.append(val)
    1004710048        if not choices: return
    10048         dlg = G2G.G2MultiChoiceDialog(G2frame,
    10049                     'Select atoms','Choose atoms to select',choices)
     10049        dlg = G2G.G2MultiChoiceDialog(G2frame,'Select atoms','Choose atoms to select',choices)
    1005010050        indx = []
    1005110051        if dlg.ShowModal() == wx.ID_OK:
    … …  
    1070610706
    1070710707        SetPhaseWindow(drawOptions,mainSizer)
     10708       
     10709####  Deformation form factor routines ################################################################
     10710
     10711    def SetDefDist(event):
     10712        generalData = data['General']
     10713        DisAglCtls = {}
     10714        if 'DisAglCtls' in generalData:
     10715            DisAglCtls = generalData['DisAglCtls']
     10716        dlg = G2G.DisAglDialog(G2frame,DisAglCtls,generalData,Angle=False)
     10717        if dlg.ShowModal() == wx.ID_OK:
     10718            generalData['DisAglCtls'] = dlg.GetData()
     10719        UpdateDeformation()
     10720        event.StopPropagation()
     10721       
     10722    def SelDeformAtom(event):
     10723        'select deformation atom using a filtered listbox'
     10724        generalData = data['General']
     10725        cx,ct,cs,cia = generalData['AtomPtrs']
     10726        choices = []
     10727        types = []
     10728        Ids = []
     10729        for atom in data['Atoms']:
     10730            if atom[ct] in atmdata.OrbFF:
     10731                choices.append(atom[ct-1])
     10732                types.append(atom[ct])
     10733                Ids.append(atom[cia+8])
     10734        if not choices: return      #no atoms in phase!
     10735        dlg = G2G.G2SingleChoiceDialog(G2frame,'Select atom','Choose atom to select',choices)
     10736        indx = -1
     10737        if dlg.ShowModal() == wx.ID_OK:
     10738            indx = dlg.GetSelection()
     10739            orbs = atmdata.OrbFF[types[indx]]
     10740            data['Deformations'][Ids[indx]] = []
     10741            for orb in orbs:
     10742                if 'core' in orb:
     10743                    continue        #skip core - has no parameters
     10744                else:
     10745                    if 'j0' in orb:
     10746                        data['Deformations'][Ids[indx]].append([orb,{'Ne':[float(orbs[orb]['Ne']),False],'kappa':[1.0,False]}])   #no sp. harm for j0 terms
     10747                    elif 'j' in orb:
     10748                        orbDict = {'kappa':[1.0,False],}
     10749                        Order = int(orb.split('>')[0][-1])
     10750                        cofNames,cofSgns = G2lat.GenRBCoeff('1','1',Order)
     10751                        cofNames = [name.replace('C','D') for name in cofNames]
     10752                        cofTerms = {name:[0.0,False] for name in cofNames if str(Order) in name}
     10753                        for name in cofNames:
     10754                            if str(Order) in name and '0' not in name:
     10755                                negname = name.replace(',',',-')
     10756                                cofTerms.update({negname:[0.0,False]})
     10757                        orbDict.update(cofTerms)
     10758                        data['Deformations'][Ids[indx]].append([orb,orbDict])
     10759        dlg.Destroy()
     10760        if indx < 0: return
     10761        drawAtoms.ClearSelection()       
     10762        drawAtoms.SelectRow(indx,True)
     10763        G2plt.PlotStructure(G2frame,data)
     10764        UpdateDeformation()
     10765        event.StopPropagation()
     10766
     10767    def UpdateDeformation():
     10768       
     10769        def OnDeformRef(event):
     10770            Obj = event.GetEventObject()
     10771            dId,oId,dkey = Indx[Obj.GetId()]
     10772            deformationData[dId][oId][1][dkey][1] = not deformationData[dId][oId][1][dkey][1]
     10773
     10774        def OnDelAtm(event):
     10775            Obj = event.GetEventObject()
     10776            dId = Indx[Obj.GetId()]
     10777            del deformationData[dId]
     10778            UpdateDeformation()
     10779       
     10780        # UpdateDeformation exectable code starts here
     10781        generalData = data['General']
     10782        cx,ct,cs,cia = generalData['AtomPtrs']
     10783        Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     10784        atomData = data['Atoms']
     10785        AtLookUp = G2mth.FillAtomLookUp(atomData,cia+8)
     10786        AtNames = [atom[ct-1] for atom in atomData]
     10787        deformationData = data['Deformations']
     10788       
     10789        if deformation.GetSizer():
     10790            deformation.GetSizer().Clear(True)
     10791        mainSizer = wx.BoxSizer(wx.VERTICAL)
     10792        topSizer = wx.BoxSizer(wx.HORIZONTAL)
     10793        topSizer.Add(wx.StaticText(deformation,label=' Atomic deformation data:'),0,WACV)
     10794        # add help button to bring up help web page - at right side of window
     10795        topSizer.Add((-1,-1),1,wx.EXPAND)
     10796        topSizer.Add(G2G.HelpButton(deformation,helpIndex=G2frame.dataWindow.helpKey))
     10797        mainSizer.Add(topSizer,0,wx.EXPAND)
     10798        Indx = {}
     10799       
     10800        for dId in deformationData:
     10801            atom = atomData[AtLookUp[dId]]
     10802            neigh = G2mth.FindAllNeighbors(data,atom[ct-1],AtNames)
     10803            lineSizer = wx.BoxSizer(wx.HORIZONTAL)
     10804            lineSizer.Add(wx.StaticText(deformation,label=' For atom %s, site sym %s:'%(atom[ct-1],atom[cs])),0,WACV)
     10805            names = []
     10806            if not len(neigh[0]):
     10807                lineSizer.Add(wx.StaticText(deformation,label=' No neighbors found; Do Set bond parms to expand search'),0,WACV)
     10808            else:
     10809                names = [item[0] for item in neigh[0]]
     10810                lineSizer.Add(wx.StaticText(deformation,label=' Neighbors: '+str(names)),0,WACV)
     10811            delAtm = wx.Button(deformation,label='Delete')
     10812            delAtm.Bind(wx.EVT_BUTTON,OnDelAtm)
     10813            Indx[delAtm.GetId()] = dId
     10814            lineSizer.Add(delAtm,0,WACV)
     10815            mainSizer.Add(lineSizer)
     10816            orbSizer = wx.FlexGridSizer(0,9,2,2)
     10817            for iorb,orb in enumerate(deformationData[dId]):
     10818                orbSizer.Add(wx.StaticText(deformation,label=orb[0]+' kappa:'))
     10819                orbSizer.Add(G2G.ValidatedTxtCtrl(deformation,orb[1]['kappa'],0,nDig=(8,3),xmin=0.,xmax=10.))
     10820                Tcheck = wx.CheckBox(deformation,-1,'Refine?')
     10821                Tcheck.SetValue(orb[1]['kappa'][1])
     10822                Tcheck.Bind(wx.EVT_CHECKBOX,OnDeformRef)
     10823                Indx[Tcheck.GetId()] = [dId,iorb,'kappa']
     10824                orbSizer.Add(Tcheck)
     10825                if '<j0>' in orb[0]:
     10826                    orbSizer.Add(wx.StaticText(deformation,label=' Ne:'))
     10827                    orbSizer.Add(G2G.ValidatedTxtCtrl(deformation,orb[1]['Ne'],0,nDig=(8,3),xmin=0.,xmax=10.))
     10828                    Tcheck = wx.CheckBox(deformation,-1,'Refine?')
     10829                    Tcheck.SetValue(orb[1]['Ne'][1])
     10830                    Tcheck.Bind(wx.EVT_CHECKBOX,OnDeformRef)
     10831                    Indx[Tcheck.GetId()] = [dId,iorb,'Ne']
     10832                    orbSizer.Add(Tcheck)
     10833                    for i in range(3): orbSizer.Add((5,5),0)
     10834                    continue
     10835                nItem = 0
     10836                for item in orb[1]:
     10837                    if 'D' in item:
     10838                        nItem += 1
     10839                        orbSizer.Add(wx.StaticText(deformation,label=item+':'))
     10840                        orbSizer.Add(G2G.ValidatedTxtCtrl(deformation,orb[1][item],0,nDig=(8,3)))
     10841                        Tcheck = wx.CheckBox(deformation,-1,'Refine?')
     10842                        Tcheck.SetValue(orb[1][item][1])
     10843                        Tcheck.Bind(wx.EVT_CHECKBOX,OnDeformRef)
     10844                        Indx[Tcheck.GetId()] = [dId,iorb,item]
     10845                        orbSizer.Add(Tcheck)
     10846                        if nItem in [2,4,6,8,10]:
     10847                            for i in range(3): orbSizer.Add((5,5),0)
     10848                for i in range(3): orbSizer.Add((5,5),0)
     10849                   
     10850            mainSizer.Add(orbSizer)   
     10851            G2G.HorizontalLine(mainSizer,deformation)
     10852
     10853        SetPhaseWindow(deformation,mainSizer)
    1070810854
    1070910855####  Texture routines ################################################################################       
    … …  
    1502615172            G2plt.PlotStructure(G2frame,data,firstCall=True)
    1502715173            UpdateDrawAtoms()
     15174        elif text == 'Deformation':
     15175            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.DeformationMenu)
     15176            G2plt.PlotStructure(G2frame,data,firstCall=True)
     15177            UpdateDeformation()
    1502815178        elif text == 'RB Models':
    1502915179            G2gd.SetDataMenuBar(G2frame,G2frame.dataWindow.RigidBodiesMenu)
    … …  
    1515715307        G2frame.Bind(wx.EVT_MENU, DrawLoadSel, id=G2G.wxID_DRAWLOADSEL)
    1515815308       
     15309        # Deformation form factors
     15310        FillSelectPageMenu(TabSelectionIdDict, G2frame.dataWindow.DeformationMenu)
     15311        G2frame.Bind(wx.EVT_MENU, SelDeformAtom, id=G2G.wxID_DEFORMSETSEL)
     15312        G2frame.Bind(wx.EVT_MENU, SetDefDist, id=G2G.wxID_DEFORMDISTSET)
     15313       
    1515915314        # RB Models
    1516015315        FillSelectPageMenu(TabSelectionIdDict, G2frame.dataWindow.RigidBodiesMenu)
    … …  
    1528515440    if 'ISODISTORT' not in data:
    1528615441        data['ISODISTORT'] = {}
     15442    if 'Deformations' not in data:
     15443        data['Deformations'] = {}
    1528715444#end patch   
    1528815445
    … …  
    1532515482    Pages.append('Draw Atoms')
    1532615483   
     15484    if any('X' in item for item in G2frame.GetHistogramTypes()):
     15485        deformation = wx.ScrolledWindow(G2frame.phaseDisplay)
     15486        G2frame.phaseDisplay.AddPage(deformation,'Deformation')
     15487       
     15488        Pages.append('Deformation')
     15489   
    1532715490    if data['General']['Type'] not in ['faulted',] and not data['General']['Modulated']:
    1532815491        RigidBodies = wx.ScrolledWindow(G2frame.phaseDisplay)

  • TabularUnified trunk/GSASIIrestrGUI.py

    r5655 r5663  
    289289        try:
    290290            bondlst = G2mth.searchBondRestr(Lists['origin'],Lists['target'],
    291                                             bond,Factor,General['Type'],
    292                                             SGData,Amat,0.01,dlg)
     291                bond,Factor,General['Type'],SGData,Amat,0.01,dlg)
    293292            for newBond in bondlst:
    294293                if newBond not in bondRestData['Bonds']:
    … …  
    405404            VectB = []
    406405            for Tid,Ttype,Tcoord in targAtoms:
    407                 result = G2spc.GenAtom(Tcoord,SGData,False,Move=False)
     406                result = G2spc.GenAtom(Tcoord,SGData,All=False,Move=False)
    408407                BsumR = (Radii[Otype][0]+Radii[Ttype][0])*Factor
    409408                AsumR = (Radii[Otype][1]+Radii[Ttype][1])*Factor

  • TabularUnified trunk/GSASIIstrIO.py

    r5639 r5663  
    202202    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
    203203    parmDict.update(rbDict)
    204     (Natoms, atomIndx, phaseVary,phaseDict, pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave) = \
     204    (Natoms, atomIndx, phaseVary,phaseDict, pawleyLookup,FFtables,EFtables,ORBtables,BLtables,MFtables,maxSSwave) = \
    205205        GetPhaseData(Phases,RestraintDict=None,seqHistName=seqHist,rbIds=rbIds,Print=False) # generates atom symmetry constraints
    206206    parmDict.update(phaseDict)
    … …  
    11671167       in G2constrGUI and for parameter impact estimates in AllPrmDerivs).
    11681168    :returns: lots of stuff: Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
    1169         FFtables,EFtables,BLtables,MFtables,maxSSwave (see code for details).
     1169        FFtables,EFtables,ORBtables,BLtables,MFtables,maxSSwave (see code for details).
    11701170    '''
    11711171           
    … …  
    11931193                pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
    11941194                    (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fa[4],fb[0],fb[1],fb[2],fb[3],fb[4]))
     1195               
     1196    def PrintORBtable(ORBtable):
     1197        if not len(ORBtable):
     1198            return
     1199        pFile.write('\n X-ray orbital scattering factors:\n')
     1200        pFile.write('   Symbol          fa                                      fb                                      fc\n')
     1201        pFile.write(103*'-'+'\n')
     1202        for Ename in ORBtable:
     1203            orbdata = ORBtable[Ename]
     1204            for item in orbdata:
     1205                fa = orbdata[item]['fa']
     1206                fb = orbdata[item]['fb']
     1207                if 'core' in item:                   
     1208                    name = ('%s %s'%(Ename,item)).ljust(13)
     1209                else:
     1210                    name = '   '+item.ljust(10)
     1211                pFile.write(' %13s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
     1212                    (name,fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],orbdata[item]['fc']))
    11951213               
    11961214    def PrintMFtable(MFtable):
    … …  
    13771395                    '%10.5f'%(at[cmx+2])
    13781396                pFile.write(line+'\n')
    1379        
     1397               
     1398    def PrintDeformations(General,Atoms,Deformations):
     1399        cx,ct,cs,cia = General['AtomPtrs']
     1400        pFile.write('\n Atomic deformation parameters:\n')
     1401        for AtDef in Deformations:
     1402            atom = G2mth.GetAtomsById(Atoms,AtLookup,[AtDef,])[0]
     1403            DefTable = Deformations[AtDef]
     1404            pFile.write('\n Atom: %s at %10.5f, %10.5f, %10.5f sytsym: %s\n'%(atom[ct-1],atom[cx],atom[cx+1],atom[cx+2],atom[cs]))
     1405            pFile.write(135*'-'+'\n')
     1406            for orb in DefTable:
     1407                names = orb[0].ljust(12)+ 'names : '
     1408                values = 12*' '+'values: '
     1409                refine = 12*' '+'refine: '
     1410                for item in orb[1]:
     1411                    names += item.rjust(10)
     1412                    values += '%10.5f'%orb[1][item][0]
     1413                    refine += '%10s'%str(orb[1][item][1])
     1414                pFile.write(names+'\n')
     1415                pFile.write(values+'\n')
     1416                pFile.write(refine+'\n')
    13801417               
    13811418    def PrintWaves(General,Atoms):
    … …  
    15891626    EFtables = {}                   #scattering factors - electrons
    15901627    MFtables = {}                   #Mag. form factors
    1591     BLtables = {}                   # neutrons
     1628    BLtables = {}                  # neutrons
     1629    ORBtables = {}
    15921630    Natoms = {}
    15931631    maxSSwave = {}
    … …  
    16021640        pId = PhaseData[name]['pId']
    16031641        pfx = str(pId)+'::'
     1642        Deformations = PhaseData[name].get('Deformations',[])
    16041643        FFtable = G2el.GetFFtable(General['AtomTypes'])
    16051644        EFtable = G2el.GetEFFtable(General['AtomTypes'])
    … …  
    16201659            cx,ct,cs,cia = General['AtomPtrs']
    16211660            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
     1661        DefAtm = []
     1662        for iAt in Deformations:
     1663            DefAtm.append(G2mth.GetAtomItemsById(Atoms,AtLookup,iAt,ct)[0])
     1664        ORBtable = G2el.GetORBtable(set(DefAtm))
     1665        ORBtables.update(ORBtable)
    16221666        PawleyRef = PhaseData[name].get('Pawley ref',[])
    16231667        cell = General['Cell']
    … …  
    18041848                                            G2mv.StoreEquivalence(name,(eqv,))
    18051849                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
     1850                           
     1851            if len(Deformations) and not General.get('doPawley'):
     1852                for iAt in Deformations:
     1853                    AtId = AtLookup[iAt]
     1854                    DefAtm = Deformations[iAt]
     1855                    for iorb,orb in enumerate(DefAtm):
     1856                        for parm in orb[1]:
     1857                            name = pfx+'A%s%d:%d'%(parm,iorb,AtId)
     1858                            phaseDict[name] = orb[1][parm][0]
     1859                            if orb[1][parm][1]:
     1860                                phaseVary.append(name)
     1861                   
    18061862            textureData = General['SH Texture']
    18071863            if textureData['Order']:    # and seqHistName is not None:
    … …  
    18221878                PrintFFtable(FFtable)
    18231879                PrintEFtable(EFtable)
     1880                PrintORBtable(ORBtable)
    18241881                PrintBLtable(BLtable)
    18251882                if General['Type'] == 'magnetic':
    … …  
    18451902                PrintRBObjects(resRBData,vecRBData,spnRBData)
    18461903                PrintAtoms(General,Atoms)
     1904                if len(Deformations):
     1905                    PrintDeformations(General,Atoms,Deformations)
     1906                   
    18471907                if General['Type'] == 'magnetic':
    18481908                    PrintMoments(General,Atoms)
    … …  
    18981958            phaseVary += pawleyVary
    18991959               
    1900     return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave
     1960    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,ORBtables,BLtables,MFtables,maxSSwave
    19011961   
    19021962def cellFill(pfx,SGData,parmDict,sigDict):
    … …  
    23182378                pFile.write(name+'\n')
    23192379                pFile.write(valstr+'\n')
     2380                pFile.write(sigstr+'\n')
     2381               
     2382    def PrintDeformationsAndSig(General,Atoms,Deformations,deformSig):
     2383        pFile.write('\n Atom deformations:\n')
     2384        cx,ct,cs,cia = General['AtomPtrs']
     2385        for AtDef in Deformations:
     2386            pFile.write(135*'-'+'\n')
     2387            atom = G2mth.GetAtomsById(Atoms,AtLookup,[AtDef,])[0]
     2388            AtId = AtLookup[AtDef]
     2389            DefTable = Deformations[AtDef]
     2390            pFile.write('\n Atom: %s at %10.5f, %10.5f, %10.5f sytsym: %s\n'%(atom[ct-1],atom[cx],atom[cx+1],atom[cx+2],atom[cs]))
     2391            pFile.write(135*'-'+'\n')
     2392            for iorb,orb in enumerate(DefTable):
     2393                names = orb[0].ljust(12)+ 'names : '
     2394                values = 12*' '+'values: '
     2395                sigstr = 12*' '+'esds  : '
     2396                for item in orb[1]:
     2397                    pName = 'A%s%d:%d'%(item,iorb,AtId)
     2398                    names += item.rjust(10)
     2399                    values += '%10.5f'%orb[1][item][0]
     2400                    if pName in deformSig:
     2401                        sigstr += '%10.5f'%deformSig[pName]
     2402                    else:
     2403                        sigstr += 10*' '
     2404                pFile.write(names+'\n')
     2405                pFile.write(values+'\n')
    23202406                pFile.write(sigstr+'\n')
    23212407           
    … …  
    25732659            sigKey = {}
    25742660            wavesSig = {}
     2661            deformSig = {}
    25752662            cx,ct,cs,cia = General['AtomPtrs']
    25762663            for i,at in enumerate(Atoms):
    … …  
    26412728                                if pfx+name in sigDict:
    26422729                                    wavesSig[name] = sigDict[pfx+name]
     2730                                   
     2731            Deformations = Phase.get('Deformations',{})
     2732            for iAt in Deformations:
     2733                AtId = AtLookup[iAt]
     2734                DefAtm = Deformations[iAt]
     2735                for iorb,orb in enumerate(DefAtm):
     2736                    for parm in orb[1]:
     2737                        name = 'A%s%d:%d'%(parm,iorb,AtId)
     2738                        orb[1][parm][0] = parmDict[pfx+name]
     2739                        if pfx+name in sigDict:
     2740                            deformSig[name] = sigDict[pfx+name]
    26432741
    26442742            if pFile: PrintAtomsAndSig(General,Atoms,sigDict,sigKey)
     2743            if pFile and len(Deformations):
     2744                PrintDeformationsAndSig(General,Atoms,Deformations,deformSig)
    26452745            if pFile and General['Type'] == 'magnetic':
    26462746                PrintMomentsAndSig(General,Atoms,atomsSig)

  • TabularUnified trunk/GSASIIstrMain.py

    r5661 r5663  
    140140        if 'msg' not in Rvals: Rvals['msg'] = ''
    141141        Rvals['msg'] += msg
    142 
    143 def IgnoredLatticePrms(Phases):
    144     ignore = []
    145     copydict = {}
    146     for p in Phases:
    147         pfx = str(Phases[p]['pId']) + '::'
    148         laue = Phases[p]['General']['SGData']['SGLaue']
    149         axis = Phases[p]['General']['SGData']['SGUniq']
    150         if laue in ['-1',]:
    151             pass
    152         elif laue in ['2/m',]:
    153             if axis == 'a':
    154                 ignore += [pfx+'A4',pfx+'A5']
    155             elif axis == 'b':
    156                 ignore += [pfx+'A3',pfx+'A5']
    157             else:
    158                 ignore += [pfx+'A3',pfx+'A4']
    159         elif laue in ['mmm',]:
    160             ignore += [pfx+'A3',pfx+'A4',pfx+'A5']
    161         elif laue in ['4/m','4/mmm']:
    162             ignore += [pfx+'A1',pfx+'A3',pfx+'A4',pfx+'A5']
    163             copydict[pfx+'A0':[pfx+'A1']]
    164         elif laue in ['6/m','6/mmm','3m1', '31m', '3']:
    165             ignore += [pfx+'A1',pfx+'A3',pfx+'A4',pfx+'A5']
    166             copydict[pfx+'A0'] = [pfx+'A1',pfx+'A3']
    167         elif laue in ['3R', '3mR']:
    168             ignore += [pfx+'A1',pfx+'A2',pfx+'A4',pfx+'A5']
    169             copydict[pfx+'A0'] = [pfx+'A1',pfx+'A2']
    170             copydict[pfx+'A3'] = [pfx+'A4',pfx+'A5']
    171         elif laue in ['m3m','m3']:
    172             ignore += [pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
    173             copydict[pfx+'A0'] = [pfx+'A1',pfx+'A2']
    174     return ignore,copydict
    175 
    176 def AllPrmDerivs(Controls,Histograms,Phases,restraintDict,rigidbodyDict,
    177                      parmDict,varyList,calcControls,pawleyLookup,symHold,
    178                      dlg=None):
    179     '''Computes the derivative of the fitting function (total Chi**2) with
    180     respect to every parameter in the parameter dictionary (parmDict)
    181     by applying shift below the parameter value as well as above.
    182    
    183     :returns: a dict with the derivatives keyed by variable number.
    184       Derivatives are a list with three values: evaluated over
    185       v-d to v; v-d to v+d; v to v+d where v is the current value for the
    186       variable and d is a small delta value chosen for that variable type.
    187     '''
    188     import re
    189     rms = lambda y: np.sqrt(np.mean(y**2))
    190     G2mv.Map2Dict(parmDict,varyList)
    191     begin = time.time()
    192     seqList = Controls.get('Seq Data',[])
    193     hId = '*'
    194     if seqList:
    195         hId = str(Histograms[seqList[0]]['hId'])
    196 #    values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
    197 #    if np.any(np.isnan(values)):
    198 #        raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
    199     latIgnoreLst,latCopyDict = IgnoredLatticePrms(Phases)
    200     HistoPhases=[Histograms,Phases,restraintDict,rigidbodyDict]
    201     origDiffs = G2stMth.errRefine([],HistoPhases,parmDict,[],calcControls,pawleyLookup,None)
    202     chiStart = rms(origDiffs)
    203     origParms = copy.deepcopy(parmDict)
    204     #print('after 1st calc',time.time()-begin)
    205     derivCalcs = {}
    206     if dlg: dlg.SetRange(len(origParms))
    207     for i,prm in enumerate(origParms):
    208         if dlg:
    209             if not dlg.Update(i)[0]:
    210                 return None
    211         parmDict = copy.deepcopy(origParms)
    212         p,h,nam = prm.split(':')[:3]
    213         if hId != '*' and h != '' and h != hId: continue
    214         if (type(parmDict[prm]) is bool or type(parmDict[prm]) is str or
    215                  type(parmDict[prm]) is int): continue
    216         if type(parmDict[prm]) is not float and type(parmDict[prm]
    217                         ) is not np.float64:
    218             print('*** unexpected type for ',prm,parmDict[prm],type(parmDict[prm]))
    219             continue
    220         if prm in latIgnoreLst: continue # remove unvaried lattice params
    221         if re.match(r'\d:\d:D[012][012]',prm): continue   # don't need Dij terms
    222         if nam in ['Vol','Gonio. radius']: continue
    223         if nam.startswith('dA') and nam[2] in ['x','y','z']: continue
    224         delta = max(abs(parmDict[prm])*0.0001,1e-6)
    225         if nam in ['Shift','DisplaceX','DisplaceY',]:
    226             delta = 0.1
    227         elif nam.startswith('AUiso'):
    228             delta = 1e-5
    229         if nam[0] == 'A' and nam[1] in ['x','y','z']:
    230             dprm = prm.replace('::A','::dA')
    231             if dprm in symHold: continue # held by symmetry
    232             delta = 1e-6
    233         else:
    234             dprm = prm
    235         #print('***',prm,type(parmDict[prm]))
    236         #origVal = parmDict[dprm]
    237         parmDict[dprm] -= delta
    238         G2mv.Dict2Map(parmDict)
    239         if dprm in latCopyDict:         # apply contraints on lattice parameters
    240             for i in latCopyDict:
    241                 parmDict[i] = parmDict[dprm]
    242         #for i in parmDict:
    243         #    if origParms[i] != parmDict[i]: print('changed',i,origParms[i],parmDict[i])
    244         chiLow = rms(G2stMth.errRefine([],HistoPhases,parmDict,[],calcControls,pawleyLookup,None))
    245         parmDict[dprm] += 2*delta
    246         G2mv.Dict2Map(parmDict)
    247         if dprm in latCopyDict:         # apply contraints on lattice parameters
    248             for i in latCopyDict:
    249                 parmDict[i] = parmDict[dprm]
    250         #for i in parmDict:
    251         #    if origParms[i] != parmDict[i]: print('changed',i,origParms[i],parmDict[i])
    252         chiHigh = rms(G2stMth.errRefine([],HistoPhases,parmDict,[],calcControls,pawleyLookup,None))
    253         #print('===>',prm,parmDict[dprm],delta)
    254         #print(chiLow,chiStart,chiHigh)
    255         #print((chiLow-chiStart)/delta,0.5*(chiLow-chiHigh)/delta,(chiStart-chiHigh)/delta)
    256         derivCalcs[prm] = ((chiLow-chiStart)/delta,0.5*(chiLow-chiHigh)/delta,(chiStart-chiHigh)/delta)
    257     print('derivative computation time',time.time()-begin)
    258     return derivCalcs
    259142
    260143def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
    … …  
    429312    return IfOK,Rvals,result,covMatrix,sig,Lastshft
    430313
    431 def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None,newLeBail=False,allDerivs=False):
     314def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None,newLeBail=False):
    432315    '''Global refinement -- refines to minimize against all histograms.
    433316    This can be called in one of three ways, from :meth:`GSASIIdataGUI.GSASII.OnRefine` in an
    … …  
    467350    if allDerivs: #=============  develop partial derivative map
    468351        symHold = []
    469     (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,BLtables,MFtables,
     352    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,ORBtables,BLtables,MFtables,
    470353         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile,symHold=symHold)
    471354    calcControls['atomIndx'] = atomIndx
    … …  
    473356    calcControls['FFtables'] = FFtables
    474357    calcControls['EFtables'] = EFtables
     358    calcControls['ORBtables'] = ORBtables
    475359    calcControls['BLtables'] = BLtables
    476360    calcControls['MFtables'] = MFtables
    … …  
    492376    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
    493377    msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
    494     if allDerivs: #=============  develop partial derivative map
    495         varyListStart = varyList
    496         varyList = None
    497378    if msg:
    498379        return False,{'msg':'Unable to interpret multiplier(s): '+msg}
    … …  
    510391    if 'FrozenList' not in Controls['parmFrozen']:
    511392        Controls['parmFrozen']['FrozenList'] = []
    512     if varyList is not None:
    513         parmFrozenList = Controls['parmFrozen']['FrozenList']
    514         frozenList = [i for i in varyList if i in parmFrozenList]
    515         if len(frozenList) != 0:
    516             varyList = [i for i in varyList if i not in parmFrozenList]
    517             G2fil.G2Print(
    518                 'Frozen refined variables (due to exceeding limits)\n\t:{}'
    519                 .format(frozenList))
     393    parmFrozenList = Controls['parmFrozen']['FrozenList']
     394    frozenList = [i for i in varyList if i in parmFrozenList]
     395    if len(frozenList) != 0:
     396        varyList = [i for i in varyList if i not in parmFrozenList]
     397        G2fil.G2Print(
     398            'Frozen refined variables (due to exceeding limits)\n\t:{}'
     399            .format(frozenList))
    520400       
    521     ifSeq = False   
     401    ifSeq = False
    522402    printFile.write('\n Refinement results:\n')
    523403    printFile.write(135*'-'+'\n')
    524404    Rvals = {}
    525405    G2mv.Dict2Map(parmDict)  # impose constraints initially
    526     if allDerivs: #=============  develop partial derivative map
    527         derivDict = AllPrmDerivs(Controls, Histograms, Phases, restraintDict,
    528                      rigidbodyDict, parmDict, varyList, calcControls,
    529                      pawleyLookup,symHold,dlg)
    530         return derivDict,varyListStart
    531406   
    532407    try:
    … …  
    904779        try:
    905780            errmsg,warnmsg,groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,
    906                                                                       seqHistNum=hId,raiseException=True)
     781                seqHistNum=hId,raiseException=True)
    907782            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
    908783            G2mv.normParms(parmDict)

  • TabularUnified trunk/GSASIIstrMath.py

    r5655 r5663  
    334334            if 'U' in RBObj['ThermalMotion'][0]:
    335335                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
    336                
    337 def MakeSpHarmFF(HKL,Bmat,SHCdict,Tdata,hType,FFtables,BLtables,FF,SQ,ifDeriv=False):
     336                               
     337def MakeSpHarmFF(HKL,Bmat,SHCdict,Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ,ifDeriv=False):
    338338    ''' Computes hkl dependent form factors & derivatives from spinning rigid bodies
    339339    :param array HKL: reflection hkl set to be considered
    340340    :param array Bmat: inv crystal to Cartesian transfomation matrix
    341     :param dict SHCdict: spin RB data
     341    :param dict SHCdict: RB spin/deformation data
    342342    :param array Tdata: atom type info
    343343    :param str hType: histogram type
    344344    :param dict FFtables: x-ray form factor tables
     345    :param dict ORBtables: x-ray orbital form factor tables
    345346    :param dict BLtables: neutron scattering lenghts
    346     :param array FF: form factors - will be modified by adding the spin RB spherical harmonics terms
     347    :param array FF: form factors - will be modified by adding the spin/deformation RB spherical harmonics terms
    347348    :param array SQ: 1/4d^2 for the HKL set
    348     :param bool ifDeriv: True if dFF/dR, dFF/dA, & dFF/dShC to be returned
     349    :param bool ifDeriv: True if dFF/dcoff to be returned
    349350   
    350351    :returns: dict dFFdS of derivatives if ifDeriv = True
    … …  
    358359    dFFdS = {}
    359360    atFlg = []
     361    Th,Ph = G2lat.H2ThPh(np.reshape(HKL,(-1,3)),Bmat,[1.,0.,0.,1.])
     362    SQR = np.repeat(SQ,HKL.shape[1])
    360363    for iAt,Atype in enumerate(Tdata):
    361364        if 'Q' in Atype:
    … …  
    375378            ThMk,PhMk = MakePolar([SHdat['Oa'],SHdat['Oi'],SHdat['Oj'],SHdat['Ok']-dp],QB)
    376379            QR = np.repeat(twopi*np.sqrt(4.*SQ),HKL.shape[1])     #refl Q for Bessel fxn
    377             SQR = np.repeat(SQ,HKL.shape[1])
    378380            FF[:,iAt] = 0.
    379381            ishl = 0
    … …  
    450452            dFFdS[Ojname] = dSHdOj
    451453            dFFdS[Okname] = dSHdOk
     454        elif iAt in SHCdict and 'X' in hType:
     455            orKeys = list(ORBtables[Atype].keys())
     456            orbs = SHCdict[iAt]
     457            atFlg.append(1.0)
     458            FFcore = G2el.ScatFac(ORBtables[Atype][orKeys[0]],SQR)    #core
     459            FFtot = np.zeros_like(FFcore)
     460            for orb in orbs:
     461                Ne = orbs[orb].get('Ne',1.0) # not there for non <j0> orbs
     462                kappa = orbs[orb]['kappa']
     463                SQk = SQR/kappa**2
     464                ff = Ne*G2el.ScatFac(ORBtables[Atype][orKeys[int(orb)+1]],SQk)
     465                dffdk = G2el.ScatFacDer(ORBtables[Atype][orKeys[int(orb)+1]],SQk)
     466                dSH = 0.0
     467                if '<j0>' in orKeys[int(orb)+1]:
     468                    dSH = 1.0
     469                for term in orbs[orb]:
     470                    if 'D(' in term:
     471                        item = term.replace('D','C')
     472                        SH = G2lat.KslCalc(item,Th,Ph)
     473                        FFtot += SH*orbs[orb][term]*ff
     474                        name = 'A%s%s:%d'%(term,orb,iAt)
     475                        dFFdS[name] = SH*ff
     476                        dSH += SH*orbs[orb][term]
     477                    elif 'Ne' in term:
     478                        name = 'ANe%s:%d'%(orb,iAt)
     479                        dFFdS[name] = ff/Ne
     480                    elif 'kappa' in term:
     481                        if 'j0' in orKeys[int(orb)+1]:
     482                            FFtot += ff
     483                name = 'Akappa%s:%d'%(orb,iAt)
     484                dFFdS[name] = -2.0*Ne*SQk*dSH*dffdk/kappa
     485            FF[:,iAt] = FFcore+FFtot
    452486        else:
    453487            atFlg.append(0.)
    … …  
    477511                cof = bits[2]
    478512                SHCdict[atid][shno][cof] = parmDict[parm]
     513        if 'AD(' in parm or 'Akappa' in parm or 'ANe' in parm:       #atom deformation parms
     514            items = parm.split(':')
     515            atid = int(items[-1])
     516            name = items[2][1:]    #strip 'A'
     517            if atid not in SHCdict:
     518                SHCdict[atid] = {}
     519            orb = name[-1]
     520            if orb not in SHCdict[atid]:
     521                SHCdict[atid][orb] = {}
     522            SHCdict[atid][orb][name[:-1]] = parmDict[parm] #[atom id][orb no.][deform. coef]
    479523    return SHCdict
    480524   
    … …  
    948992    EFtables = calcControls['EFtables']
    949993    BLtables = calcControls['BLtables']
     994    ORBtables = calcControls['ORBtables']
    950995    Amat,Bmat = G2lat.Gmat2AB(G)
    951996    Flack = 1.0
    … …  
    10361081        #this must be done here. NB: same place for non-spherical atoms; same math except no Bessel part.
    10371082        if len(SHCdict):
    1038             MakeSpHarmFF(Uniq,Bmat,SHCdict,Tdata,hType,FFtables,BLtables,FF,SQ)     #Not Amat!
     1083            MakeSpHarmFF(Uniq,Bmat,SHCdict,Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ)     #Not Amat!
    10391084        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    10401085        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
    … …  
    10891134    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    10901135    FFtables = calcControls['FFtables']
     1136    ORBtables = calcControls['ORBtables']
    10911137    BLtables = calcControls['BLtables']
    10921138    hType = calcControls[hfx+'histType']
    … …  
    11441190        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    11451191        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
    1146         Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
     1192        Uniq = np.inner(H,SGMT)             # array(nSGOp,3,3)
    11471193        Phi = np.inner(H,SGT)
    11481194        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    … …  
    11581204        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
    11591205        if len(SHCdict):
    1160             dffdsh,atFlg = MakeSpHarmFF(Uniq,Bmat,SHCdict,Tdata,hType,FFtables,BLtables,FF,SQ,True)
     1206            dffdsh,atFlg = MakeSpHarmFF(Uniq,Bmat,SHCdict,Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ,True)
    11611207            if len(dffdSH):
    11621208                for item in dffdSH:
    … …  
    12251271#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    12261272        #loop over atoms - each dict entry is list of derivatives for all the reflections
    1227     for i in range(len(Mdata)):
     1273    for i in range(len(Tdata)):         #loop over atoms
    12281274        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    12291275        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
    … …  
    12381284        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
    12391285        for item in dffdSH:
    1240             i = int(item.split(':')[1])
    1241             dFdvDict[pfx+'RBS'+item] = np.sum(dFdff[:,:,i]*np.reshape(dffdSH[item],(nRef,-1)),axis=-1)
     1286            if 'SH' in item or 'O' in item:
     1287                i = int(item.split(':')[1])
     1288                dFdvDict[pfx+'RBS'+item] = np.sum(dFdff[:,:,i]*np.reshape(dffdSH[item],(nRef,-1)),axis=-1)
     1289            else:
     1290                dFdvDict[pfx+item] = np.sum(dFdff[:,:,i]*np.reshape(dffdSH[item],(nRef,-1)),axis=-1)
    12421291    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
    12431292    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
Note: See TracChangeset for help on using the changeset viewer.