Changeset 4399 for trunk


Ignore:
Timestamp:
Apr 12, 2020 8:46:57 PM (3 years ago)
Author:
toby
Message:

major refactoring & corrections to ISODISTORT mode import and computations

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIO.py

    r4374 r4399  
    22292229def postURL(URL,postdict):
    22302230    '''Posts a set of values as from a web form. If access fails to an https
    2231     site the access is retried with http
     2231    site the access is retried with http.
     2232
    22322233    :param str URL: the URL to post; typically something
    22332234       like 'http://www.../dir/page?'

  • trunk/GSASIIconstrGUI.py

    r4379 r4399  
    1717from __future__ import division, print_function
    1818import sys
     19import copy
     20import os.path
    1921import wx
    2022import wx.grid as wg
     23import wx.lib.scrolledpanel as wxscroll
    2124import random as ran
    2225import numpy as np
    2326import numpy.ma as ma
    2427import numpy.linalg as nl
    25 import os.path
    2628import GSASIIpath
    2729GSASIIpath.SetVersionNumber("$Revision$")
    … …  
    3739import GSASIIobj as G2obj
    3840import GSASIIspc as G2spc
     41import GSASIIpy3 as G2py3
    3942VERY_LIGHT_GREY = wx.Colour(235,235,235)
    4043
    … …  
    822825                        varMean = G2obj.fmtVarDescr(var)
    823826                        helptext += "\n" + var + " ("+ varMean + ")"
     827                    # Add ISODISTORT help items
    824828                    if '_Explain' in data:
    825                         if data['_Explain'].get(item[-3]):
     829                        # this ignores the phase number. TODO: refactor that in
     830                        hlptxt = None
     831                        try:
     832                            hlptxt = data['_Explain'].get(item[-3])
     833                        except TypeError:
     834                            # note that phase RanId is item[-3].phase
     835                            hlptxt = data['_Explain'].get(str(item[-3].phase)+item[-3].name)
     836                        if hlptxt:
    826837                            helptext += '\n\n'
    827                             helptext += data['_Explain'][item[-3]]
     838                            helptext += hlptxt
    828839                    # typeString = 'NEWVAR'
    829840                    # if item[-3]:
    … …  
    832843                    #     eqString[-1] += ' = New Variable'
    833844                    if item[-3]:
    834                         typeString = item[-3] + ' = '
     845                        typeString = str(item[-3]) + ' = '
    835846                    else:
    836847                        typeString = 'New Variable = '
    … …  
    937948            constType = 'New Variable'
    938949            if data[name][Id][-3]:
    939                 varname = data[name][Id][-3]
     950                varname = str(data[name][Id][-3])
    940951            else:
    941952                varname = ""
    … …  
    945956        elif data[name][Id][-1] == 'c':
    946957            items = data[name][Id][:-3]+[
    947                 [data[name][Id][-3],'fixed value =']]
     958                [str(data[name][Id][-3]),'fixed value =']]
    948959            constType = 'Constraint'
    949960            lbl = 'Edit value for each term in constant constraint sum'
    … …  
    957968            return
    958969        try:
    959             prev = data[name][Id][:]
     970            prev = copy.deepcopy(data[name][Id])
    960971            if dlg.ShowModal() == wx.ID_OK:
    961972                result = dlg.GetData()
    … …  
    967978                    data[name][Id][-3] = str(result[-1][0])
    968979                elif data[name][Id][-1] == 'f':
    969                     # process the variable name to put in global form (::var)
    970                     varname = str(dlg.newvar[0]).strip().replace(' ','_')
    971                     if varname.startswith('::'):
    972                         varname = varname[2:]
    973                     varname = varname.replace(':',';')
    974                     if varname:
    975                         data[name][Id][-3] = varname
    976                     else:
    977                         data[name][Id][-3] = ''
    978980                    data[name][Id][-2] = dlg.newvar[1]
     981                    if type(data[name][Id][-3]) is str:
     982                        # process the variable name to put in global form (::var)
     983                        varname = str(dlg.newvar[0]).strip().replace(' ','_')
     984                        if varname.startswith('::'):
     985                            varname = varname[2:]
     986                        varname = varname.replace(':',';')
     987                        if varname:
     988                            data[name][Id][-3] = varname
     989                        else:
     990                            data[name][Id][-3] = ''
    979991                if not CheckChangedConstraint():
    980992                    data[name][Id] = prev
     993            else:
     994                data[name][Id] = prev
    981995        except:
    982996            import traceback
    … …  
    10861100        return
    10871101    G2obj.IndexAllIds(Histograms,Phases)
     1102    for p in Phases:
     1103        if 'ISODISTORT' in Phases[p]:
     1104            G2frame.dataWindow.ConstraintEdit.Enable(G2G.wxID_SHOWISO,True)
     1105            break
    10881106    ##################################################################################
    10891107    # patch: convert old-style (str) variables in constraints to G2VarObj objects
    … …  
    11771195    G2frame.Bind(wx.EVT_MENU, OnAddAtomEquiv, id=G2G.wxID_EQUIVALANCEATOMS)
    11781196#    G2frame.Bind(wx.EVT_MENU, OnAddRiding, id=G2G.wxID_ADDRIDING)
     1197    def OnShowISODISTORT(event):
     1198        ShowIsoDistortCalc(G2frame)
     1199    G2frame.Bind(wx.EVT_MENU, OnShowISODISTORT, id=G2G.wxID_SHOWISO)
    11791200    # tab commands
    11801201    for id in (G2G.wxID_CONSPHASE,
    … …  
    23782399    G2frame.rbBook.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, OnPageChanged)
    23792400    wx.CallAfter(OnPageChanged,None)
     2401
     2402def ShowIsoDistortCalc(G2frame,phase=None):
     2403    '''Compute the ISODISTORT mode values from the current coordinates.
     2404    Called in response to the (Phase/Atoms tab) AtomCompute or
     2405    Constraints/Edit Constr. "Show ISODISTORT modes" menu item, which
     2406    should be enabled only when Phase['ISODISTORT'] is defined.
     2407    '''
     2408    def _onClose(event):
     2409        dlg.EndModal(wx.ID_CANCEL)
     2410    def fmtHelp(item,fullname):
     2411        helptext = "A new variable"
     2412        if item[-3]:
     2413            helptext += " named "+str(item[-3])
     2414        helptext += " is a linear combination of the following parameters:\n"
     2415        first = True
     2416        for term in item[:-3]:
     2417            line = ''
     2418            var = str(term[1])
     2419            m = term[0]
     2420            if first:
     2421                first = False
     2422                line += ' = '
     2423            else:
     2424                if m >= 0:
     2425                    line += ' + '
     2426                else:
     2427                    line += ' - '
     2428                m = abs(m)
     2429            line += '%.3f*%s '%(m,var)
     2430            varMean = G2obj.fmtVarDescr(var)
     2431            helptext += "\n" + line + " ("+ varMean + ")"
     2432        helptext += '\n\nISODISTORT full name: '+str(fullname)
     2433        return helptext
     2434
     2435    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() # init for constraint
     2436    isophases = [p for p in Phases if 'ISODISTORT' in Phases[p]]
     2437   
     2438    if not isophases:
     2439        G2frame.ErrorDialog('no ISODISTORT phases',
     2440                            'Unexpected error: no ISODISTORT phases')
     2441        return
     2442    if phase and phase not in isophases:
     2443        G2frame.ErrorDialog('not ISODISTORT phase',
     2444                            'Unexpected error: selected phase is not an ISODISTORT phase')
     2445        print('Unexpected error: selected phase is not an ISODISTORT phase',
     2446                  phase,isophases)
     2447    elif not phase and len(isophases) == 1:
     2448        phase = isophases[0]
     2449    elif not phase:
     2450        dlg = wx.SingleChoiceDialog(G2frame,'Select phase from ISODISTORT phases',
     2451                                        'Select Phase',isophases)
     2452        if dlg.ShowModal() == wx.ID_OK:
     2453            sel = dlg.GetSelection()
     2454            phase = isophases[sel]
     2455        else:
     2456            return
     2457    # if len(data.get('Histograms',[])) == 0:
     2458    #     G2frame.ErrorDialog(
     2459    #         'No data',
     2460    #         'Sorry, this computation requires that a histogram first be added to the phase'
     2461    #         )
     2462    #     return
     2463   
     2464    covdata = G2frame.GPXtree.GetItemPyData(
     2465        G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Covariance'))
     2466    # make a lookup table for named NewVar Phase constraints
     2467    sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints')
     2468    Constraints = G2frame.GPXtree.GetItemPyData(sub)
     2469    constDict = {}
     2470    for c in Constraints['Phase']:
     2471        if c[-1] != 'f' or not c[-3]: continue
     2472        constDict[str(c[-3])] = c
     2473
     2474    parmDict,varyList = G2frame.MakeLSParmDict()
     2475
     2476    dlg = wx.Dialog(G2frame,wx.ID_ANY,'ISODISTORT mode values',#size=(630,400),
     2477                       style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
     2478    mainSizer = wx.BoxSizer(wx.VERTICAL)
     2479    if 'ISODISTORT' not in Phases[phase]:
     2480        G2frame.ErrorDialog('not ISODISTORT phase',
     2481                            'Unexpected error: selected phase is not an ISODISTORT phase')
     2482        return
     2483    data = Phases[phase]
     2484    ISO = data['ISODISTORT']
     2485    mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
     2486                                'ISODISTORT mode computation for cordinates in phase '+
     2487                                str(data['General'].get('Name'))))
     2488    aSizer = wx.BoxSizer(wx.HORIZONTAL)
     2489    panel1 = wxscroll.ScrolledPanel(
     2490        dlg, wx.ID_ANY,#size=(100,200),
     2491        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
     2492    subSizer1 = wx.FlexGridSizer(cols=2,hgap=5,vgap=2)
     2493    panel2 = wxscroll.ScrolledPanel(
     2494        dlg, wx.ID_ANY,#size=(100,200),
     2495        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
     2496    subSizer2 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
     2497    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'Parameter name  '))
     2498    subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
     2499    subSizer2.Add((-1,-1))
     2500    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'Mode name  '))
     2501    subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
     2502    # ISODISTORT displacive modes
     2503    if 'G2VarList' in ISO:
     2504        deltaList = []
     2505        for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
     2506            dvar = gv.varname()
     2507            var = dvar.replace('::dA','::A')
     2508            albl = Ilbl[:Ilbl.rfind('_')]
     2509            v = Ilbl[Ilbl.rfind('_')+1:]
     2510            pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
     2511            if var in parmDict:
     2512                cval = parmDict[var][0]
     2513            else:
     2514                dlg.EndModal(wx.ID_CANCEL)
     2515                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
     2516                return
     2517            deltaList.append(cval-pval)
     2518        modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
     2519        for lbl,xyz,var,val,norm,G2mode in zip(
     2520                ISO['IsoVarList'],deltaList,
     2521                ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
     2522            #GSASIIpath.IPyBreak()
     2523            if str(G2mode) in constDict:
     2524                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
     2525                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
     2526            else:
     2527                subSizer2.Add((-1,-1))
     2528            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
     2529            try:
     2530                value = G2py3.FormatSigFigs(xyz)
     2531            except TypeError:
     2532                value = str(xyz)           
     2533            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
     2534            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
     2535            try:
     2536                value = G2py3.FormatSigFigs(val/norm)
     2537                if 'varyList' in covdata:
     2538                    if str(G2mode) in covdata['varyList']:
     2539                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
     2540                        value = G2mth.ValEsd(val/norm,sig/norm)
     2541            except TypeError:
     2542                value = '?'
     2543            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
     2544            #GSASIIpath.IPyBreak()
     2545    # ISODISTORT occupancy modes
     2546    if 'G2OccVarList' in ISO:
     2547        deltaList = []
     2548        for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
     2549            var = gv.varname()
     2550            albl = Ilbl[:Ilbl.rfind('_')]
     2551            pval = ISO['BaseOcc'][albl]
     2552            if var in parmDict:
     2553                cval = parmDict[var][0]
     2554            else:
     2555                dlg.EndModal(wx.ID_CANCEL)
     2556                G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
     2557                return
     2558            deltaList.append(cval-pval)
     2559        modeVals = np.inner(ISO['Var2OccMatrix'],deltaList)
     2560        for lbl,delocc,var,val,norm,G2mode in zip(
     2561                ISO['OccVarList'],deltaList,
     2562                ISO['OccModeList'],modeVals,ISO['OccNormList'],ISO['G2OccModeList']):
     2563            if str(G2mode) in constDict:
     2564                ch = G2G.HelpButton(panel2,fmtHelp(constDict[str(G2mode)],var))
     2565                subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|wx.ALIGN_CENTER_VERTICAL|wx.ALIGN_CENTER,1)
     2566            else:
     2567                subSizer2.Add((-1,-1))
     2568            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
     2569            try:
     2570                value = G2py3.FormatSigFigs(delocc)
     2571            except TypeError:
     2572                value = str(delocc)
     2573            subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
     2574            #subSizer.Add((10,-1))
     2575            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
     2576            try:
     2577                value = G2py3.FormatSigFigs(val/norm)
     2578                if 'varyList' in covdata:
     2579                    if str(G2mode) in covdata['varyList']:
     2580                        sig = covdata['sig'][covdata['varyList'].index(str(G2mode))]
     2581                        value = G2mth.ValEsd(val/norm,sig/norm)
     2582            except TypeError:
     2583                value = '?'
     2584            subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
     2585
     2586    # finish up ScrolledPanel
     2587    panel1.SetSizer(subSizer1)
     2588    panel2.SetSizer(subSizer2)
     2589    panel1.SetAutoLayout(1)
     2590    panel1.SetupScrolling()
     2591    panel2.SetAutoLayout(1)
     2592    panel2.SetupScrolling()
     2593    # Allow window to be enlarged but not made smaller
     2594    dlg.SetSizer(mainSizer)
     2595    w1,l1 = subSizer1.GetSize()
     2596    w2,l2 = subSizer2.GetSize()
     2597    panel1.SetMinSize((w1+10,200))
     2598    panel2.SetMinSize((w2+20,200))
     2599    aSizer.Add(panel1,1, wx.ALL|wx.EXPAND,1)
     2600    aSizer.Add(panel2,2, wx.ALL|wx.EXPAND,1)
     2601    mainSizer.Add(aSizer,1, wx.ALL|wx.EXPAND,1)
     2602
     2603    # make OK button
     2604    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
     2605    btn = wx.Button(dlg, wx.ID_CLOSE)
     2606    btn.Bind(wx.EVT_BUTTON,_onClose)
     2607    btnsizer.Add(btn)
     2608    mainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER|wx.ALL, 5)
     2609
     2610    mainSizer.Fit(dlg)
     2611    dlg.SetMinSize(dlg.GetSize())
     2612    dlg.ShowModal()
     2613    dlg.Destroy()

  • trunk/GSASIIdataGUI.py

    r4396 r4399  
    10411041                sub = GetGPXtreeItemId(self,self.root,'Constraints') # was created in CheckNotebook if needed
    10421042                Constraints = self.GPXtree.GetItemPyData(sub)               
    1043                 # TODO: make sure that NEWVAR names are unique here?
    10441043                for i in rd.Constraints:
    10451044                    if type(i) is dict:
    1046                         #for j in i: print j,' --> ',i[j]
    10471045                        if '_Explain' not in Constraints: Constraints['_Explain'] = {}
    10481046                        Constraints['_Explain'].update(i)
    … …  
    53095307#            help='Add H atom riding constraints between atom parameter values')
    53105308#        self.ConstraintEdit.Enable(G2G.wxID_ADDRIDING,False)
     5309        G2G.Define_wxId('wxID_SHOWISO')
     5310        self.ConstraintEdit.Append(G2G.wxID_SHOWISO,'Show ISODISTORT modes',
     5311                'Show ISODISTORT mode values for all phases')
     5312        self.ConstraintEdit.Enable(G2G.wxID_SHOWISO,False)
     5313        # for DEBUG only
     5314        #def OnShowISODISTORT(event):
     5315        #    import imp
     5316        #    imp.reload(G2cnstG)  # for testing changes to GSASIIconstrGUI
     5317        #    G2cnstG.ShowIsoDistortCalc(G2frame)
     5318        #G2frame.Bind(wx.EVT_MENU, OnShowISODISTORT, id=G2G.wxID_SHOWISO)
     5319        # end DEBUG
     5320       
    53115321        self.PostfillDataMenu()
    53125322

  • trunk/GSASIImapvars.py

    r4379 r4399  
    162162A "Hold" constraint is stored as described for type "h" in the
    163163:ref:`constraint definitions table <Constraint_definitions_table>`.
     164
     165.. _Constraints_processing:
    164166
    165167Constraint Processing
    … …  
    14451447            s3 = ' sig   :'
    14461448        s1 += '%15s' % (name)
    1447         s2 += '%15.4f' % (val)
     1449        s2 += '%15.5f' % (val)
    14481450        if esd is None:
    14491451            s3 += '%15s' % ('n/a')
    1450         else:   
    1451             s3 += '%15.4f' % (esd)
     1452        else:
     1453            s3 += '%15.5f' % (esd)
    14521454
    14531455def ComputeDepESD(covMatrix,varyList,parmDict):

  • trunk/GSASIIobj.py

    r4370 r4399  
    10001000letter or string flag value (such as I or A for iso/anisotropic).
    10011001
     1002ISODISTORT implementation
     1003------------------------------
     1004
     1005CIFs prepared with the ISODISTORT web site
     1006https://stokes.byu.edu/iso/isodistort_version5.6.1/isodistort.php
     1007[B. J. Campbell, H. T. Stokes, D. E. Tanner, and D. M. Hatch, "ISODISPLACE: An Internet Tool for Exploring Structural Distortions." J. Appl. Cryst. 39, 607-614 (2006).]
     1008can be read into GSAS-II using import CIF. This will cause constraints to be established for structural distortion modes read from the CIF.
     1009At present, of the five types of modes  only displacive(``_iso_displacivemode``...) and occupancy (``_iso_occupancymode``...) are processed. Not yet processed: ``_iso_magneticmode``..., ``_iso_rotationalmode``... & ``_iso_strainmode``...
     1010
     1011The CIF importer :mod:`G2phase_CIF` implements class :class:`G2phase_CIF.CIFPhaseReader` which offers two methods associated with ISODISTORT (ID) input. Method :meth:`G2phase_CIF.CIFPhaseReader.ISODISTORT_test`
     1012checks to see if a CIF block contains the loops with ``_iso_displacivemode_label`` or  ``_iso_occupancymode_label`` items.
     1013If so, method :meth:`G2phase_CIF.CIFPhaseReader.ISODISTORT_proc` is called to read and interpret them.
     1014The results are placed into the reader object's ``.Phase`` class variable as a dict item with key ``'ISODISTORT'``.
     1015
     1016Note that each mode ID has a long label with a name such as  Pm-3m[1/2,1/2,1/2]R5+(a,a,0)[La:b:dsp]T1u(a). Function :func:`G2phase_CIF.ISODISTORT_shortLbl` is used to create a short name for this,
     1017such as R5_T1u(a) which is made unique by addition of _n if the short name is duplicated.
     1018As each mode is processed, a constraint corresponding to that mode is
     1019created and is added to list in the reader object's ``.Constraints`` class variable. Items placed into that list can either be a list, which corresponds to a function (new var) type :ref:`constraint definition <Constraints_table>` entry, or an item can be a dict, which provides help information for each constraint.
     1020
     1021------------------------------
     1022Displacive modes
     1023------------------------------
     1024
     1025The coordinate variables, as named by ISODISTORT, are placed in
     1026``.Phase['ISODISTORT']['IsoVarList']`` and the corresponding :class:`GSASIIobj.G2VarObj` objects for each are placed in ``.Phase['ISODISTORT']['G2VarList']``. The mode variables,
     1027as named by ISODISTORT, are placed in ``.Phase['ISODISTORT']['IsoModeList']`` and the corresponding :class:`GSASIIobj.G2VarObj` objects for each are placed in ``.Phase['ISODISTORT']['G2ModeList']``.
     1028[Use ``str(G2VarObj)`` to get the variable name from the G2VarObj object, but note that the phase number, *n*, for the prefix "*n*::" cannot be determined as the phase number is not yet assigned.]
     1029
     1030Displacive modes are a bit complex in that they relate to delta displacements, relative to an offset value for each coordinate, and because the modes are normalized, while GSAS-II also uses displacements, but these are added to the coordinates after each refinement cycle and the delta values are set to zero. ISODISTORT uses fixed offsets (subtracted from the actual position to obtain the delta values) that are taken from ``_iso_coordinate_formula`` and these are placed in ``.Phase['ISODISTORT']['ParentStructure]`` (keyed by atom label). The normalization factors (which the delta values are divided by) are taken from ``_iso_displacivemodenorm_value`` and are placed in ``.Phase['ISODISTORT']['NormList']`` in the same order as as ``...['IsoModeList']`` and
     1031``...['G2ModeList']``.
     1032
     1033The CIF contains a sparse matrix, from the ``loop_`` containing
     1034``_iso_displacivemodematrix_value`` which provides the equations for determining the mode values from the coordinates, that matrix is placed in ``.Phase['ISODISTORT']['Mode2VarMatrix']``. The matrix is inverted to produce
     1035``.Phase['ISODISTORT']['Var2ModeMatrix']``, which determines how to compute the
     1036mode values from the delta coordinate values. These values are used for the
     1037in :func:`GSASIIconstrGUI.ShowIsoDistortCalc` which shows coordinate and mode
     1038values, the latter with s.u. values.
     1039
     1040
     1041------------------------------
     1042Occupancy modes
     1043------------------------------
     1044
     1045
     1046The delta occupancy variables, as named by ISODISTORT, are placed in
     1047``.Phase['ISODISTORT']['OccVarList']`` and the corresponding :class:`GSASIIobj.G2VarObj` objects for each are placed in ``.Phase['ISODISTORT']['G2OccVarList']``. The mode variables, as named by ISODISTORT, are placed in ``.Phase['ISODISTORT']['OccModeList']`` and the corresponding :class:`GSASIIobj.G2VarObj` objects for each are placed in ``.Phase['ISODISTORT']['G2OccModeList']``.
     1048
     1049Occupancy modes, like Displacive modes, are also refined as delta values.  However, GSAS-II directly refines the fractional occupancies.
     1050Offset values for each atom, are taken from ``_iso_occupancy_formula`` and are placed in
     1051``.Phase['ISODISTORT']['ParentOcc]``.
     1052(Offset values are subtracted from the actual position to obtain the delta values.)
     1053Modes are normalized (where the mode values are divided by the normalization factor)
     1054are taken from ``_iso_occupancymodenorm_value`` and are placed in ``.Phase['ISODISTORT']['OccNormList']`` in the same order as as ``...['OccModeList']`` and
     1055``...['G2OccModeList']``.
     1056
     1057The CIF contains a sparse matrix, from the ``loop_`` containing
     1058``_iso_occupancymodematrix_value``, which provides the equations for determining the mode values from the coordinates. That matrix is placed in ``.Phase['ISODISTORT']['Occ2VarMatrix']``. The matrix is inverted to produce
     1059``.Phase['ISODISTORT']['Var2OccMatrix']``, which determines how to compute the
     1060mode values from the delta coordinate values.
     1061
     1062
     1063------------------------------
     1064Mode Computations
     1065------------------------------
     1066
     1067Constraints are processed after the CIF has been read in
     1068:meth:`GSASIIdataGUI.GSASII.OnImportPhase` or 
     1069:meth:`GSASIIscriptable.G2Project.add_phase` by moving them from the reader
     1070object's ``.Constraints`` class variable to
     1071the Constraints tree entry's ['Phase'] list (for list items defining constraints) or
     1072the Constraints tree entry's ['_Explain'] dict (for dict items defining constraint help information)
     1073
     1074The information in ``.Phase['ISODISTORT']`` is used
     1075in :func:`GSASIIconstrGUI.ShowIsoDistortCalc` which shows coordinate and mode
     1076values, the latter with s.u. values. This can be called from the Constraints and Phase/Atoms tree items.
     1077
     1078Before each refinement, constraints are processed as
     1079:ref:`described elsewhere <Constraints_processing>`. After a refinement
     1080is complete, :func:`GSASIImapvars.PrintIndependentVars` shows the
     1081shifts and s.u.'s on the refined modes, using GSAS-II values, but
     1082:func:`GSASIIstrIO.PrintISOmodes` prints the ISODISTORT modes as computed
     1083in the web site.
     1084
    10021085
    10031086*Classes and routines*

  • trunk/GSASIIphsGUI.py

    r4398 r4399  
    5858import GSASIImath as G2mth
    5959import GSASIIpwd as G2pwd
    60 import GSASIIpy3 as G2py3
    6160import GSASIIobj as G2obj
    6261import GSASIIctrlGUI as G2G
    … …  
    40934092
    40944093    def OnIsoDistortCalc(event):
    4095         '''Compute the ISODISTORT mode values from the current coordinates.
    4096         Called in response to the (Phase/Atoms tab) AtomCompute
    4097         "Compute ISODISTORT mode values" menu item, which should be enabled
    4098         only when Phase['ISODISTORT'] is defined.
    4099         '''
    4100         def _onClose(event):
    4101             dlg.EndModal(wx.ID_CANCEL)
    4102         def fmtHelp(item,fullname):
    4103             helptext = "A new variable"
    4104             if item[-3]:
    4105                 helptext += " named "+str(item[-3])
    4106             helptext += " is a linear combination of the following parameters:\n"
    4107             first = True
    4108             for term in item[:-3]:
    4109                 line = ''
    4110                 var = str(term[1])
    4111                 m = term[0]
    4112                 if first:
    4113                     first = False
    4114                     line += ' = '
    4115                 else:
    4116                     if m >= 0:
    4117                         line += ' + '
    4118                     else:
    4119                         line += ' - '
    4120                     m = abs(m)
    4121                 line += '%.3f*%s '%(m,var)
    4122                 varMean = G2obj.fmtVarDescr(var)
    4123                 helptext += "\n" + line + " ("+ varMean + ")"
    4124             helptext += '\n\nISODISTORT full name: '+str(fullname)
    4125             return helptext
    4126 
    4127         if 'ISODISTORT' not in data:
    4128             raise Exception("Should not happen: 'ISODISTORT' not in data")
    4129         if len(data.get('Histograms',[])) == 0:
    4130             G2frame.ErrorDialog(
    4131                 'No data',
    4132                 'Sorry, this computation requires that a histogram first be added to the phase'
    4133                 )
    4134             return
    4135         Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() # init for constraint
    4136         # make a lookup table for constraints
    4137         sub = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Constraints')
    4138         Constraints = G2frame.GPXtree.GetItemPyData(sub)
    4139         constDict = {}
    4140         for item in Constraints:
    4141             if item.startswith('_'): continue
    4142             for c in Constraints[item]:
    4143                 if c[-1] != 'f' or not c[-3]: continue
    4144                 constDict[c[-3]] = c
    4145 
    4146         ISO = data['ISODISTORT']
    4147         parmDict,varyList = G2frame.MakeLSParmDict()
    4148            
    4149         dlg = wx.Dialog(G2frame,wx.ID_ANY,'ISODISTORT mode values',#size=(630,400),
    4150                            style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
    4151         mainSizer = wx.BoxSizer(wx.VERTICAL)
    4152         mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
    4153                                     'ISODISTORT mode computation for cordinates in phase '+
    4154                                     str(data['General'].get('Name'))))
    4155         aSizer = wx.BoxSizer(wx.HORIZONTAL)
    4156         panel1 = wxscroll.ScrolledPanel(
    4157             dlg, wx.ID_ANY,#size=(100,200),
    4158             style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
    4159         subSizer1 = wx.FlexGridSizer(cols=2,hgap=5,vgap=2)
    4160         panel2 = wxscroll.ScrolledPanel(
    4161             dlg, wx.ID_ANY,#size=(100,200),
    4162             style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
    4163         subSizer2 = wx.FlexGridSizer(cols=3,hgap=5,vgap=2)
    4164         subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,'Parameter name  '))
    4165         subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
    4166         subSizer2.Add((-1,-1))
    4167         subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,'Mode name  '))
    4168         subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,' value'),0,wx.ALIGN_RIGHT)
    4169         if 'G2VarList' in ISO:
    4170             deltaList = []
    4171             for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
    4172                 dvar = gv.varname()
    4173                 var = dvar.replace('::dA','::A')
    4174                 albl = Ilbl[:Ilbl.rfind('_')]
    4175                 v = Ilbl[Ilbl.rfind('_')+1:]
    4176                 pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
    4177                 if var in parmDict:
    4178                     cval = parmDict[var][0]
    4179                 else:
    4180                     dlg.EndModal(wx.ID_CANCEL)
    4181                     G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
    4182                     return
    4183                 deltaList.append(cval-pval)
    4184             modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
    4185             for lbl,xyz,var,val,G2var in zip(ISO['IsoVarList'],deltaList,
    4186                                              ISO['IsoModeList'],modeVals,ISO['G2ModeList']):
    4187                 if G2var in constDict:
    4188                     ch = G2G.HelpButton(panel2,fmtHelp(constDict[G2var],var))
    4189                     subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
    4190                 else:
    4191                     subSizer2.Add((-1,-1))
    4192                 subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
    4193                 try:
    4194                     value = G2py3.FormatSigFigs(xyz)
    4195                 except TypeError:
    4196                     value = str(xyz)           
    4197                 subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
    4198                 subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
    4199                 try:
    4200                     value = G2py3.FormatSigFigs(val)
    4201                 except TypeError:
    4202                     value = str(val)           
    4203                 subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
    4204         if 'G2OccVarList' in ISO:
    4205             deltaList = []
    4206             for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
    4207                 var = gv.varname()
    4208                 albl = Ilbl[:Ilbl.rfind('_')]
    4209                 #v = Ilbl[Ilbl.rfind('_')+1:]
    4210                 pval = ISO['BaseOcc'][albl]
    4211                 if var in parmDict:
    4212                     cval = parmDict[var][0]
    4213                 else:
    4214                     dlg.EndModal(wx.ID_CANCEL)
    4215                     G2frame.ErrorDialog('Atom not found',"No value found for parameter "+str(var))
    4216                     return
    4217                 deltaList.append(cval-pval)
    4218             modeVals = np.inner(ISO['Var2OccMatrix'],deltaList)
    4219             for lbl,xyz,var,val,G2var in zip(ISO['OccVarList'],deltaList,
    4220                                              ISO['OccModeList'],modeVals,ISO['G2OccModeList']):
    4221                 if G2var in constDict:
    4222                     ch = G2G.HelpButton(panel2,fmtHelp(constDict[G2var],var))
    4223                     subSizer2.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
    4224                 else:
    4225                     subSizer2.Add((-1,-1))
    4226                 subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,str(lbl)))
    4227                 try:
    4228                     value = G2py3.FormatSigFigs(xyz)
    4229                 except TypeError:
    4230                     value = str(xyz)           
    4231                 subSizer1.Add(wx.StaticText(panel1,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
    4232                 #subSizer.Add((10,-1))
    4233                 subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,str(var)))
    4234                 try:
    4235                     value = G2py3.FormatSigFigs(val)
    4236                 except TypeError:
    4237                     value = str(val)           
    4238                 subSizer2.Add(wx.StaticText(panel2,wx.ID_ANY,value),0,wx.ALIGN_RIGHT)
    4239 
    4240         # finish up ScrolledPanel
    4241         panel1.SetSizer(subSizer1)
    4242         panel2.SetSizer(subSizer2)
    4243         panel1.SetAutoLayout(1)
    4244         panel1.SetupScrolling()
    4245         panel2.SetAutoLayout(1)
    4246         panel2.SetupScrolling()
    4247         # Allow window to be enlarged but not made smaller
    4248         dlg.SetSizer(mainSizer)
    4249         w1,l1 = subSizer1.GetSize()
    4250         w2,l2 = subSizer2.GetSize()
    4251         panel1.SetMinSize((w1+10,200))
    4252         panel2.SetMinSize((w2+20,200))
    4253         aSizer.Add(panel1,1, wx.ALL|wx.EXPAND,1)
    4254         aSizer.Add(panel2,2, wx.ALL|wx.EXPAND,1)
    4255         mainSizer.Add(aSizer,1, wx.ALL|wx.EXPAND,1)
    4256 
    4257         # make OK button
    4258         btnsizer = wx.BoxSizer(wx.HORIZONTAL)
    4259         btn = wx.Button(dlg, wx.ID_CLOSE)
    4260         btn.Bind(wx.EVT_BUTTON,_onClose)
    4261         btnsizer.Add(btn)
    4262         mainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER|wx.ALL, 5)
    4263 
    4264         mainSizer.Fit(dlg)
    4265         dlg.SetMinSize(dlg.GetSize())
    4266         dlg.ShowModal()
    4267         dlg.Destroy()
    4268        
     4094        G2cnstG.ShowIsoDistortCalc(G2frame,data['General']['Name'])
     4095           
    42694096    def OnReImport(event):
    42704097        generalData = data['General']

  • trunk/GSASIIscriptable.py

    r4334 r4399  
    22402240        self.data['Phases'][phasename] = phasereader.Phase
    22412241
     2242        # process constraints, currently generated only from ISODISTORT CIFs
    22422243        if phasereader.Constraints:
    22432244            Constraints = self.data['Constraints']

  • trunk/GSASIIstrIO.py

    r4371 r4399  
    4040import GSASIImath as G2mth
    4141import GSASIIfiles as G2fil
     42import GSASIIpy3 as G2py3
    4243
    4344sind = lambda x: np.sin(x*np.pi/180.)
    … …  
    211212            if len(D) > 1:
    212213                # add extra dict terms for input variable name and vary flag
    213                 if varname is not None:                   
    214                     if varname.startswith('::'):
    215                         varname = varname[2:].replace(':',';')
     214                if varname is not None:
     215                    varname = str(varname) # in case this is a G2VarObj
     216                    if ':' in varname:
     217                        D['_name'] = varname
    216218                    else:
    217                         varname = varname.replace(':',';')
    218                     D['_name'] = '::' + varname
     219                        D['_name'] = '::' + varname
    219220                D['_vary'] = varyFlag == True # force to bool
    220221                constDict.append(D)
    … …  
    367368    pFile.write(FPstr+'\n')
    368369    pFile.write(FPPstr+'\n')
     370
     371def PrintISOmodes(pFile,Phases,parmDict,sigDict):
     372    '''Prints the values for the ISODISTORT modes into the project's
     373    .lst file after a refinement.
     374    '''
     375    for phase in Phases:
     376        if 'ISODISTORT' not in Phases[phase]: continue
     377        data = Phases[phase]
     378        ISO = data['ISODISTORT']
     379        if 'G2VarList' in ISO:
     380            deltaList = []
     381            for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
     382                dvar = gv.varname()
     383                var = dvar.replace('::dA','::A')
     384                albl = Ilbl[:Ilbl.rfind('_')]
     385                v = Ilbl[Ilbl.rfind('_')+1:]
     386                pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
     387                if var in parmDict:
     388                    cval = parmDict[var]
     389                else:
     390                    print('PrintISOmodes Error: Atom not found',"No value found for parameter "+str(var))
     391                    return
     392                deltaList.append(cval-pval)
     393            modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
    369394           
     395        if 'G2OccVarList' in ISO:
     396            deltaOccList = []
     397            for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
     398                var = gv.varname()
     399                albl = Ilbl[:Ilbl.rfind('_')]
     400                pval = ISO['BaseOcc'][albl]
     401                if var in parmDict:
     402                    cval = parmDict[var]
     403                else:
     404                    print('PrintISOmodes Error: Atom not found',"No value found for parameter "+str(var))
     405                    return
     406                deltaOccList.append(cval-pval)
     407            modeOccVals = np.inner(ISO['Var2OccMatrix'],deltaOccList)
     408           
     409        if 'G2VarList' in ISO:
     410            pFile.write('\n ISODISTORT Displacive Modes for phase {}\n'.format(data['General'].get('Name','')))
     411            l = str(max([len(i) for i in ISO['IsoModeList']])+3)
     412            fmt = '  {:'+l+'}{}'
     413            for var,val,norm,G2mode in zip(
     414                    ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
     415                try:
     416                    value = G2py3.FormatSigFigs(val/norm)
     417                    if str(G2mode) in sigDict:
     418                        value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm)
     419                except TypeError:
     420                    value = '?'
     421                pFile.write(fmt.format(var,value)+'\n')
     422               
     423        if 'G2OccVarList' in ISO:
     424            pFile.write('\n ISODISTORT Occupancy Modes for phase {}\n'.format(data['General'].get('Name','')))
     425            l = str(max([len(i) for i in ISO['OccModeList']])+3)
     426            fmt = '  {:'+l+'}{}'
     427            for var,val,norm,G2mode in zip(
     428                    ISO['OccModeList'],modeOccVals,ISO['OccNormList'],ISO['G2OccModeList'] ):
     429                try:
     430                    value = G2py3.FormatSigFigs(val/norm)
     431                    if str(G2mode) in sigDict:
     432                        value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm)
     433                except TypeError:
     434                    value = '?'
     435                pFile.write(fmt.format(var,value)+'\n')
     436   
    370437def GetPhaseNames(GPXfile):
    371438    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file

  • trunk/GSASIIstrMain.py

    r4370 r4399  
    275275            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
    276276            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
     277            G2stIO.PrintISOmodes(printFile,Phases,parmDict,sigDict)
    277278            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
    278279            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)

  • trunk/imports/G2phase_CIF.py

    r4379 r4399  
    3434GSASIIpath.SetVersionNumber("$Revision$")
    3535import CifFile as cif # PyCifRW from James Hester
    36 #debug = GSASIIpath.GetConfigValue('debug')
    37 debug = False
     36debug = GSASIIpath.GetConfigValue('debug')
     37#debug = False
    3838
    3939class CIFPhaseReader(G2obj.ImportPhase):
    … …  
    648648        '''Process ISODISTORT items to create constraints etc.
    649649        Constraints are generated from information extracted from
    650         loop_'s beginning with _iso_ and are placed into
     650        loops beginning with _iso_ and are placed into
    651651        self.Constraints, which contains a list of
    652652        :ref:`constraints tree items <Constraint_definitions_table>`
    … …  
    660660        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
    661661        'Maps ISODISTORT parm names to GSAS-II names'
     662        # used for all types of modes
     663        self.Constraints = []
     664        explaination = {}
    662665        #----------------------------------------------------------------------
    663666        # read in the ISODISTORT displacement modes
    664667        #----------------------------------------------------------------------
    665         self.Constraints = []
    666         explaination = {}
    667668        if blk.get('_iso_displacivemode_label'):
    668669            modelist = []
    669670            shortmodelist = []
    670             for lbl in blk.get('_iso_displacivemode_label'):
     671            idlist = []
     672            for id,lbl in zip(
     673                blk.get('_iso_displacivemode_ID'),
     674                blk.get('_iso_displacivemode_label')):
     675                idlist.append(int(id))
    671676                modelist.append(lbl)
    672677                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
     678            # just in case the items are not ordered increasing by id, sort them here
     679            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
     680            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
    673681            # read in the coordinate offset variables names and map them to G2 names/objects
    674682            coordVarLbl = []
    675             G2varLbl = []
    676683            G2varObj = []
     684            idlist = []
    677685            error = False
    678             for lbl in blk.get('_iso_deltacoordinate_label'):
     686            for id,lbl in zip(
     687                blk.get('_iso_deltacoordinate_ID'),
     688                blk.get('_iso_deltacoordinate_label') ):
     689                idlist.append(int(id))
    679690                coordVarLbl.append(lbl)
    680691                if '_' in lbl:
    … …  
    696707                    error = True
    697708                    continue
    698                 G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
    699709                G2varObj.append(G2obj.G2VarObj(
    700710                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
    … …  
    702712            if error:
    703713                raise Exception("Error decoding variable labels")
     714            # just in case the items are not ordered increasing by id, sort them here
     715            coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])]
     716            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
    704717
    705718            if len(G2varObj) != len(modelist):
    … …  
    711724            for lbl,exp in zip(
    712725                blk.get('_iso_coordinate_label'),
    713                 blk.get('_iso_coordinate_formula'),
    714                 ):
     726                blk.get('_iso_coordinate_formula') ):
    715727                if '_' in lbl:
    716728                    albl = lbl[:lbl.rfind('_')]
    … …  
    750762            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
    751763            # create the constraints
     764            modeVarList = []
    752765            for i,row in enumerate(displacivemodeInvmatrix):
    753766                constraint = []
    … …  
    755768                    if k == 0: continue
    756769                    constraint.append([k,G2varObj[j]])
    757                 constraint += [shortmodelist[i],False,'f']
     770                modeVar = G2obj.G2VarObj(
     771                    (self.Phase['ranId'],None,shortmodelist[i],None))
     772                modeVarList.append(modeVar)               
     773                constraint += [modeVar,False,'f']
    758774                self.Constraints.append(constraint)
     775            # normilization constants
     776            normlist = []
     777            idlist = []
     778            for id,exp in zip(
     779                blk.get('_iso_displacivemodenorm_ID'),
     780                blk.get('_iso_displacivemodenorm_value'),
     781                ):
     782                idlist.append(int(id))
     783                normlist.append(float(exp))
     784            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
    759785            #----------------------------------------------------------------------
    760786            # save the ISODISTORT info for "mode analysis"
    761787            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
    762788            self.Phase['ISODISTORT'].update({
    763                 'IsoModeList' : modelist,
    764                 'G2ModeList' : shortmodelist,
     789                # coordinate items
    765790                'IsoVarList' : coordVarLbl,
    766791                'G2VarList' : G2varObj,
    767792                'ParentStructure' : ParentCoordinates,
     793                # mode items
     794                'IsoModeList' : modelist,
     795                'G2ModeList' : modeVarList,
     796                'NormList' : normlist,
     797                # transform matrices
    768798                'Var2ModeMatrix' : displacivemodeInvmatrix,
    769799                'Mode2VarMatrix' : displacivemodematrix,
    … …  
    771801            # make explaination dictionary
    772802            for mode,shortmode in zip(modelist,shortmodelist):
    773                 explaination[shortmode] = "ISODISTORT full name "+str(mode)
     803                explaination[str(self.Phase['ranId'])+shortmode] = (
     804                    "ISODISTORT full name "+str(mode))
    774805        #----------------------------------------------------------------------
    775806        # now read in the ISODISTORT occupancy modes
    … …  
    778809            modelist = []
    779810            shortmodelist = []
    780             for lbl in blk.get('_iso_occupancymode_label'):
     811            idlist = []
     812            for id,lbl in zip(
     813                blk.get('_iso_occupancymode_ID'),
     814                blk.get('_iso_occupancymode_label')):
     815                idlist.append(int(id))
    781816                modelist.append(lbl)
    782817                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
     818            # just in case the items are not ordered increasing by id, sort them here
     819            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
     820            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
    783821            # read in the coordinate offset variables names and map them to G2 names/objects
    784822            occVarLbl = []
    785             G2varLbl = []
    786823            G2varObj = []
     824            idlist = []
    787825            error = False
    788             for lbl in blk.get('_iso_deltaoccupancy_label'):
     826            for id,lbl in zip(
     827                blk.get('_iso_deltaoccupancy_ID'),
     828                blk.get('_iso_deltaoccupancy_label') ):
     829                idlist.append(int(id))
    789830                occVarLbl.append(lbl)
    790831                if '_' in lbl:
    … …  
    806847                    error = True
    807848                    continue
    808                 G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
    809849                G2varObj.append(G2obj.G2VarObj(
    810850                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
    … …  
    812852            if error:
    813853                raise Exception("Error decoding variable labels")
     854            # just in case the items are not ordered increasing by id, sort them here
     855            occVarLbl = [i for i,j in sorted(zip(occVarLbl,idlist),key=lambda k:k[1])]
     856            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
    814857
    815858            if len(G2varObj) != len(modelist):
    … …  
    818861
    819862            error = False
    820             ParentCoordinates = {}
     863            ParentOcc = {}
    821864            for lbl,exp in zip(
    822865                blk.get('_iso_occupancy_label'),
    823                 blk.get('_iso_occupancy_formula'),
    824                 ):
     866                blk.get('_iso_occupancy_formula') ):
    825867                if '_' in lbl:
    826868                    albl = lbl[:lbl.rfind('_')]
    … …  
    841883                        error = True
    842884                        continue
    843                     ParentCoordinates[albl] = val
     885                    ParentOcc[albl] = val
    844886            if error:
    845887                raise Exception("Error decoding occupancy labels")
    … …  
    854896            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
    855897            # create the constraints
     898            modeVarList = []
    856899            for i,row in enumerate(occupancymodeInvmatrix):
    857900                constraint = []
    … …  
    859902                    if k == 0: continue
    860903                    constraint.append([k,G2varObj[j]])
    861                 constraint += [shortmodelist[i],False,'f']
     904                modeVar = G2obj.G2VarObj(
     905                    (self.Phase['ranId'],None,shortmodelist[i],None))
     906                modeVarList.append(modeVar)               
     907                constraint += [modeVar,False,'f']
    862908                self.Constraints.append(constraint)
     909            # normilization constants
     910            normlist = []
     911            idlist = []
     912            for id,exp in zip(
     913                blk.get('_iso_occupancymodenorm_ID'),
     914                blk.get('_iso_occupancymodenorm_value'),
     915                ):
     916                idlist.append(int(id))
     917                normlist.append(float(exp))
     918            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
    863919            #----------------------------------------------------------------------
    864920            # save the ISODISTORT info for "mode analysis"
    865921            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
    866922            self.Phase['ISODISTORT'].update({
    867                 'OccModeList' : modelist,
    868                 'G2OccModeList' : shortmodelist,
     923                # coordinate items
    869924                'OccVarList' : occVarLbl,
    870925                'G2OccVarList' : G2varObj,
    871                 'BaseOcc' : ParentCoordinates,
     926                'BaseOcc' : ParentOcc,
     927                # mode items
     928                'OccModeList' : modelist,
     929                'G2OccModeList' : modeVarList,
     930                'OccNormList' : normlist,
     931                # transform matrices
    872932                'Var2OccMatrix' : occupancymodeInvmatrix,
    873933                'Occ2VarMatrix' : occupancymodematrix,
    … …  
    875935            # make explaination dictionary
    876936            for mode,shortmode in zip(modelist,shortmodelist):
    877                 explaination[shortmode] = "ISODISTORT full name "+str(mode)
     937                explaination[str(self.Phase['ranId'])+shortmode] = (
     938                    "ISODISTORT full name "+str(mode))
    878939        #----------------------------------------------------------------------
    879940        # done with read
    … …  
    883944        if not debug: return
    884945
    885         # debug: show the mode var to mode relations
    886         print("shortname ==> ISODISTORT full name")
    887         for mode,shortmode in zip(modelist,shortmodelist):
    888             print(shortmode,' ==>', mode)
    889         print()
    890         for i,row in enumerate(displacivemodeInvmatrix):
    891             l = ''
    892             for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
    893                 if k == 0: continue
    894                 if l: l += ' + '
    895                 #l += lbl+' * '+str(k)
    896                 l += G2varLbl[j]+' * '+str(k)
    897             print( str(i) + ': '+shortmodelist[i]+' = '+l)
    898         print( 70*'=')
    899 
    900         # debug: Get the ISODISTORT offset values
    901         coordVarDelta = {}
    902         for lbl,val in zip(
    903             blk.get('_iso_deltacoordinate_label'),
    904             blk.get('_iso_deltacoordinate_value'),):
    905             coordVarDelta[lbl] = float(val)
    906         modeVarDelta = {}
    907         for lbl,val in zip(
    908             blk.get('_iso_displacivemode_label'),
    909             blk.get('_iso_displacivemode_value'),):
    910             modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
    911 
    912         print( 70*'=')
    913         # compute the mode values from the reported coordinate deltas
    914         for i,row in enumerate(displacivemodeInvmatrix):
    915             l = ''
    916             sl = ''
    917             s = 0.
    918             for lbl,k in zip(coordVarLbl,row):
    919                 if k == 0: continue
    920                 if l: l += ' + '
    921                 l += lbl+' * '+str(k)
    922                 if sl: sl += ' + '
    923                 sl += str(coordVarDelta[lbl])+' * '+str(k)
    924                 s += coordVarDelta[lbl] * k
    925             print('a'+str(i)+' = '+l)
    926             print('\t= '+sl)
    927             print(modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s)
    928             print()
    929 
    930         print( 70*'=')
    931         # compute the coordinate displacements from the reported mode values
    932         for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
    933             l = ''
    934             sl = ''
    935             s = 0.0
    936             for j,k in enumerate(row):
    937                 if k == 0: continue
    938                 if l: l += ' + '
    939                 l += 'a'+str(j+1)+' * '+str(k)
    940                 if sl: sl += ' + '
    941                 sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
    942                 s += modeVarDelta[modelist[j]] * k
    943             print( lbl+' = '+l)
    944             print( '\t= '+sl)
    945             print( lbl,G2varLbl[i],coordVarDelta[lbl],s)
    946             print()
    947 
    948         # determine the coordinate delta values from deviations from the parent structure
    949         for atmline in self.Phase['Atoms']:
    950             lbl = atmline[0]
    951             x,y,z = atmline[3:6]
    952             if lbl not in ParentCoordinates:
     946        # debug: show displacive mode var to mode relations
     947        if 'IsoVarList' in self.Phase['ISODISTORT']:
     948            # coordinate items
     949            #coordVarLbl = self.Phase['ISODISTORT']['IsoVarList']
     950            G2varObj = self.Phase['ISODISTORT']['G2VarList']
     951            #ParentCoordinates = self.Phase['ISODISTORT']['ParentStructure']
     952            # mode items
     953            modelist = self.Phase['ISODISTORT']['IsoModeList']
     954            modeVarList = self.Phase['ISODISTORT']['G2ModeList']
     955            normlist = self.Phase['ISODISTORT']['NormList']
     956            # transform matrices
     957            #displacivemodeInvmatrix = self.Phase['ISODISTORT']['Var2ModeMatrix']
     958            #displacivemodematrix = self.Phase['ISODISTORT']['Mode2VarMatrix']
     959       
     960            print( 70*'=')
     961            print('\nVar2ModeMatrix' ,'IsoVarList' )
     962            for i,row in enumerate(displacivemodeInvmatrix):
     963                l = ''
     964                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
     965                    if k == 0: continue
     966                    if l: l += ' + '
     967                    #l += lbl+' * '+str(k)
     968                    l += str(G2varObj[j])+' * '+str(k)
     969                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
     970
     971            # Get the ISODISTORT offset values
     972            coordVarDelta = {}
     973            for lbl,val in zip(
     974                blk.get('_iso_deltacoordinate_label'),
     975                blk.get('_iso_deltacoordinate_value'),):
     976                coordVarDelta[lbl] = float(val)
     977            modeVarDelta = {}
     978            for lbl,val in zip(
     979                blk.get('_iso_displacivemode_label'),
     980                blk.get('_iso_displacivemode_value'),):
     981                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
     982
     983            print( 70*'=')
     984            print('\nInverse relations using Var2ModeMatrix, NormList, IsoVarList')
     985            # compute the mode values from the reported coordinate deltas
     986            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
     987                l = ''
     988                for lbl,k in zip(coordVarLbl,row):
     989                    if k == 0: continue
     990                    if l: l += ' + '
     991                    l += lbl+' * '+str(k)
     992                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
     993            print('\nCalculation checks\n')
     994            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
     995                #l = ''
     996                sl = ''
     997                s = 0.
     998                for lbl,k in zip(coordVarLbl,row):
     999                    if k == 0: continue
     1000                    #if l: l += ' + '
     1001                    #l += lbl+' * '+str(k)
     1002                    if sl: sl += ' + '
     1003                    sl += str(coordVarDelta[lbl])+' * '+str(k)
     1004                    s += coordVarDelta[lbl] * k
     1005                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
     1006                print(' ?= ',modeVarDelta[modelist[i]])
     1007                print()
     1008
     1009            print( 70*'=')
     1010            print('\nDirect relations using Mode2VarMatrix, NormList, IsoVarList')
     1011            # compute the coordinate displacements from the reported mode values
     1012            DeltaCoords = {}
     1013            for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
     1014                l = ''
     1015                s = 0.0
     1016                at,d = lbl.split('_')
     1017                if at not in DeltaCoords:
     1018                    DeltaCoords[at] = [0,0,0]
     1019                for j,(k,n) in enumerate(zip(row,normlist)):
     1020                    if k == 0: continue
     1021                    if l: l += ' + '
     1022                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
     1023                    s += n * modeVarDelta[modelist[j]] * k
     1024                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
     1025                if d == 'dx':
     1026                    DeltaCoords[at][0] = s
     1027                elif d == 'dy':
     1028                    DeltaCoords[at][1] = s
     1029                elif d == 'dz':
     1030                    DeltaCoords[at][2] = s
     1031                else:
     1032                    print('unexpected',d)
     1033
     1034            print( 70*'=')
     1035            print('\nCoordinate checks')
     1036            print('\nComputed')
     1037            for at in sorted(DeltaCoords):
     1038                s = at
     1039                for i in range(3):
     1040                    s += '  '
     1041                    s += str(ParentCoordinates[at][i]+DeltaCoords[at][i])
     1042                print(s)
     1043
     1044            # determine the coordinate delta values from deviations from the parent structure
     1045            print('\nFrom CIF')
     1046            for atmline in self.Phase['Atoms']:
     1047                lbl = atmline[0]
     1048                x,y,z = atmline[3:6]
    9531049                print( lbl,x,y,z)
    954                 continue
    955             px,py,pz = ParentCoordinates[lbl]
    956             print( lbl,x,y,z,x-px,y-py,z-pz)
    957         print('\nGenerated constraints')
    958         for i in self.Constraints:
    959             if type(i) is dict:
    960                 print('\n  constraint help')
    961                 for key in i:
    962                     print('\t',key,i[key])
    963             elif i[-1] == 'f':
    964                 print('\n\tNewvar',i[-3],' =')
    965                 for m,j in i[:-3]:
    966                     print('\t\t+',m,' * ',j)
    967             else:
    968                 print('  unexpected: ',repr(i))
    969                    
     1050
     1051            print( 70*'=')
     1052            print('\nGenerated constraints')
     1053            for i in self.Constraints:
     1054                if type(i) is dict:
     1055                    print('\nconstraint help dict')
     1056                    for key in i:
     1057                        print('\t',key,':',i[key])
     1058                elif i[-1] == 'f':
     1059                    print('\n\t',i[-3],' =')
     1060                    for m,j in i[:-3]:
     1061                        print('\t\t+',m,' * ',j,'   ',repr(j))
     1062                else:
     1063                    print('  unexpected: ',repr(i))
     1064            print("\nG2name ==> ISODISTORT full name",
     1065                      ".Phase['ISODISTORT']['IsoModeList']",
     1066                      ".Phase['ISODISTORT']['G2ModeList']")
     1067            for mode,G2mode in zip(modelist,modeVarList):
     1068                print("  ?::"+str(G2mode),' ==>', mode)
     1069
     1070        #======================================================================
     1071        # debug: show occupancy mode var to mode relations
     1072        if 'G2OccVarList' in self.Phase['ISODISTORT']:
     1073            # coordinate items
     1074            #occVarLbl = self.Phase['ISODISTORT']['OccVarList']
     1075            G2varObj = self.Phase['ISODISTORT']['G2OccVarList']
     1076            #ParentOcc = self.Phase['ISODISTORT']['BaseOcc']
     1077            # mode items
     1078            modelist = self.Phase['ISODISTORT']['OccModeList']
     1079            modeVarList = self.Phase['ISODISTORT']['G2OccModeList']
     1080            normlist = self.Phase['ISODISTORT']['OccNormList']
     1081            # transform matrices
     1082            #occupancymodeInvmatrix = self.Phase['ISODISTORT']['Var2OccMatrix']
     1083            #occupancymodematrix = self.Phase['ISODISTORT']['Occ2VarMatrix']
     1084           
     1085            print( 70*'=')
     1086            print('\nVar2OccMatrix' ,'OccVarList' )
     1087            for i,row in enumerate(occupancymodeInvmatrix):
     1088                l = ''
     1089                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
     1090                    if k == 0: continue
     1091                    if l: l += ' + '
     1092                    #l += lbl+' * '+str(k)
     1093                    l += str(G2varObj[j])+' * '+str(k)
     1094                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
     1095
     1096            # Get the ISODISTORT offset values
     1097            occVarDelta = {}
     1098            for lbl,val in zip(
     1099                blk.get('_iso_deltaoccupancy_label'),
     1100                blk.get('_iso_deltaoccupancy_value'),):
     1101                occVarDelta[lbl] = float(val)
     1102            modeVarDelta = {}
     1103            for lbl,val in zip(
     1104                blk.get('_iso_occupancymode_label'),
     1105                blk.get('_iso_occupancymode_value'),):
     1106                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
     1107
     1108            print( 70*'=')
     1109            print('\nInverse relations using Var2OccModeMatrix, OccNormList, OccVarList')
     1110            # compute the mode values from the reported coordinate deltas
     1111            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
     1112                l = ''
     1113                for lbl,k in zip(occVarLbl,row):
     1114                    if k == 0: continue
     1115                    if l: l += ' + '
     1116                    l += lbl+' * '+str(k)
     1117                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
     1118            print('\nCalculation checks\n')
     1119            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
     1120                #l = ''
     1121                sl = ''
     1122                s = 0.
     1123                for lbl,k in zip(occVarLbl,row):
     1124                    if k == 0: continue
     1125                    #if l: l += ' + '
     1126                    #l += lbl+' * '+str(k)
     1127                    if sl: sl += ' + '
     1128                    sl += str(occVarDelta[lbl])+' * '+str(k)
     1129                    s += occVarDelta[lbl] * k
     1130                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
     1131                print(' ?= ',modeVarDelta[modelist[i]])
     1132                print()
     1133
     1134            print( 70*'=')
     1135            print('\nDirect relations using Occ2VarMatrix, OccNormList, OccVarList')
     1136            # compute the coordinate displacements from the reported mode values
     1137            Occ = {}
     1138            for i,lbl,row in zip(range(len(occVarLbl)),occVarLbl,occupancymodematrix):
     1139                l = ''
     1140                s = 0.0
     1141                for j,(k,n) in enumerate(zip(row,normlist)):
     1142                    if k == 0: continue
     1143                    if l: l += ' + '
     1144                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
     1145                    s += n * modeVarDelta[modelist[j]] * k
     1146                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
     1147                j = lbl.split('_')[0]
     1148                Occ[j] = ParentOcc[j]+s
     1149               
     1150            # determine the coordinate delta values from deviations from the parent structure
     1151            print('\nOccupancy from CIF vs computed')
     1152            for atmline in self.Phase['Atoms']:
     1153                lbl = atmline[0]
     1154                print( lbl,atmline[6],Occ[lbl])
     1155
     1156            print( 70*'=')
     1157            print('\nGenerated constraints')
     1158            for i in self.Constraints:
     1159                if type(i) is dict:
     1160                    print('\nconstraint help dict')
     1161                    for key in i:
     1162                        print('\t',key,':',i[key])
     1163                elif i[-1] == 'f':
     1164                    print('\n\t',i[-3],' =')
     1165                    for m,j in i[:-3]:
     1166                        print('\t\t+',m,' * ',j,'   ',repr(j))
     1167                else:
     1168                    print('  unexpected: ',repr(i))
     1169            print("\nG2name ==> ISODISTORT full name",
     1170                      ".Phase['ISODISTORT']['OccModeList']",
     1171                      ".Phase['ISODISTORT']['G2OccModeList']")
     1172            for mode,G2mode in zip(modelist,modeVarList):
     1173                print("  ?::"+str(G2mode),' ==>', mode)
     1174               
    9701175def ISODISTORT_shortLbl(lbl,shortmodelist):
    9711176    '''Shorten model labels and remove special characters
Note: See TracChangeset for help on using the changeset viewer.