Ignore:
Timestamp:
Jul 21, 2021 3:04:47 PM (5 months ago)
Author:
toby
Message:

partial fix of magnetic constraints (still need to deal w/unvaried items; fullrmc minor stuff

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIconstrGUI.py

    r4983 r4998  
    17071707            Opr *= -1
    17081708        XOpr = np.inner(Opr,Trans)
     1709        invOpr = nl.inv(XOpr)
    17091710        for i,ix in enumerate(list(CSX[0])):
    17101711            if not ix:
     
    17131714            IndpCon = [1.0,G2obj.G2VarObj('%d::%s:%d'%(npId,name,ia))]
    17141715            DepCons = []
    1715             for iop,opval in enumerate(XOpr[i]):
    1716                 if opval:
     1716            for iop,opval in enumerate(invOpr[i]):
     1717                if abs(opval) > 1e-6:
    17171718                    DepCons.append([opval,G2obj.G2VarObj('%d::%s:%s'%(opId,xnames[iop],iat))])
    17181719            if len(DepCons) == 1:
    17191720                constraints['Phase'].append([DepCons[0],IndpCon,None,None,'e'])
    17201721            elif len(DepCons) > 1:
    1721                 for Dep in DepCons:
    1722                     Dep[0] *= -1
     1722                IndpCon[0] = -1.
    17231723                constraints['Phase'].append([IndpCon]+DepCons+[0.0,None,'c'])
    17241724        for name in ['Afrac','AUiso']:
     
    17611761
    17621762    # constraints on lattice parameters between phases
    1763 #    T = nl.inv(Trans).T
    1764     # T = Trans.T
    1765     # conMat = [
    1766     #     [T[0,0]**2,T[0,1]**2,T[0,2]**2,T[0,0]*T[0,1],T[0,0]*T[0,2],T[0,1]*T[0,2]],
    1767     #     [T[1,0]**2,T[1,1]**2,T[1,2]**2,T[1,0]*T[1,1],T[1,0]*T[1,2],T[1,1]*T[1,2]],
    1768     #     [T[2,0]**2,T[2,1]**2,T[2,2]**2,T[2,0]*T[2,1],T[2,0]*T[2,2],T[2,1]*T[2,2]],
    1769     #     [2.*T[0,0]*T[1,0],2.*T[0,1]*T[1,1],2.*T[0,2]*T[1,2],T[0,0]*T[1,1]+T[0,1]*T[1,0],T[0,0]*T[1,2]+T[0,2]*T[1,0],T[0,1]*T[1,2]+T[0,2]*T[1,1]],
    1770     #     [2.*T[0,0]*T[2,0],2.*T[0,1]*T[2,1],2.*T[0,2]*T[2,2],T[0,0]*T[2,1]+T[0,1]*T[2,0],T[0,0]*T[2,2]+T[0,2]*T[2,0],T[0,1]*T[2,2]+T[0,2]*T[2,1]],
    1771     #     [2.*T[1,0]*T[2,0],2.*T[1,1]*T[2,1],2.*T[1,2]*T[2,2],T[1,0]*T[2,1]+T[1,1]*T[2,0],T[1,0]*T[2,2]+T[1,2]*T[2,0],T[1,1]*T[2,2]+T[1,2]*T[2,1]]
    1772     #     ]
    1773     # Gnew = conMat * A:
    1774 #         T00**2*a0  T01**2*a1 T02**2*a2 T00*T01*a3    T00*T02*a4    T01*T02*a5
    1775 #         T10**2*a0  T11**2*a1 T12**2*a2 T10*T11*a3    T10*T12*a4    T11*T12*a5
    1776 #         T20**2*a0  T21**2*a1 T22**2*a2 T20*T21*a3    T20*T22*a4    T21*T22*a5
    1777 #         2*T00*T10*a0      2*T01*T11*a1     2*T02*T12*a2     (T00*T11 + T01*T10)*a3      (T00*T12 + T02*T10)*a4      (T01*T12 + T02*T11)*a5
    1778 #         2*T00*T20*a0      2*T01*T21*a1     2*T02*T22*a2     (T00*T21 + T01*T20)*a3      (T00*T22 + T02*T20)*a4      (T01*T22 + T02*T21)*a5
    1779 #         2*T10*T20*a0      2*T11*T21*a1     2*T12*T22*a2     (T10*T21 + T11*T20)*a3      (T10*T22 + T12*T20)*a4      (T11*T22 + T12*T21)*a5
    1780     # Generated as symbolic code using:
    1781     # import sym
    1782     # A0, A1, A2, A3, A4, A5 = sym.symbols('A0, A1, A2, A3, A4, A5')
    1783     # G = sym.Matrix([ [A0,    A3/2.,  A4/2.], [A3/2.,  A1,    A5/2.], [A4/2., A5/2.,    A2]])
    1784     # transformation matrix
    1785     # T00, T10, T20, T01, T11, T21, T02, T12, T22 = sym.symbols('T00, T10, T20, T01, T11, T21, T02, T12, T22')
    1786     # Tr = sym.Matrix([ [T00, T10, T20], [T01, T11, T21], [T02, T12, T22],])
    1787     # Gnew = (Tr.T*G)*Tr
    1788    
    1789     #print('old A',G2lat.cell2A(oldPhase['General']['Cell'][1:7]))
    1790     #print('new A',G2lat.cell2A(newPhase['General']['Cell'][1:7]))
    1791    
    1792 #this is still incorrect for hex/trig/ortho/tetragonal --> monoclinic
    1793    
    1794 #    for iAnew,Asi in enumerate(['A0','A1','A2','A3','A4','A5']): # loop through A[i] for new cell
    1795 #        Nparm = str(npId) + '::' + Asi
    1796 #        if Nparm != SetUniqAj(npId,iAnew,nSGData):
    1797 #            continue # skip: Ai constrained from Aj or must be zero
    1798 #        multDict = {}
    1799 #        for iAorg in range(6):
    1800 #            cA = conMat[iAnew][iAorg] # coeff for A[i] in constraint matrix
    1801 #            if abs(cA) < 1.e-8: continue
    1802 #            parm = SetUniqAj(opId,iAorg,oSGData) # translate to unique A[i] in original cell
    1803 #            if not parm: continue # must be zero
    1804 #            # sum coeff
    1805 #            if parm in multDict:
    1806 #                multDict[parm] += cA
    1807 #            else:
    1808 #                multDict[parm] = cA
    1809 #        # any non-zero multipliers?
    1810 #        maxMult = 0
    1811 #        for i in multDict:
    1812 #            maxMult = max(maxMult,abs(multDict[i]))
    1813 #        if maxMult <= 0:  # Nparm computes as zero; Fix this parameter
    1814 #            constraints['Phase'] += [[
    1815 #                [0.0,G2obj.G2VarObj(Nparm)],
    1816 #                None,None,'h']]
    1817 #        elif len(multDict) == 1:        # create equivalence
    1818 #            key = list(multDict.keys())[0]
    1819 #            constraints['Phase'] += [[
    1820 #                [1.0,G2obj.G2VarObj(key)],
    1821 #                [multDict[key],G2obj.G2VarObj(Nparm)],
    1822 #                None,None,'e']]
    1823 #        else:                           # create constraint
    1824 #            constr = [[-1.0,G2obj.G2VarObj(Nparm)]]
    1825 #            for key in multDict:
    1826 #                constr += [[multDict[key],G2obj.G2VarObj(key)]]
    1827 #            constr += [0.0,None,'c']
    1828 #            constraints['Phase'] += [constr]
    1829    
     1763    Aold = G2lat.cell2A(oldPhase['General']['Cell'][1:7])
     1764    if True: # debug
     1765        constraints['Phase'] += G2lat.GenCellConstraints(Trans,opId,npId,Aold,True)
     1766        print('old A*',G2lat.cell2A(oldPhase['General']['Cell'][1:7]))
     1767        print('new A*',G2lat.cell2A(newPhase['General']['Cell'][1:7]))
     1768        print('old cell',oldPhase['General']['Cell'][1:7])
     1769        print('new cell',newPhase['General']['Cell'][1:7])
     1770    else:
     1771        constraints['Phase'] += G2lat.GenCellConstraints(Trans,opId,npId,Aold)
     1772
    18301773    # constraints on HAP Scale, etc.
    18311774    for hId,hist in enumerate(UseList):    #HAP - seems OK
Note: See TracChangeset for help on using the changeset viewer.