Changeset 1948


Ignore:
Timestamp:
Aug 3, 2015 8:33:08 AM (7 years ago)
Author:
vondreele
Message:

modify SS special position code

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIspc.py

    r1917 r1948  
    14721472        return A
    14731473       
     1474    def DoFrac():
     1475        delt2 = np.eye(2)*0.001
     1476        FSC = np.ones(2,dtype='i')
     1477        VFSC = np.ones(2)
     1478        CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
     1479        if 'Crenel' in waveType:
     1480            dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze()
     1481        else:
     1482            dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
     1483        dFT = np.zeros_like(dF)
     1484        dFTP = []
     1485        for i in SdIndx:
     1486            sop = Sop[i]
     1487            ssop = SSop[i]
     1488            ssopinv = nl.inv(ssop[0])
     1489            mst = ssopinv[3][:3]
     1490            epsinv = ssopinv[3][3]
     1491            sdet = nl.det(sop[0])
     1492            ssdet = nl.det(ssop[0])
     1493            dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
     1494            dT = 1.0
     1495            if np.any(dtau%.5):
     1496                dT = np.tan(np.pi*np.sum(dtau%.5))
     1497            tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1498            fsc = np.ones(2,dtype='i')
     1499            if 'Crenel' in waveType:
     1500                dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze()
     1501                fsc = [1,1]
     1502            else:   #Fourier
     1503                dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
     1504                dFT = nl.det(sop[0])*dFT
     1505                dFT = dFT[:,np.argsort(tauT)]
     1506                dFT[0] *= ssdet
     1507                dFT[1] *= sdet
     1508                dFTP.append(dFT)
     1509           
     1510                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
     1511                    fsc = [1,1]
     1512                    CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
     1513                    FSC = np.zeros(2,dtype='i')
     1514                    return CSI,dF,dFTP
     1515                else:
     1516                    for i in range(2):
     1517                        if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
     1518                            fsc[i] = 1
     1519                        else:
     1520                            fsc[i] = 0
     1521                    FSC &= fsc
     1522                    print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc
     1523        n = -1
     1524        for i,F in enumerate(FSC):
     1525            if F:
     1526                n += 1
     1527                CSI[0][i] = n+1
     1528                CSI[1][i] = 1.0
     1529       
     1530        return CSI,dF,dFTP
     1531       
     1532    def DoXYZ():
     1533        delt6 = np.eye(6)*0.001
     1534        if 'Fourier' in waveType:
     1535            dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
     1536              #3x6x12 modulated position array (X,Spos,tau)& force positive
     1537            CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
     1538        elif waveType == 'Sawtooth':
     1539            CSI = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])]
     1540        elif waveType == 'ZigZag':
     1541            CSI = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])]
     1542        XSC = np.ones(6,dtype='i')
     1543        dXTP = []
     1544        for i in SdIndx:
     1545            sop = Sop[i]
     1546            ssop = SSop[i]
     1547            xsc = np.ones(6,dtype='i')
     1548            ssopinv = nl.inv(ssop[0])
     1549            mst = ssopinv[3][:3]
     1550            epsinv = ssopinv[3][3]
     1551            sdet = nl.det(sop[0])
     1552            ssdet = nl.det(ssop[0])
     1553            dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
     1554            dT = 1.0
     1555            if np.any(dtau%.5):
     1556                dT = np.tan(np.pi*np.sum(dtau%.5))
     1557            tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1558            if waveType == 'Fourier':
     1559                dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
     1560            elif waveType == 'Sawtooth':
     1561                dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
     1562            elif waveType == 'ZigZag':
     1563                dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]           
     1564            dXT = np.inner(sop[0],dXT.T)
     1565            dXT = np.swapaxes(dXT,1,2)
     1566            dXT[:,:3,:] *= ssdet
     1567            dXTP.append(dXT)
     1568            if waveType == 'Fourier':
     1569                if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
     1570                    xsc[3:6] = 0
     1571                    CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
     1572                        [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]                   
     1573                    if '(x)' in siteSym:
     1574                        CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
     1575                        if 'm' in siteSym and len(SdIndx) == 1:
     1576                            CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
     1577                    elif '(y)' in siteSym:
     1578                        CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
     1579                        if 'm' in siteSym and len(SdIndx) == 1:
     1580                            CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
     1581                    elif '(z)' in siteSym:
     1582                        CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
     1583                        if 'm' in siteSym and len(SdIndx) == 1:
     1584                            CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
     1585                    elif '(xy)' in siteSym:
     1586                        CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
     1587                        CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
     1588                elif np.allclose(dX[0,0,:],dXT[0,1,:]*sdet):
     1589                    if 'xy' in siteSym or '+-0' in siteSym:
     1590                        CSI[0][3:5] = [[12,0,0],[12,0,0]]
     1591                        CSI[1][3:5] = [[1.,0,0],[-sdet,0,0]]
     1592                        xsc[3:5] = 0
     1593                else:
     1594                    for i in range(3):
     1595                        if np.allclose(dX[i,i,:],dXT[i,i,:]*sdet):
     1596                            xsc[i] = 1
     1597                        else:
     1598                            xsc[i] = 0
     1599                        if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
     1600                            xsc[i+3] = 1
     1601                        else:
     1602                            xsc[i+3] = 0
     1603                XSC &= xsc
     1604                print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc
     1605        n = -1
     1606        print XSC
     1607        for i,X in enumerate(XSC):
     1608            if X:
     1609                n += 1
     1610                CSI[0][i][0] = n+1
     1611                CSI[1][i][0] = 1.0
     1612       
     1613        return CSI,dX,dXTP
     1614       
     1615    def DoUij():
     1616        delt12 = np.eye(12)*0.0001
     1617        dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
     1618        CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
     1619        USC = np.ones(12,dtype='i')
     1620        dUTP = []
     1621        for i in SdIndx:
     1622            sop = Sop[i]
     1623            ssop = SSop[i]
     1624            ssopinv = nl.inv(ssop[0])
     1625            mst = ssopinv[3][:3]
     1626            epsinv = ssopinv[3][3]
     1627            sdet = nl.det(sop[0])
     1628            ssdet = nl.det(ssop[0])
     1629            dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
     1630            dT = 1.0
     1631            if np.any(dtau%.5):
     1632                dT = np.tan(np.pi*np.sum(dtau%.5))
     1633            tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1634            usc = np.ones(12,dtype='i')
     1635            dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
     1636            dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
     1637            dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3)
     1638            dUT = np.array(U2Uij(dUijT))
     1639            dUT = dUT[:,:,np.argsort(tauT)]
     1640            dUTP.append(dUT)
     1641            if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
     1642                CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
     1643                [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
     1644                [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
     1645                [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
     1646                if 'mm2(x)' in siteSym:
     1647                    CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
     1648                    USC = [1,1,1,0,1,0,1,1,1,0,1,0]
     1649                elif '(xy)' in siteSym:
     1650                    CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
     1651                        [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
     1652                    CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
     1653                    USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
     1654                elif '(x)' in siteSym:
     1655                    CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
     1656                elif '(y)' in siteSym:
     1657                    CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
     1658                elif '(z)' in siteSym:
     1659                    CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
     1660                for i in range(6):
     1661                    if not USC[i]:
     1662                        CSI[0][i] = [0,0,0]
     1663                        CSI[1][i] = [0.,0.,0.]
     1664                        CSI[0][i+6] = [0,0,0]
     1665                        CSI[1][i+6] = [0.,0.,0.]
     1666            else:                       
     1667                for i in range(6):
     1668                    if np.allclose(dU[i,i,:],-dUT[i,i,:]):
     1669                        usc[i] = 1
     1670                    else:
     1671                        usc[i] = 0
     1672                    if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):
     1673                        usc[i+6] = 1
     1674                    else:
     1675                        usc[i+6] = 0
     1676                if np.any(dUT[0,1,:]):
     1677                    if '4/m' in siteSym:
     1678                        CSI[0][6:8] = [[12,0,0],[12,0,0]]
     1679                        if ssop[1][3]:
     1680                            CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
     1681                            usc[9] = 1
     1682                        else:
     1683                            CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
     1684                            usc[9] = 0
     1685                    elif '4' in siteSym:
     1686                        CSI[0][6:8] = [[12,0,0],[12,0,0]]
     1687                        CSI[0][:2] = [[11,0,0],[11,0,0]]
     1688                        if ssop[1][3]:
     1689                            CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
     1690                            CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
     1691                            usc[2] = 0
     1692                            usc[8] = 0
     1693                            usc[3] = 1
     1694                            usc[9] = 1
     1695                        else:
     1696                            CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
     1697                            CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
     1698                            usc[2] = 1
     1699                            usc[8] = 1
     1700                            usc[3] = 0               
     1701                            usc[9] = 0
     1702                    elif 'xy' in siteSym or '+-0' in siteSym:
     1703                        if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
     1704                            print np.allclose(dU[0,0,:],dUT[0,1,:]*sdet),sdet
     1705                            CSI[0][4:6] = [[12,0,0],[12,0,0]]
     1706                            CSI[0][6:8] = [[11,0,0],[11,0,0]]
     1707                            CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
     1708                            CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
     1709                            usc[4:6] = 0
     1710                            usc[6:8] = 0
     1711                       
     1712                print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc
     1713            USC &= usc
     1714        if not np.any(dtau%.5):
     1715            n = -1
     1716            for i,U in enumerate(USC):
     1717                if U:
     1718                    n += 1
     1719                    CSI[0][i][0] = n+1
     1720                    CSI[1][i][0] = 1.0
     1721
     1722        return CSI,dU,dUTP
     1723       
     1724    def DoUiso():
     1725        print 'Uiso'
     1726        delt4 = np.eye(4)*0.001
     1727       
    14741728    print 'super space group: ',SSGData['SSpGrp']
    14751729    CSI = {'Sfrac':[[[1,0],[2,0]],[[1.,0.],[1.,0.]]],
    … …  
    15071761            SGOps.append([-op[0],-op[1]%1.])
    15081762            SSGOps.append([-sop[0],-sop[1]%1.])
    1509     #build set of sym ops around special poasition       
     1763    #build set of sym ops around special position       
    15101764    SSop = []
    15111765    Sop = []
    … …  
    15251779    print 'special pos super operators: ',SSOpText
    15261780    #setup displacement arrays
    1527     tau = np.linspace(0,1,49,True)
    1528     delt2 = np.eye(2)*0.001
    1529     delt4 = np.eye(4)*0.001
    1530     delt6 = np.eye(6)*0.001
    1531     delt12 = np.eye(12)*0.0001
     1781    tau = np.linspace(-1,1,49,True)
    15321782    #make modulation arrays - one parameter at a time
    15331783    #site fractions
    1534     CSI['Sfrac'] = [np.zeros((2),dtype='i'),np.ones(2)]
    1535     if 'Crenel' in waveType:
    1536         dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze()
    1537     else:
    1538         dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
    1539     dFT = np.zeros_like(dF)
    1540     #positions       
    1541     if 'Fourier' in waveType:
    1542         dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
    1543           #3x6x12 modulated position array (X,Spos,tau)& force positive
    1544         CSI['Spos'] = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
    1545     elif waveType == 'Sawtooth':
    1546         CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])]
    1547     elif waveType == 'ZigZag':
    1548         CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])]
     1784    CSI['Sfrac'],dF,dFTP = DoFrac()
     1785    #positions
     1786    CSI['Spos'],dX,dXTP = DoXYZ()       
    15491787    #anisotropic thermal motion
    1550     dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
    1551     CSI['Sadp'] = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
    1552        
    1553     FSC = np.ones(2,dtype='i')
    1554     VFSC = np.ones(2)
    1555     XSC = np.ones(6,dtype='i')
    1556     USC = np.ones(12,dtype='i')
    1557     dFTP = []
    1558     dXTP = []
    1559     dUTP = []
    1560     for i in SdIndx:
    1561         sop = Sop[i]
    1562         ssop = SSop[i]
    1563         fsc = np.ones(2,dtype='i')
    1564         xsc = np.ones(6,dtype='i')
    1565         ssopinv = nl.inv(ssop[0])
    1566         mst = ssopinv[3][:3]
    1567         epsinv = ssopinv[3][3]
    1568         sdet = nl.det(sop[0])
    1569         ssdet = nl.det(ssop[0])
    1570         dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
    1571         dT = 1.0
    1572         if np.any(dtau%.5):
    1573             dT = np.tan(np.pi*np.sum(dtau%.5))
    1574         tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
    1575         if waveType == 'Fourier':
    1576             dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
    1577         elif waveType == 'Sawtooth':
    1578             dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
    1579         elif waveType == 'ZigZag':
    1580             dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]           
    1581         dXT = np.inner(sop[0],dXT.T)
    1582         dXT = np.swapaxes(dXT,1,2)
    1583         dXT[:,:3,:] *= ssdet
    1584         dXTP.append(dXT)
    1585         if waveType == 'Fourier':
    1586             if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
    1587                 CSI['Spos'] = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
    1588                     [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]                   
    1589                 if '(x)' in siteSym:
    1590                     CSI['Spos'][1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
    1591                     if 'm' in siteSym and len(SdIndx) == 1:
    1592                         CSI['Spos'][1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
    1593                 elif '(y)' in siteSym:
    1594                     CSI['Spos'][1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
    1595                     if 'm' in siteSym and len(SdIndx) == 1:
    1596                         CSI['Spos'][1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
    1597                 elif '(z)' in siteSym:
    1598                     CSI['Spos'][1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
    1599                     if 'm' in siteSym and len(SdIndx) == 1:
    1600                         CSI['Spos'][1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
    1601                 elif '(xy)' in siteSym:
    1602                     CSI['Spos'][0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
    1603                     CSI['Spos'][1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
    1604                 for i in range(3):
    1605                     if not XSC[i]:
    1606                         CSI['Spos'][0][i] = [0,0,0]
    1607                         CSI['Spos'][1][i] = [0.,0.,0.]
    1608                         CSI['Spos'][0][i+3] = [0,0,0]
    1609                         CSI['Spos'][1][i+3] = [0.,0.,0.]
    1610             elif np.allclose(dX[0,0,:],dXT[0,1,:]*sdet):
    1611                 if 'xy' in siteSym or '+-0' in siteSym:
    1612                     CSI['Spos'][0][3:5] = [[12,0,0],[12,0,0]]
    1613                     CSI['Spos'][1][3:5] = [[1.,0,0],[-sdet,0,0]]
    1614                     xsc[3:5] = 0
    1615             else:
    1616                 for i in range(3):
    1617                     if np.allclose(dX[i,i,:],dXT[i,i,:]*sdet):
    1618                         xsc[i] = 1
    1619                     else:
    1620                         xsc[i] = 0
    1621                     if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
    1622                         xsc[i+3] = 1
    1623                     else:
    1624                         xsc[i+3] = 0
    1625             XSC &= xsc
    1626            
    1627         fsc = np.ones(2,dtype='i')
    1628         if 'Crenel' in waveType:
    1629             dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze()
    1630             fsc = [1,1]
    1631         else:
    1632             dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
    1633             dFT = nl.det(sop[0])*dFT
    1634             dFT = dFT[:,np.argsort(tauT)]
    1635             dFT[0] *= ssdet
    1636             dFT[1] *= sdet
    1637             dFTP.append(dFT)
    1638        
    1639             if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
    1640                 fsc = [1,1]
    1641                 CSI['Sfrac'] = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
    1642                 for i in range(2):
    1643                     if not FSC[i]:
    1644                         CSI['Sfrac'][0][i] = [0,0]
    1645                         CSI['Sfrac'][1][i] = [0.,0.]
    1646             else:
    1647                 for i in range(2):
    1648                     if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
    1649                         fsc[i] = 1
    1650                     else:
    1651                         fsc[i] = 0
    1652         FSC &= fsc
    1653            
    1654         usc = np.ones(12,dtype='i')
    1655         dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
    1656         dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x49x3x3
    1657         dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3)
    1658         dUT = np.array(U2Uij(dUijT))
    1659         dUT = dUT[:,:,np.argsort(tauT)]
    1660         dUTP.append(dUT)
    1661         if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
    1662             CSI['Sadp'] = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
    1663             [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
    1664             [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
    1665             [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
    1666             if 'mm2(x)' in siteSym:
    1667                 CSI['Sadp'][1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
    1668                 USC = [1,1,1,0,1,0,1,1,1,0,1,0]
    1669             elif '(xy)' in siteSym:
    1670                 CSI['Sadp'][0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
    1671                     [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
    1672                 CSI['Sadp'][1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
    1673                 USC = [1,1,1,1,1,1,1,1,1,1,1,1]                             
    1674             elif '(x)' in siteSym:
    1675                 CSI['Sadp'][1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
    1676             elif '(y)' in siteSym:
    1677                 CSI['Sadp'][1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
    1678             elif '(z)' in siteSym:
    1679                 CSI['Sadp'][1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
    1680             for i in range(6):
    1681                 if not USC[i]:
    1682                     CSI['Sadp'][0][i] = [0,0,0]
    1683                     CSI['Sadp'][1][i] = [0.,0.,0.]
    1684                     CSI['Sadp'][0][i+6] = [0,0,0]
    1685                     CSI['Sadp'][1][i+6] = [0.,0.,0.]
    1686         else:                       
    1687             for i in range(6):
    1688                 if np.allclose(dU[i,i,:],-dUT[i,i,:]):
    1689                     usc[i] = 1
    1690                 else:
    1691                     usc[i] = 0
    1692                 if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):
    1693                     usc[i+6] = 1
    1694                 else:
    1695                     usc[i+6] = 0
    1696             if np.any(dUT[0,1,:]):
    1697                 if '4/m' in siteSym:
    1698                     CSI['Sadp'][0][6:8] = [[12,0,0],[12,0,0]]
    1699                     if ssop[1][3]:
    1700                         CSI['Sadp'][1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
    1701                         usc[9] = 1
    1702                     else:
    1703                         CSI['Sadp'][1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
    1704                         usc[9] = 0
    1705                 elif '4' in siteSym:
    1706                     CSI['Sadp'][0][6:8] = [[12,0,0],[12,0,0]]
    1707                     CSI['Sadp'][0][:2] = [[11,0,0],[11,0,0]]
    1708                     if ssop[1][3]:
    1709                         CSI['Sadp'][1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
    1710                         CSI['Sadp'][1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
    1711                         usc[2] = 0
    1712                         usc[8] = 0
    1713                         usc[3] = 1
    1714                         usc[9] = 1
    1715                     else:
    1716                         CSI['Sadp'][1][:2] = [[1.,0.,0.],[1.,0.,0.]]
    1717                         CSI['Sadp'][1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
    1718                         usc[2] = 1
    1719                         usc[8] = 1
    1720                         usc[3] = 0               
    1721                         usc[9] = 0
    1722                 elif 'xy' in siteSym or '+-0' in siteSym:
    1723                     if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
    1724                         print np.allclose(dU[0,0,:],dUT[0,1,:]*sdet),sdet
    1725                         CSI['Sadp'][0][4:6] = [[12,0,0],[12,0,0]]
    1726                         CSI['Sadp'][0][6:8] = [[11,0,0],[11,0,0]]
    1727                         CSI['Sadp'][1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
    1728                         CSI['Sadp'][1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
    1729                         usc[4:6] = 0
    1730                         usc[6:8] = 0
    1731                    
    1732             print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc
    1733         USC &= usc
    1734     if not np.any(dtau%.5):
    1735         n = -1
    1736         for i,U in enumerate(USC):
    1737             if U:
    1738                 n += 1
    1739                 CSI['Sadp'][0][i][0] = n+1
    1740                 CSI['Sadp'][1][i][0] = 1.0
    1741         if waveType == 'Fourier':
    1742             n = -1
    1743             for i,X in enumerate(XSC):
    1744                 if X:
    1745                     n += 1
    1746                     CSI['Spos'][0][i][0] = n+1
    1747                     CSI['Spos'][1][i][0] = 1.0
    1748         n = -1
    1749         for i,F in enumerate(FSC):
    1750             if F:
    1751                 n += 1
    1752                 CSI['Sfrac'][0][i] = n+1
    1753                 CSI['Sfrac'][1][i] = 1.0
    1754             else:
    1755                 CSI['Sfrac'][0][i] = 0
    1756                 CSI['Sfrac'][1][i] = 0.
     1788    CSI['Sadp'],dU,dUTP = DoUij()
    17571789    CSI['Spos'][0] = orderParms(CSI['Spos'][0])
    17581790    CSI['Sadp'][0] = orderParms(CSI['Sadp'][0])           
    17591791    if debug:
    1760         return CSI,[tau,tauT],[dF,dFTP],[dX,dXTP],[dU,dUTP]
     1792        return CSI,[tau,tau],[dF,dFTP],[dX,dXTP],[dU,dUTP]
    17611793    else:
    17621794        return CSI
Note: See TracChangeset for help on using the changeset viewer.