implement mag constraints w/CIF import from ISO; remove duplicate constraints; improve constraint error & warning msgs

Author: briantoby

GSASII/GSASIIconstrGUI.py

@@ -1916,14 +1916,16 @@ def SetUniqAj(pId,iA,SGData):
                                 G2lat.cell2A(oldPhase['General']['Cell'][1:7]),
                                 oldPhase['General']['SGData'],
                                 newPhase['General']['SGData'],debug=False)
+    detTrans = np.abs(nl.det(Trans))  # volume ratio
+    #Vratio = newPhase['General']['Cell'][7]/oldPhase['General']['Cell'][7]
 
     # create constraints on HAP Scale, Isotropic size and mustrain
     # (anisotropic constraints would be more complex)
     UseList = newPhase['Histograms']
     Histograms, Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
     for hId,hist in enumerate(UseList):
         hRanId = Histograms[hist]['ranId']
-        constrList = G2lat.GenHAPConstraints(Trans,oRanId,nRanId,hRanId)
+        constrList = G2lat.GenHAPConstraints(detTrans,oRanId,nRanId,hRanId)
         constraints['HAP'] += constrList
 
 #### Rigid bodies #############################################################

GSASII/GSASIIdataGUI.py

@@ -17,6 +17,7 @@
 import re
 import numpy as np
 import numpy.ma as ma
+import numpy.linalg as nl
 import matplotlib as mpl
 try:
     import OpenGL as ogl
@@ -1364,6 +1365,52 @@ def OnImportPhase(self,event):
                         Constraints['_Explain'].update(i)
                     else:
                         Constraints['Phase'].append(i)
+            # make ISODISTORT magnetic phase constraints here
+            Phases = self.GetPhaseData()
+            prevPhases = list(Phases.keys())
+            parentPhase = None
+            Vratio = None
+            if ('ISODISTORT' in rd.Phase and
+                    'XformInfo' in rd.Phase['ISODISTORT'] and 
+                    PhaseName in prevPhases and len(prevPhases) > 1):
+                Trans = rd.Phase['ISODISTORT']['XformInfo'].get('Trans',np.eye(3))
+                Vec = rd.Phase['ISODISTORT']['XformInfo'].get('offset',[0,0,0])
+                prevPhases.pop(prevPhases.index(PhaseName))
+                if len(prevPhases) == 1:
+                    parentPhase = prevPhases[0]
+                else:
+                    dlg = G2G.G2SingleChoiceDialog(self,
+                                f'Select parent phase to {PhaseName}',
+                                'Select phase',prevPhases)
+                    dlg.CenterOnParent()
+                    if dlg.ShowModal() == wx.ID_OK:
+                        sel = dlg.GetSelection()
+                        parentPhase = prevPhases[sel]
+                        dlg.Destroy()
+                    else:
+                        dlg.Destroy()
+                if parentPhase:
+                    oldPhase = Phases[parentPhase]
+                    newPhase = Phases[PhaseName]
+                    oRanId = oldPhase['ranId']
+                    nRanId = newPhase['ranId']
+                    Vratio = np.abs(nl.det(Trans))  # volume ratio
+                    #Vratio = newPhase['General']['Cell'][7]/oldPhase['General']['Cell'][7]
+                    debug = True
+                    # create constraints relating two unit cells
+                    Constraints['Phase'] += G2lat.GenCellConstraints(Trans,oRanId,nRanId,
+                            G2lat.cell2A(oldPhase['General']['Cell'][1:7]),
+                            oldPhase['General']['SGData'],
+                            newPhase['General']['SGData'],debug)
+                    # create constraints relating the atoms in newPhase to their positions in the chemical cell
+                    constr,message = G2lat.MatchGenAtomConstraints(
+                        oldPhase,newPhase,Trans,[0,0,0],Vec)
+                    if message:
+                        wx.MessageDialog(self,message,
+                            'Constraint Gen. Problem',style=wx.ICON_INFORMATION).ShowModal()
+
+                    Constraints['Phase'] += constr
+
         if not newPhaseList: return # somehow, no new phases
         # get a list of existing histograms
         PWDRlist = []
@@ -1434,7 +1481,12 @@ def OnImportPhase(self,event):
                     data['Histograms'][histoName] = G2mth.SetDefaultDData(Inst['Type'][0],histoName,NShkl=NShkl,NDij=NDij)
                 else:
                     raise Exception('Unexpected histogram '+histoName)
+
         Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree() # reindex
+        if Vratio: # make HAP constraints for ISODISTORT mag phase
+            for hist in data['Histograms']:
+                hRanId = Histograms[hist]['ranId']
+                Constraints['HAP'] += G2lat.GenHAPConstraints(Vratio,oRanId,nRanId,hRanId)
         wx.EndBusyCursor()
         self.EnableRefineCommand()
 

GSASII/GSASIIlattice.py

@@ -878,25 +878,24 @@ def GenCellConstraints(Trans,oRanId,nRanId,origA,oSGLaue,nSGLaue,debug=False):
         print('xfm cell',fmt4(A2cell(Anew)))
     return constrList
 
-def GenHAPConstraints(Trans,oRanId,nRanId,hRanId):
-    '''Generate the HAP constraints between two phase constants 
-    for a phase transformed by matrix Trans.
+def GenHAPConstraints(Vratio,oRanId,nRanId,hRanId):
+    '''Generate the HAP constraints between two phases: the original 
+    chemical (nuclear) phase and the new, magnetic phase. 
 
-    :param np.array Trans: a 3x3 direct cell transformation matrix where,
-       Trans = np.array([ [2/3, 4/3, 1/3], [-1, 0, 0], [-1/3, -2/3, 1/3] ])
-       (for a' = 2/3a + 4/3b + 1/3c; b' = -a; c' = -1/3, -2/3, 1/3)
+    :param float Vratio: the ratio of the volume of the magnetic phase
+      divided by the original phase. This should be 1 or >1 [often 2 
+      or sqrt(2) etc.]
     :param int oRanId: Random id for the original phase
-    :param int nRanId: Random id for the transformed phase to be constrained from
-      original phase
+    :param int nRanId: Random id for the transformed (magnetic) phase to 
+      be constrained to match the original phase
     :param str hRanId: histogram Random Id
     :returns: a list of generated constraints
     '''
     from . import GSASIIobj as G2obj
     constrList = []
-    detTrans = np.abs(nl.det(Trans))  # volume ratio
     # Phase fractions (note volume term)
     IndpCon = [1.0,G2obj.G2VarObj([oRanId,hRanId,'Scale',None])]
-    DepCons = [detTrans,G2obj.G2VarObj([nRanId,hRanId,'Scale',None])]
+    DepCons = [Vratio,G2obj.G2VarObj([nRanId,hRanId,'Scale',None])]
     constrList.append([IndpCon,DepCons,None,None,'e'])
     # Isotropic size and mustrain (anisotropic constraints are messy)
     for name in ['Size;i','Mustrain;i']:

GSASII/GSASIImapvars.py

@@ -254,10 +254,24 @@ def warn(msg,cdict=None,val=None,prefix=None):
             if '_vary' in cdict:
                 warninfo['msg'] += prefix + 'new var expression: ' + _FormatConstraint(cdict,cdict.get('_name','New Var'))
             else:
-                warninfo['msg'] += prefix + 'constraint equation: ' + _FormatConstraint(cdict,val)
+                warninfo['msg'] += prefix + 'constraint equation:\n\t' + _FormatConstraint(cdict,val)
         if warninfo['msg']: warninfo['msg'] += '\n'
         warninfo['msg'] += '  ' + msg
-    
+    def showChanged(group):
+        errmsg = ''
+        changed = False
+        for rel in group:
+            orel = OrigSeq[rel]
+            if constrDict[rel] != OrigConstrDict[orel]: changed = True
+        if changed:
+            errmsg += '\n\nNote that these constraints were changed (due\n'
+            errmsg += 'to unrefined parameters?) Original constraints:\n'
+            for rel in group:
+                orel = OrigSeq[rel]
+                errmsg += '\n\t' + _FormatConstraint(OrigConstrDict[orel],OrigFixedList[orel])
+            errmsg += '\n\nSee the Warnings for details on unrefined parameters.\n'
+        return errmsg
+
     global dependentParmList,arrayList,invarrayList,indParmList,consNum
     # lists of parameters used for error reporting
     global undefinedVars # parameters that are used in equivalences but are not defined
@@ -284,6 +298,9 @@ def warn(msg,cdict=None,val=None,prefix=None):
     # Hold, Unvaried & Undefined parameters
     skipList = []
     invalidParms = []
+    OrigConstrDict = copy.deepcopy(constrDict)
+    OrigFixedList = copy.deepcopy(fixedList)
+    OrigSeq = list(range(len(constrDict)))
     for cnum,(cdict,fixVal) in enumerate(zip(constrDict,fixedList)):
         #constrVarList += [i for i in cdict if i not in constrVarList and not i.startswith('_')]
         valid = 0          # count of good parameters
@@ -334,11 +351,11 @@ def warn(msg,cdict=None,val=None,prefix=None):
             notDefList = []
 
         if noVaryList and cdict.get('_vary',True): # true for constr eq & varied New Var
-            msg = "parameter(s) not varied: "
+            noVarMsg = "parameter(s) not varied: "
             for i,v in enumerate(noVaryList):
-                if i != 0: msg += ', '
-                msg += v
-            warn(msg,cdict,fixVal,prefix='\nUnused ')
+                if i != 0: noVarMsg += ', '
+                noVarMsg += v
+            warn(noVarMsg,cdict,fixVal,prefix='\nUnused parameters in ')
         for l,m in ((zeroList,"have zero multipliers"), # show warning
                       (holdList,'set as "Hold"'),
                       #(noVaryList,"not varied"),
@@ -352,20 +369,19 @@ def warn(msg,cdict=None,val=None,prefix=None):
                 warn(msg,cdict,fixVal)
         if valid == 0: # no valid entries
             if seqHistNum is None:
-                warn('Ignoring constraint, contains no usable parameters',cdict,prefix='\nUnused ')
+                warn('Ignoring this constraint; contains no refined parameters',cdict,prefix='\nUnused ')
             skipList.append(cnum)
         elif problem: # mix of valid & refined and undefined items, cannot use this
             if cdict.get('_vary',True): # true for constr eq & varied New Var
                 warn('New Var constraint will be ignored',cdict,prefix='\nUnused ')
                 skipList.append(cnum)
             invalidParms += VarKeys(cdict)
         elif len(dropList) > 0: # mix of valid and problematic items, drop problem vars, but keep rest
-            if GSASIIpath.GetConfigValue('debug'): 
-                msg = ''
-                for v in dropList:
-                    if msg: msg += ' ,'
-                    msg += v
-                warn('removing: '+msg,cdict)
+            msg = ''
+            for v in dropList:
+                if msg: msg += ' ,'
+                msg += v
+            warn('removing unvaried: '+msg,cdict)
             value = fixedList[cnum]
             for var in dropList:   # do cleanup
                 # NB expressions in constraint multipliers have already been evaluated
@@ -375,19 +391,36 @@ def warn(msg,cdict=None,val=None,prefix=None):
                     value = float(value) - cdict[var]*parmDict[var]
                 del cdict[var]
             if float(value) != float(fixedList[cnum]): fixedList[cnum] = float(np.round(value,12))
-            if GSASIIpath.GetConfigValue('debug'):
-                warn('revised as: '+_FormatConstraint(constrDict[cnum],fixedList[cnum]))
+            warn('constraint revised as: '+_FormatConstraint(constrDict[cnum],fixedList[cnum]))
     for i in list(range(len(constrDict)-1,-1,-1)): # remove the dropped constraints
         if i in skipList:
             del constrDict[i]
             del fixedList[i]
+            del OrigSeq[i] # keep pointer to original lists
 
     for i in invalidParms: StoreHold(i,"Used in invalid constraint")
     if warning: warning += '\n'
     warning += warninfo['msg']
 
     groups,parmlist = GroupConstraints(constrDict)
 
+    # Look for duplicate constraints and eliminate one
+    for group,depPrmList in zip(groups,parmlist):
+        droplist = []
+        for i in group:
+            for j in group:
+                if i == j: continue
+                if i in droplist: continue
+                if j in droplist: continue
+                if constrDict[i] == constrDict[j] and fixedList[i] == fixedList[j]:
+                    droplist.append(j)
+        if droplist and warning: warning += '\n'
+        for j in sorted(droplist,reverse=True):
+            warning += '\nIgnoring duplicated constraint\n\t' 
+            warning += _FormatConstraint(constrDict[j],fixedList[j])
+            warning += '\n'
+            del group[group.index(j)]
+                    
     # now process each group and create the relations that are needed to form
     # a non-singular square matrix
     # Now check that all parameters are varied (probably do not need to do this
@@ -402,6 +435,7 @@ def warn(msg,cdict=None,val=None,prefix=None):
             errmsg += ") than parameters (" + str(len(depPrmList)) + ")\nin these constraints:"
             for rel in group:
                 errmsg += '\n\t'+ _FormatConstraint(constrDict[rel],fixedList[rel])
+            errmsg += showChanged(group)
             groupErrors += depPrmList
             continue # go on to next group
 
@@ -416,6 +450,7 @@ def warn(msg,cdict=None,val=None,prefix=None):
                 errmsg += "\nError expanding matrix with these constraints:"
                 for rel in group:
                     errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
+            errmsg += showChanged(group)
             groupErrors += depPrmList
             continue
 
@@ -426,6 +461,7 @@ def warn(msg,cdict=None,val=None,prefix=None):
             errmsg += "\nUnexpected singularity with constraints group (in Gram-Schmidt)"
             for rel in group:
                 errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
+            errmsg += showChanged(group)
             groupErrors += depPrmList
             continue
 
@@ -440,6 +476,7 @@ def warn(msg,cdict=None,val=None,prefix=None):
             errmsg += 'This is unexpected. Please report this (toby@anl.gov)'
             for rel in group:
                 errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
+            errmsg += showChanged(group)
             groupErrors += depPrmList
             continue
 
@@ -559,12 +596,12 @@ def CheckEquivalences(constrDict,varyList,fixedList,parmDict=None,seqHistNum=Non
     '''
 
     warninfo = {'msg':'', 'shown':-1}
-    def warnEqv(msg,cnum=None):
+    def warnEqv(msg,cnum=None,prefix='Problem with equivalence'):
         if cnum is not None and cnum != warninfo['shown']:
             warninfo['shown'] = cnum
-            if warninfo['msg']: warninfo['msg'] += '\n'
-            warninfo['msg'] += '\nProblem with equivalence: ' + _showEquiv(
-                dependentParmList[cnum],indParmList[cnum],invarrayList[cnum])
+            if warninfo['msg']: warninfo['msg'] += '\n\n'
+            warninfo['msg'] += f'{prefix}: '
+            warninfo['msg'] += _showEquiv(dependentParmList[cnum],indParmList[cnum],invarrayList[cnum])
         if warninfo['msg']: warninfo['msg'] += '\n'
         warninfo['msg'] += '  ' + msg
 
@@ -634,23 +671,27 @@ def warnEqv(msg,cnum=None):
                 if v in constrVarList+convVarList:
                     changed = True
                     msg = True
-                    warnEqv("Independent parameter "+str(v)+' used in constraint',cnum)
+                    warnEqv("Independent parameter "+str(v)+' used in constraint',cnum,
+                            prefix='Converting equivalence to constraint')
                     if cnum not in convertList: convertList.append(cnum)
             for v in varlist:
                 if v in multdepVarList:
                     changed = True
                     msg = True
-                    warnEqv("Dependent parameter "+str(v)+' repeated',cnum)
+                    warnEqv("Dependent parameter "+str(v)+' repeated',cnum,
+                            prefix='Converting equivalence to constraint')
                     if cnum not in convertList: convertList.append(cnum)
                 elif v in indepVarList:
                     changed = True
                     msg = True
-                    warnEqv("Dependent parameter "+str(v)+' used elsewhere as independent',cnum)
+                    warnEqv("Dependent parameter "+str(v)+' used elsewhere as independent',cnum,
+                            prefix='Converting equivalence to constraint')
                     if cnum not in convertList: convertList.append(cnum)
                 elif v in constrVarList+convVarList:
                     changed = True
                     msg = True
-                    warnEqv("Dependent parameter "+str(v)+' used in constraint',cnum)
+                    warnEqv("Dependent parameter "+str(v)+' used in constraint',cnum,
+                            prefix='Converting equivalence to constraint')
                     if cnum not in convertList: convertList.append(cnum)
             if msg:
                 warnEqv('Converting to "Constr"',cnum)
@@ -685,7 +726,8 @@ def warnEqv(msg,cnum=None):
             else:
                 msg = '  All parameters set as "Hold" '
             msg += "\n   Will ignore equivalence"
-            warnEqv(msg,cnum)
+            warnEqv(msg,cnum,prefix='Unused equivalence')
+
             removeList.append(cnum)
             continue
 
@@ -711,9 +753,10 @@ def warnEqv(msg,cnum=None):
                 removeList.append(cnum)
                 continue
             else:
-                msg = 'No parameters varied '
-            msg += "\n  Will ignore equivalence"
-            warnEqv(msg,cnum)
+                msg = 'Contains no varied parameters'
+            msg += "\n  Will ignore this equivalence"
+            warnEqv(msg,cnum,prefix='Unused equivalence')
+
             removeList.append(cnum)
             continue