Changeset 5067


Ignore:
Timestamp:
Nov 5, 2021 10:11:25 PM (2 years ago)
Author:
toby
Message:

ISODISTORT mode/displacement fixes; allow indexing of all Phases where there are no used histograms (needs to make IndexAllIds? conditional elsewhere?)

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIconstrGUI.py

    r5065 r5067  
    11361136                    helptext += " is created from a linear combination of the following variables:\n"
    11371137                    for term in item[:-3]:
     1138                        m = term[0]
     1139                        if np.isclose(m,0): continue
    11381140                        #var = str(term[1])
    11391141                        var,explain,note,warnmsg = term[1].fmtVarByMode(seqmode,note,warnmsg)
    … …  
    11411143                        if len(eqString[-1]) > maxlen:
    11421144                            eqString.append(' ')
    1143                         m = term[0]
    11441145                        if eqString[-1] != '':
    11451146                            if m >= 0:
    … …  
    11511152                            eqString[-1] += '{:} '.format(var)
    11521153                        else:
    1153                             eqString[-1] += '{:g}*{:} '.format(m,var)
     1154                            eqString[-1] += '{:.3g}*{:} '.format(m,var)
    11541155                        varMean = G2obj.fmtVarDescr(var)
    1155                         helptext += '\n  {:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
    1156                     # Add ISODISTORT help items
     1156                        helptext += '\n  {:.5g} * {:} '.format(m,var) + " ("+ varMean + ")"
     1157                    # Add extra notes about this constraint (such as from ISODISTORT)
    11571158                    if '_Explain' in data:
    1158                         # this ignores the phase number. TODO: refactor that in?
    11591159                        hlptxt = None
    11601160                        try:
    11611161                            hlptxt = data['_Explain'].get(item[-3])
    11621162                        except TypeError:
    1163                             # note that phase RanId is item[-3].phase
     1163                            # Patch fixed Nov 2021. Older projects have phase RanId in
     1164                            # item[-3].phase rather than a properly formed G2VarObj
    11641165                            hlptxt = data['_Explain'].get(str(item[-3].phase)+item[-3].name)
    11651166                        if hlptxt:
    1166                             if helptext: helptext += '\n\nNote warning:\n'
    1167                             helptext += hlptxt
     1167                            helptext += '\n\n'+ hlptxt
    11681168                    if item[-3]:
    11691169                        typeString = '(NEWVAR) ' + str(item[-3]) + ' = '
    … …  
    11901190                            eqString[-1] += '{:} '.format(var)
    11911191                        else:
    1192                             eqString[-1] += '{:g}*{:} '.format(m,var)
     1192                            eqString[-1] += '{:.3g}*{:} '.format(m,var)
    11931193                        varMean = G2obj.fmtVarDescr(var)
    1194                         helptext += '\n  {:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
     1194                        helptext += '\n  {:.5g} * {:} '.format(m,var) + " ("+ varMean + ")"
    11951195                        helptext += explain
    11961196                    typeString = 'CONST'
    … …  
    12071207                    #var1 = str(item[1][1])
    12081208                    var1,explain,note,warnmsg = item[1][1].fmtVarByMode(seqmode,note,warnmsg)
    1209                     helptext += '\n  {:3g} * {:} '.format(m,var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
     1209                    helptext += '\n  {:.5g} * {:} '.format(m,var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
    12101210                    eqString[-1] += '{:} = {:}'.format(var1,var)
    12111211                    if m != 1:
    … …  
    12351235                            eqString[-1] += '{:}'.format(var)
    12361236                        else:
    1237                             eqString[-1] += '{:g}*{:} '.format(m,var)
    1238                         helptext += '\n  {:3g} * {:} '.format(m,var) + " ("+ varMean + ")"
     1237                            eqString[-1] += '{:.3g}*{:} '.format(m,var)
     1238                        helptext += '\n  {:.5g} * {:} '.format(m,var) + " ("+ varMean + ")"
    12391239                    eqString[-1] += ' = {:} '.format(indepterm)
    12401240                    typeString = 'EQUIV'
    … …  
    15291529    #     dlg.Destroy()
    15301530    #     return
    1531     G2obj.IndexAllIds(Histograms,Phases)
     1531    if Histograms and Phases:
     1532        G2obj.IndexAllIds(Histograms,Phases)
    15321533    # for p in Phases:
    15331534    #     if 'ISODISTORT' in Phases[p] and 'G2VarList' in Phases[p]['ISODISTORT']:
    … …  
    38303831        else:
    38313832            return
    3832    
     3833
    38333834    covdata = G2frame.GPXtree.GetItemPyData(
    38343835        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Covariance'))
    … …  
    38403841        if c[-1] != 'f' or not c[-3]: continue
    38413842        constrDict[str(c[-3])] = c
    3842 
    3843     parmDict,varyList = G2frame.MakeLSParmDict()
    38443843
    38453844    dlg = wx.Dialog(G2frame,wx.ID_ANY,'ISODISTORT mode values',#size=(630,400),
    … …  
    38543853        dlg, wx.ID_ANY,#size=(100,200),
    38553854        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
    3856     subSizer1 = wx.FlexGridSizer(cols=2,hgap=5,vgap=2)
     3855    subSizer1 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
    38573856    panel2 = wxscroll.ScrolledPanel(
    38583857        dlg, wx.ID_ANY,#size=(100,200),
    38593858        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
    3860     subSizer2 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
    3861     subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'Parameter name  '))
     3859    subSizer2 = wx.FlexGridSizer(cols=4,hgap=5,vgap=2)
     3860    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'ISODISTORT\nname'))
     3861    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'GSAS-II\nname'))
    38623862    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
     3863    for i in range(3): subSizer1.Add((-1,5)) # spacer
    38633864    subSizer2.Add((-1,-1))
    3864     subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'Mode name  '))
     3865    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'ISODISTORT\nMode name'))
     3866    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'GSAS-II\nname'))
    38653867    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
     3868    for i in range(4): subSizer2.Add((-1,5))
    38663869    # ISODISTORT displacive modes
    38673870    if 'G2VarList' in ISO:
    3868         deltaList = []
    3869         for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
    3870             dvar = gv.varname()
    3871             var = dvar.replace('::dA','::A')
    3872             albl = Ilbl[:Ilbl.rfind('_')]
    3873             v = Ilbl[Ilbl.rfind('_')+1:]
    3874             pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
    3875             if var in parmDict:
    3876                 cval = parmDict[var][0]
    3877             else:
    3878                 cval = 0.0
    3879             deltaList.append(cval-pval)
    3880         modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
    3881         for lbl,xyz,var,val,norm,G2mode in zip(
    3882                 ISO['IsoVarList'],deltaList,
    3883                 ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
     3871        dispVals,dispSUs,modeVals,modeSUs = G2mth.CalcIsoDisp(Phases[phase],covdata=covdata)
     3872        for (lbl,xyz,xyzsig,G2var,
     3873             var,mval,msig,G2mode) in zip(
     3874                ISO['IsoVarList'],dispVals,dispSUs,ISO['G2VarList'],
     3875                ISO['IsoModeList'],modeVals,modeSUs,ISO['G2ModeList'] ):
    38843876            if str(G2mode) in constrDict:
    38853877                ch = G2G.HelpButton(panel2,fmtHelp(constrDict[str(G2mode)],var))
    … …  
    38883880                subSizer2.Add((-1,-1))
    38893881            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
     3882            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(G2var)))
    38903883            try:
    3891                 value = G2py3.FormatSigFigs(xyz)
     3884                value = G2mth.ValEsd(xyz,xyzsig)
    38923885            except TypeError:
    38933886                value = str(xyz)           
    38943887            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
     3888
    38953889            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
     3890            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(G2mode)))
    38963891            try:
    3897                 value = G2py3.FormatSigFigs(val/norm)
    3898                 if 'varyList' in covdata:
    3899                     if str(G2mode) in covdata['varyList']:
    3900                         sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
    3901                         value = G2mth.ValEsd(val/norm,sig/norm)
     3892                value = G2mth.ValEsd(mval,msig)
    39023893            except TypeError:
    3903                 value = '?'
     3894                value = str(mval)
    39043895            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
    39053896    # ISODISTORT occupancy modes
    39063897    if 'G2OccVarList' in ISO:
    39073898        deltaList = []
     3899        parmDict,varyList = G2frame.MakeLSParmDict()
    39083900        for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
    39093901            var = gv.varname()

  • trunk/GSASIIdataGUI.py

    r5066 r5067  
    44384438                self.EnableRefineCommand()
    44394439                self.init_vars()
     4440                G2obj.IndexAllIds({},{}) # clear old index info
    44404441        finally:
    44414442            dlg.Destroy()
    … …  
    50465047        G2obj.IndexAllIds(Histograms=Histograms,Phases=phaseData)
    50475048        return Histograms,Phases
    5048        
     5049   
    50495050    def MakeLSParmDict(self,seqHist=None):
    50505051        '''Load all parameters used for computation from the tree into a
    … …  
    64396440       
    64406441        # Phase/ ISODISTORT tab
    6441         G2G.Define_wxId('wxID_ISODNEWPHASE',)       
    64426442        self.ISODData = wx.MenuBar()
    64436443        self.PrefillDataMenu(self.ISODData)
    … …  
    64456445        self.ISODDataEdit = wx.Menu(title='')
    64466446        self.ISODData.Append(menu=self.ISODDataEdit, title='Operations')
    6447         self.ISODDataEdit.Append(G2G.wxID_ISODNEWPHASE,'Make cif file','From ISODISTORT selection')
     6447        G2G.Define_wxId('wxID_ISODNEWPHASE')
     6448        self.ISODDataEdit.Append(G2G.wxID_ISODNEWPHASE,'Make CIF file','From ISODISTORT selection')
     6449        G2G.Define_wxId('wxID_SHOWISO1')
     6450        self.ISODDataEdit.Append(G2G.wxID_SHOWISO1,'Show ISODISTORT modes',
     6451                'Show ISODISTORT mode values for all phases')
    64486452        self.PostfillDataMenu()
    64496453

  • trunk/GSASIImath.py

    r5060 r5067  
    20702070    return drawAtoms,Fade
    20712071
     2072def patchIsoDisp(ISO):
     2073    '''patch: look for older ISODISTORT imports (<Nov 2021)'''
     2074    print('''
     2075======================================================================
     2076Warning: The ISODISTORT modes were read before the importer
     2077was corrected to save displacement offsets. Will attempt to correct
     2078from ParentStructure (correct only if displacivemode values are all
     2079zero in initial CIF.) Reimporting is suggested.
     2080======================================================================
     2081''')
     2082    ISO['G2coordOffset'] = []
     2083    for Ilbl in ISO['IsoVarList']:
     2084        albl = Ilbl[:Ilbl.rfind('_')]
     2085        v = Ilbl[Ilbl.rfind('_')+1:]
     2086        ISO['G2coordOffset'].append(
     2087                ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
     2088                ) # this will be wrong if the _iso_deltacoordinate_value are not zero
     2089
     2090def CalcIsoDisp(Phase,parmDict={},covdata={}):
     2091    '''Compute the ISODISTORT displacement variable values from the
     2092    atomic coordinates, applying the p::dA?:n displacements if parmDict
     2093    is supplied. Uncertainties are computed if covdata is supplied.
     2094    '''
     2095    ISO = Phase['ISODISTORT']
     2096    if 'G2coordOffset' not in ISO: patchIsoDisp(ISO) # patch Nov 2021
     2097    atmIndex = {a[-1]:i for i,a in enumerate(Phase['Atoms'])}
     2098    cx,ct,cs,ci = getAtomPtrs(Phase,False)
     2099    coordOffset = {xyz:cx+i for i,xyz in enumerate(('x','y','z'))}
     2100    # get uncertainties on modes and compute them for displacements
     2101    if 'varyList' in covdata:
     2102        modes = [i.name for i in ISO['G2ModeList']]
     2103        covDict = dict(zip(covdata.get('varyList',[]),covdata.get('sig',[])))
     2104        modeEsds = [covDict.get(str(g2),-1) for g2 in modes]
     2105        vcov = getVCov(modes,covdata['varyList'],covdata['covMatrix'])
     2106        normMode2Var = ISO['Mode2VarMatrix']*ISO['NormList']
     2107        # to get displacements from modes use np.dot(normMode2Var,vector-of-mode-values)
     2108        sigMat = np.inner(normMode2Var,np.inner(normMode2Var,vcov))
     2109        var = np.diag(sigMat)
     2110        dispEsds = list(np.where(var>0.,np.sqrt(var),-0.0001))
     2111    else:
     2112        dispEsds = len(ISO['G2VarList'])*[-0.0001]
     2113        modeEsds = len(ISO['G2ModeList'])*[-0.0001]
     2114
     2115    dispValues = []
     2116    for iso,g2,off in zip(ISO['IsoVarList'],ISO['G2VarList'],ISO['G2coordOffset']):
     2117        if g2.atom not in atmIndex:
     2118            print('Atom not found in atom list',g2)
     2119            return [],[],[],[]
     2120        atm = Phase['Atoms'][atmIndex[g2.atom]]
     2121        pos = atm[coordOffset[g2.name[-1]]] + parmDict.get(str(g2),0.0) - off
     2122        dispValues.append(pos)
     2123
     2124    modeValues = np.dot(ISO['Var2ModeMatrix'],dispValues) / ISO['NormList']
     2125       
     2126    return dispValues,dispEsds,modeValues,modeEsds
     2127
     2128def CalcIsoCoords(Phase,parmDict,covdata={}):
     2129    '''Compute the coordinate positions from ISODISTORT displacement mode values
     2130    Uncertainties are computed if covdata is supplied.
     2131
     2132    :returns: modeDict,posDict where modeDict contains pairs of mode values
     2133    and mode s.u. values; posDict contains pairs of displacement values
     2134    and their s.u. values.
     2135    '''
     2136    ISO = Phase['ISODISTORT']
     2137    if 'G2coordOffset' not in ISO: patchIsoDisp(ISO) # patch Nov 2021
     2138    modes = [i.name for i in ISO['G2ModeList']]  # modes from the parmDict
     2139    normMode2Var = ISO['Mode2VarMatrix']*ISO['NormList']
     2140    modeVals = []
     2141    for i in modes:
     2142        if i not in parmDict:
     2143            print('Mode ',i,'not defined in the parameter dictionary')
     2144            return {},{}
     2145        try:
     2146            modeVals.append(float(parmDict[i][0]))
     2147        except:
     2148            modeVals.append(float(parmDict[i]))
     2149    if 'varyList' in covdata:
     2150        covDict = dict(zip(covdata.get('varyList',[]),covdata.get('sig',[])))
     2151        modeEsds = [covDict.get(str(g2),-1) for g2 in modes]
     2152        vcov = getVCov(modes,covdata['varyList'],covdata['covMatrix'])
     2153        sigMat = np.inner(normMode2Var,np.inner(normMode2Var,vcov))
     2154        var = np.diag(sigMat)
     2155        dispEsds = list(np.where(var>0.,np.sqrt(var),-0.0001))
     2156    else:
     2157        modeEsds = len(ISO['G2ModeList'])*[-0.0001]
     2158        dispEsds = len(ISO['G2VarList'])*[-0.0001]       
     2159    dispValues =   np.dot(normMode2Var,modeVals)
     2160    modeDict = {str(g2):([val,esd]) for val,g2,esd in
     2161                        zip(modeVals,ISO['G2ModeList'],modeEsds)}
     2162    posDict = {str(g2).replace('::dA','::A'):([dsp,esd]) for dsp,g2,off,esd in
     2163                        zip(dispValues,ISO['G2VarList'],ISO['G2coordOffset'],dispEsds)}
     2164    return modeDict,posDict
     2165
    20722166def ApplyModeDisp(data):
    2073     ''' Applies ISODISTORT mode displacements to drawing atoms
     2167    ''' Applies ISODISTORT mode displacements to drawing atoms.
     2168    This changes the contents of the Draw Atoms positions but not
     2169    the Atoms positions.
     2170
     2171    :param dict data: the contents of the Phase data tree item for a
     2172      particular phase
    20742173    '''
    20752174    generalData = data['General']

  • trunk/GSASIIobj.py

    r5062 r5067  
    11721172
    11731173Displacive modes are a bit complex in that they relate to delta displacements, relative to an offset value for each coordinate,
    1174 and because the modes are normalized, while GSAS-II also uses displacements, but these are added to the coordinates after
    1175 each refinement cycle and the delta values are set to zero. ISODISTORT uses fixed offsets (subtracted from the actual position
    1176 to obtain the delta values) that are taken from ``_iso_coordinate_formula`` and these are placed in
    1177 ``.Phase['ISODISTORT']['ParentStructure]`` (keyed by atom label). The normalization factors (which the delta values are divided by)
    1178 are taken from ``_iso_displacivemodenorm_value`` and are placed in ``.Phase['ISODISTORT']['NormList']`` in the same order as as
    1179 ``...['IsoModeList']`` and ``...['G2ModeList']``.
     1174and because the modes are normalized. While GSAS-II also uses displacements,  these are added to the coordinates after
     1175each refinement cycle and then the delta values are set to zero.
     1176ISODISTORT uses fixed offsets (subtracted from the actual position
     1177to obtain the delta values) that are taken from the parent structure coordinate and the initial offset value
     1178(in ``_iso_deltacoordinate_value``) and these are placed in
     1179``.Phase['ISODISTORT']['G2coordOffset']`` in the same order as ``.Phase['ISODISTORT']['G2ModeList']`` and
     1180``.Phase['ISODISTORT']['IsoVarList']``.
     1181
     1182The normalization factors (which the delta values are divided by)
     1183are taken from ``_iso_displacivemodenorm_value`` and are placed in ``.Phase['ISODISTORT']['NormList']`` in the same
     1184order as as ``...['IsoModeList']`` and ``...['G2ModeList']``.
    11801185
    11811186The CIF contains a sparse matrix, from the ``loop_`` containing ``_iso_displacivemodematrix_value`` which provides the equations
    11821187for determining the mode values from the coordinates, that matrix is placed in ``.Phase['ISODISTORT']['Mode2VarMatrix']``.
    11831188The matrix is inverted to produce ``.Phase['ISODISTORT']['Var2ModeMatrix']``, which determines how to compute the
    1184 mode values from the delta coordinate values. These values are used for the in :func:`GSASIIconstrGUI.ShowIsoDistortCalc` 
     1189mode values from the delta coordinate values. These values are used for the in :func:`GSASIIconstrGUI.ShowIsoDistortCalc`,
    11851190which shows coordinate and mode values, the latter with s.u. values.
    1186 
    11871191
    11881192------------------------------

  • trunk/GSASIIphsGUI.py

    r5066 r5067  
    44284428            thresh=[[8.0,6.0],[17.191,11.527]],pickHandler=pickHandler)
    44294429
    4430     def OnIsoDistortCalc(event):
     4430    def OnShowIsoDistortCalc(event):
     4431        Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
     4432        if Histograms and Phases:
     4433            G2obj.IndexAllIds(Histograms,Phases)
    44314434        G2cnstG.ShowIsoDistortCalc(G2frame,data['General']['Name'])
    44324435           
    … …  
    67376740           
    67386741       
    6739 #### ISOTDISTORT results tab ###############################################################################
     6742#### ISODISTORT results tab ###############################################################################
    67406743
    67416744    def UpdateISODISTORT(Scroll=0):
    … …  
    67736776               
    67746777        def OnDispl(event):
     6778            '''Respond to movement of distortion mode slider'''
    67756779            Obj = event.GetEventObject()
    67766780            idsp,dispVal = Indx[Obj.GetId()]
    … …  
    67846788           
    67856789        def OnDispVal(invalid,value,tc):
     6790            '''Respond to entry of a value into a distortion mode entry widget'''
    67866791            idsp,displ = Indx[tc.GetId()]
    67876792            displ.SetValue(int(value*1000)+100)
    … …  
    67936798           
    67946799        def OnReset(event):
     6800            '''Reset all distortion mode values to zero'''
    67956801            data['ISODISTORT']['modeDispl'] = np.zeros(len(data['ISODISTORT']['G2ModeList']))
    67966802            err = G2mth.ApplyModeDisp(data)
    … …  
    68026808       
    68036809        Indx = {}     
     6810        G2frame.dataWindow.ISODDataEdit.Enable(G2G.wxID_ISODNEWPHASE,
     6811                                'rundata' in data['ISODISTORT'])
     6812        G2frame.dataWindow.ISODDataEdit.Enable(G2G.wxID_SHOWISO1,
     6813                    ('G2VarList' in data['ISODISTORT']) or
     6814                    ('G2OccVarList' in data['ISODISTORT']))
    68046815        if 'radio' not in data['ISODISTORT']:
    68056816            if not data['ISODISTORT']:
    … …  
    69076918       
    69086919    def OnNewISOPhase(event):
    6909         ''' Make cif file with ISODISTORT
     6920        ''' Make CIF file with ISODISTORT
    69106921        '''
    69116922        import ISODISTORT as ISO
    69126923
    6913         if data['ISODISTORT']:
     6924        if 'rundata' in data['ISODISTORT'] and data['ISODISTORT']['selection'] is not None:
    69146925            CIFfile = ISO.GetISODISTORTcif(data)
    69156926            G2G.G2MessageBox(G2frame,'ISODISTORT generated cif file %s has been created.'%CIFfile)
     6927        elif 'rundata' in data['ISODISTORT']:
     6928            G2G.G2MessageBox(G2frame,'Need to select an ISODISTORTdistortion model first before creating a CIF')
    69166929        else:
    69176930            G2G.G2MessageBox(G2frame,'ERROR - need to run ISODISTORT first - see General/Compute menu')       
    … …  
    1356313576        G2frame.Bind(wx.EVT_MENU, OnFracSplit, id=G2G.wxID_ATOMFRACSPLIT)
    1356413577        G2frame.Bind(wx.EVT_MENU, OnDensity, id=G2G.wxID_ATOMSDENSITY)
    13565         G2frame.Bind(wx.EVT_MENU, OnIsoDistortCalc, id=G2G.wxID_ISODISP)
     13578        G2frame.Bind(wx.EVT_MENU, OnShowIsoDistortCalc, id=G2G.wxID_ISODISP)
    1356613579        if 'HydIds' in data['General']:
    1356713580            G2frame.dataWindow.AtomEdit.Enable(G2G.wxID_UPDATEHATOM,True)
    … …  
    1364813661        FillSelectPageMenu(TabSelectionIdDict, G2frame.dataWindow.ISODData)
    1364913662        G2frame.Bind(wx.EVT_MENU, OnNewISOPhase, id=G2G.wxID_ISODNEWPHASE)
     13663        G2frame.Bind(wx.EVT_MENU, OnShowIsoDistortCalc, id=G2G.wxID_SHOWISO1)
    1365013664        # MC/SA
    1365113665        FillSelectPageMenu(TabSelectionIdDict, G2frame.dataWindow.MCSAMenu)

  • trunk/GSASIIstrIO.py

    r5051 r5067  
    321321                albl = Ilbl[:Ilbl.rfind('_')]
    322322                v = Ilbl[Ilbl.rfind('_')+1:]
    323                 pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
     323                #pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
    324324                if var in parmDict:
    325325                    cval = parmDict[var]

  • trunk/imports/G2phase_CIF.py

    r5065 r5067  
    709709            coordVarLbl = []
    710710            G2varObj = []
     711            G2coordOffset = []
    711712            idlist = []
    712             error = False
    713             for id,lbl in zip(
    714                 blk.get('_iso_deltacoordinate_ID'),
    715                 blk.get('_iso_deltacoordinate_label') ):
    716                 idlist.append(int(id))
    717                 coordVarLbl.append(lbl)
    718                 if '_' in lbl:
    719                     albl = lbl[:lbl.rfind('_')]
    720                     vlbl = lbl[lbl.rfind('_')+1:]
    721                 else:
    722                     self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
    723                     error = True
    724                     continue
    725                 if albl not in atomlbllist:
    726                     self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
    727                     error = True
    728                     continue
    729                 # else:
    730                 #     anum = atomlbllist.index(albl)
    731                 var = varLookup.get(vlbl)
    732                 if not var:
    733                     self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
    734                     error = True
    735                     continue
    736                 G2varObj.append(G2obj.G2VarObj(
    737                     (self.Phase['ranId'],None,var,ranIdlookup[albl])
    738                     ))
    739             if error:
    740                 raise Exception("Error decoding variable labels")
    741             # just in case the items are not ordered increasing by id, sort them here
    742             coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])]
    743             G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
    744 
    745             if len(G2varObj) != len(modelist):
    746                 print ("non-square input")
    747                 raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
    748 
    749713            error = False
    750714            ParentCoordinates = {}
    … …  
    779743                print (self.warnings)
    780744                raise Exception("Error decoding variable labels")
     745           
     746            error = False
     747            coordOffset = {xyz:i for i,xyz in enumerate(('dx','dy','dz'))}
     748            for id,lbl,val in zip(
     749                blk.get('_iso_deltacoordinate_ID'),
     750                blk.get('_iso_deltacoordinate_label'),
     751                blk.get('_iso_deltacoordinate_value')
     752                ):
     753                idlist.append(int(id))
     754                coordVarLbl.append(lbl)
     755                if '_' in lbl:
     756                    albl = lbl[:lbl.rfind('_')]
     757                    vlbl = lbl[lbl.rfind('_')+1:]
     758                else:
     759                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
     760                    error = True
     761                    continue
     762                if albl not in atomlbllist:
     763                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
     764                    error = True
     765                    continue
     766                # else:
     767                #     anum = atomlbllist.index(albl)
     768                var = varLookup.get(vlbl)
     769                if not var:
     770                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
     771                    error = True
     772                    continue
     773                G2coordOffset.append(ParentCoordinates[albl][coordOffset[vlbl]] - float(val))
     774                G2varObj.append(G2obj.G2VarObj(
     775                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
     776                    ))
     777            if error:
     778                raise Exception("Error decoding variable labels")
     779            # just in case the items are not ordered increasing by id, sort them here
     780            coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])]
     781            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
     782
     783            if len(G2varObj) != len(modelist):
     784                print ("non-square input")
     785                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
     786
    781787            # get mapping of modes to atomic coordinate displacements
    782788            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
    … …  
    787793                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
    788794            # Invert to get mapping of atom displacements to modes
    789             displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
     795            Var2ModeMatrix = np.linalg.inv(displacivemodematrix)
    790796            # create the constraints
    791797            modeVarList = []
    792798            modeDispl = []
    793             for i,row in enumerate(displacivemodeInvmatrix):
     799            for i,row in enumerate(Var2ModeMatrix):
    794800                constraint = []
    795801                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
    … …  
    802808                self.Constraints.append(constraint)
    803809                modeDispl.append(0.0)
    804             # normilization constants
     810            # normalization constants
    805811            normlist = []
    806812            idlist = []
    … …  
    819825                'IsoVarList' : coordVarLbl,
    820826                'G2VarList' : G2varObj,
    821                 'ParentStructure' : ParentCoordinates,
     827                'G2coordOffset' : G2coordOffset,
    822828                # mode items
    823829                'IsoModeList' : modelist,
    … …  
    825831                'NormList' : normlist,
    826832                # transform matrices
    827                 'Var2ModeMatrix' : displacivemodeInvmatrix,
     833                'Var2ModeMatrix' : Var2ModeMatrix,
    828834                'Mode2VarMatrix' : displacivemodematrix,
    829835                'modeDispl' : modeDispl
    … …  
    831837            # make explaination dictionary
    832838            for mode,shortmode in zip(modelist,shortmodelist):
    833                 explaination[str(self.Phase['ranId'])+shortmode] = (
    834                     "ISODISTORT full name "+str(mode))
     839                modeVar = G2obj.G2VarObj(
     840                    (self.Phase['ranId'],None,shortmode,None))
     841                explaination[modeVar] = ("Full ISODISTORT name for " +
     842                                            shortmode + " is " + str(mode))
    835843        #----------------------------------------------------------------------
    836844        # now read in the ISODISTORT occupancy modes
    … …  
    965973            # make explaination dictionary
    966974            for mode,shortmode in zip(modelist,shortmodelist):
    967                 explaination[str(self.Phase['ranId'])+shortmode] = (
    968                     "ISODISTORT full name "+str(mode))
     975                modeVar = G2obj.G2VarObj(
     976                    (self.Phase['ranId'],None,shortmode,None))
     977                explaination[modeVar] = ("Full ISODISTORT name for " +
     978                                            shortmode + " is " + str(mode))
     979        if explaination: self.Constraints.append(explaination)
    969980        #----------------------------------------------------------------------
    970981        # done with read
    971982        #----------------------------------------------------------------------
    972         if explaination: self.Constraints.append(explaination)
    973 
    974         if not debug: return
     983
     984        def fmtEqn(i,head,l,var,k):
     985            'format a section of a row of variables and multipliers'
     986            if k == 0: return head,l
     987            if len(head) + len(l) > 65:
     988                print(head+l)
     989                head = 20*' '
     990                l = ''
     991            if k < 0 and i > 0:
     992                l += ' - '
     993                k = -k
     994            elif i > 0:
     995                l += ' + '
     996            if k == 1:
     997                l += '%s ' % str(var)
     998            else:
     999                l += '%.3f * %s' % (k,str(var))
     1000            return head,l
    9751001
    9761002        # debug: show displacive mode var to mode relations
    977         if 'IsoVarList' in self.Phase['ISODISTORT']:
    978             # coordinate items
    979             #coordVarLbl = self.Phase['ISODISTORT']['IsoVarList']
    980             G2varObj = self.Phase['ISODISTORT']['G2VarList']
    981             #ParentCoordinates = self.Phase['ISODISTORT']['ParentStructure']
    982             # mode items
    983             modelist = self.Phase['ISODISTORT']['IsoModeList']
    984             modeVarList = self.Phase['ISODISTORT']['G2ModeList']
    985             normlist = self.Phase['ISODISTORT']['NormList']
    986             # transform matrices
    987             #displacivemodeInvmatrix = self.Phase['ISODISTORT']['Var2ModeMatrix']
    988             #displacivemodematrix = self.Phase['ISODISTORT']['Mode2VarMatrix']
    989        
    990             print( 70*'=')
    991             print('\nVar2ModeMatrix' ,'IsoVarList' )
    992             def fmtConstr(i,head,l,var,k):
    993                 if len(head) + len(l) > 65:
    994                     print(head+l)
    995                     head = 20*' '
    996                     l = ''
    997 #                if k == 0: return
    998                 if k < 0 and i > 0:
    999                     l += '  -  '
    1000                     k = -k
    1001                 elif i > 0:
    1002                     l += '  +  '
    1003                 if k == 1:
    1004                     l += '%s ' % str(var)
    1005                 else:
    1006                     l += '%.3f * %s' % (k,str(var))
    1007                 return head,l
    1008             for i,row in enumerate(displacivemodeInvmatrix):
    1009                 head = str(i) + ': ' + str(modeVarList[i]) + ' = '
     1003        if debug and 'IsoVarList' in self.Phase['ISODISTORT']:
     1004            print('\n' + 70*'=')
     1005            print('ISO modes from Iso coordinate vars (using Var2ModeMatrix, IsoVarList, G2VarList & G2ModeList)' )
     1006            for i,row in enumerate(self.Phase['ISODISTORT']['Var2ModeMatrix']):
     1007                head = '  ' + str(self.Phase['ISODISTORT']['G2ModeList'][i]) + ' = '
    10101008                line = ''
    10111009                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
    1012                     head,line = fmtConstr(j,head,line,G2varObj[j].name,k)
     1010                    var = self.Phase['ISODISTORT']['IsoVarList'][j]
     1011                    head,line = fmtEqn(j,head,line,var,k)
     1012                print(head+line)
     1013            print('\nConstraints')
     1014            for c in self.Constraints:
     1015                if type(c) is dict: continue
     1016                if c[-1] != 'f': continue
     1017                line = ''
     1018                head = '  ' + str(c[-3]) + ' = '
     1019                for j,(k,var) in enumerate(c[:-3]):
     1020                    head,line = fmtEqn(j,head,line,var,k)
    10131021                print(head+line)
    10141022
    … …  
    10251033                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
    10261034
     1035            print('\n' + 70*'=')
     1036            print('Confirming inverse mode relations computed from displacement values',
     1037                      '\nusing Var2ModeMatrix, NormList, IsoVarList')
     1038            # compute the mode values from the reported coordinate deltas
     1039            for i,(row,n) in enumerate(zip(self.Phase['ISODISTORT']['Var2ModeMatrix'],
     1040                                           self.Phase['ISODISTORT']['NormList'])):
     1041                line = ''
     1042                head = str(self.Phase['ISODISTORT']['G2ModeList'][i])+' = ('
     1043                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
     1044                    head,line = fmtEqn(j,head,line,lbl,k)
     1045                print(head+line+') / '+('%.3f'%n))
     1046                line = ''
     1047                head = str(self.Phase['ISODISTORT']['G2ModeList'][i])+' = ('
     1048                vsum = 0.
     1049                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
     1050                    val = "{:3g}".format(coordVarDelta[lbl])
     1051                    head,line = fmtEqn(j,head,line,val,k)
     1052                    vsum += coordVarDelta[lbl] * k
     1053                print(head+line+') / '+('%.3f'%n))
     1054                fileval = modeVarDelta[self.Phase['ISODISTORT']['IsoModeList'][i]]
     1055                print("{} = {:4g} (value read from CIF = {:4g})\n".format(
     1056                    self.Phase['ISODISTORT']['G2ModeList'][i], vsum, fileval))
     1057
    10271058            print( 70*'=')
    1028             print('\nInverse relations using Var2ModeMatrix, NormList, IsoVarList')
    1029             # compute the mode values from the reported coordinate deltas
    1030             for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
    1031                 line = ''
    1032                 head = 'a'+str(i)+': '+str(modeVarList[i])+' = ('
    1033                 for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
    1034                     head,line = fmtConstr(j,head,line,lbl,k)
    1035                 print(head+line+') / '+('%.3f'%n))
    1036             print('\nCalculation checks\n')
    1037             for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
    1038                 #l = ''
    1039                 sl = ''
    1040                 s = 0.
    1041                 for lbl,k in zip(coordVarLbl,row):
     1059            print('Direct displacement relations computed from ISO modes in CIF',
     1060                      '\nusing Mode2VarMatrix, NormList, IsoModeList, IsoVarList',)
     1061            # compute the coordinate displacements from the reported mode values
     1062            for lbl,row in zip(self.Phase['ISODISTORT']['IsoVarList'],
     1063                               self.Phase['ISODISTORT']['Mode2VarMatrix']):
     1064                l = ''
     1065                s = 0.0
     1066                head = lbl+' ='
     1067                for j,(k,n) in enumerate(zip(row,self.Phase['ISODISTORT']['NormList'])):
    10421068                    if k == 0: continue
    1043                     #if l: l += ' + '
    1044                     #l += lbl+' * '+str(k)
    1045                     if sl: sl += ' + '
    1046                     sl += str(coordVarDelta[lbl])+' * '+str(k)
    1047                     s += coordVarDelta[lbl] * k
    1048                 print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
    1049                 print(' ?= ',modeVarDelta[modelist[i]])
    1050                 print()
    1051 
    1052             print( 70*'=')
    1053             print('\nDirect relations using Mode2VarMatrix, NormList, IsoVarList')
    1054             # compute the coordinate displacements from the reported mode values
     1069                    if len(l) > 65:
     1070                        print(head,l)
     1071                        head = 20*' '
     1072                        l = ''
     1073                    l1 = ''
     1074                    k1 = k
     1075                    if j > 0 and k < 0:
     1076                        k1 = -k
     1077                        l1 = ' - '
     1078                    elif j > 0:
     1079                        l1 += ' + '
     1080                    l += '{:} {:3g} * {:4g} * {:}'.format(
     1081                        l1, k1, n, self.Phase['ISODISTORT']['G2ModeList'][j])
     1082                   
     1083                    s += n * modeVarDelta[self.Phase['ISODISTORT']['IsoModeList'][j]] * k
     1084                print(head,l)
     1085                print(lbl,'=',s)
     1086                print(lbl,'==>',str(self.Phase['ISODISTORT']['G2VarList'][i]),'\n')
    10551087            DeltaCoords = {}
    10561088            for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
    1057                 l = ''
    10581089                s = 0.0
     1090                for j,(k,n) in enumerate(zip(row,self.Phase['ISODISTORT']['NormList'])):
     1091                    s += k * n * modeVarDelta[self.Phase['ISODISTORT']['IsoModeList'][j]]
    10591092                at,d = lbl.rsplit('_',1)
    10601093                if at not in DeltaCoords:
    10611094                    DeltaCoords[at] = [0,0,0]
    1062                 for j,(k,n) in enumerate(zip(row,normlist)):
    1063                     if k == 0: continue
    1064                     if l: l += ' + '
    1065                     l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
    1066                     s += n * modeVarDelta[modelist[j]] * k
    1067                 print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
    10681095                if d == 'dx':
    10691096                    DeltaCoords[at][0] = s
    … …  
    10721099                elif d == 'dz':
    10731100                    DeltaCoords[at][2] = s
    1074                 else:
    1075                     print('unexpected',d)
     1101                #else:
     1102                #    print('unexpected',d)
    10761103
    10771104            print( 70*'=')
    1078             print('\nCoordinate checks')
    1079             print('\nComputed')
     1105            print('Coordinate checks')
     1106            print("\nxyz's Computed from ISO mode values, as above")
    10801107            for at in sorted(DeltaCoords):
    10811108                s = at
    … …  
    10861113
    10871114            # determine the coordinate delta values from deviations from the parent structure
    1088             print('\nFrom CIF')
     1115            print("\nxyz Values read directly from CIF")
    10891116            for atmline in self.Phase['Atoms']:
    10901117                lbl = atmline[0]
    … …  
    10921119                print( lbl,x,y,z)
    10931120
    1094             print( 70*'=')
    1095             print('\nGenerated constraints')
     1121            print('\n' + 70*'=')
     1122            print("G2 short name ==> ISODISTORT full name",
     1123                      " (from IsoModeList and G2ModeList)")
     1124            for mode,G2mode in zip(self.Phase['ISODISTORT']['IsoModeList'],
     1125                                   self.Phase['ISODISTORT']['G2ModeList']):
     1126                print("  "+str(G2mode),' ==>', mode)
     1127            print('\nConstraint help dict info')
    10961128            for i in self.Constraints:
    10971129                if type(i) is dict:
    1098                     print('\nconstraint help dict')
    10991130                    for key in i:
    11001131                        print('\t',key,':',i[key])
    1101                 elif i[-1] == 'f':
    1102                     print('\n\t',i[-3],' =')
    1103                     for m,j in i[:-3]:
    1104                         print('\t\t+',m,' * ',j,'   ',repr(j))
    1105                 else:
    1106                     print('  unexpected: ',repr(i))
    1107             print("\nG2name ==> ISODISTORT full name",
    1108                       ".Phase['ISODISTORT']['IsoModeList']",
    1109                       ".Phase['ISODISTORT']['G2ModeList']")
    1110             for mode,G2mode in zip(modelist,modeVarList):
    1111                 print("  ?::"+str(G2mode),' ==>', mode)
     1132            print( 70*'=')
    11121133
    11131134        #======================================================================
    11141135        # debug: show occupancy mode var to mode relations
    1115         if 'G2OccVarList' in self.Phase['ISODISTORT']:
     1136        if debug and 'G2OccVarList' in self.Phase['ISODISTORT']:
    11161137            # coordinate items
    11171138            #occVarLbl = self.Phase['ISODISTORT']['OccVarList']
Note: See TracChangeset for help on using the changeset viewer.