Changeset 1635


Ignore:
Timestamp:
Feb 4, 2015 4:38:12 PM (9 years ago)
Author:
vondreele
Message:

modulation functions & constraints

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1630 r1635  
    709709def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata):
    710710    import scipy.special as sp
    711     for iwt,wt in enumerate(waveTypes):    #atom loop!
    712         if wt == 'Fourier':
    713             A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])])
    714             HdotA = twopi*np.inner(A.T,SSUniq.T[:3].T)
    715             m = SSUniq.T[3]
    716             Gp = sp.jn(-m,HdotA)
    717     return np.real(Gp),np.imag(Gp)
     711    import scipy.integrate as si
     712    m = SSUniq.T[3]
     713    nh = np.zeros(1)
     714    if XSSdata.ndim > 2:
     715        nh = np.arange(XSSdata.shape[1])       
     716    M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
     717    A = np.array(XSSdata[:3])
     718    B = np.array(XSSdata[3:])
     719    HdotA = (np.inner(A.T,SSUniq.T[:3].T)+SSPhi)
     720    HdotB = (np.inner(B.T,SSUniq.T[:3].T)+SSPhi)
     721    GpA = sp.jn(M[:,np.newaxis],twopi*HdotA)
     722    GpB = sp.jn(M[:,np.newaxis],twopi*HdotB)*(1.j)**M
     723    Gp = np.sum(GpA+GpB,axis=0)
     724    return np.real(Gp).T,np.imag(Gp).T
     725   
     726def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata):
     727    import scipy.special as sp
     728    m = SSUniq.T[3]
     729    nh = np.zeros(1)
     730    if XSSdata.ndim > 2:
     731        nh = np.arange(XSSdata.shape[1])       
     732    M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
     733    A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])])
     734    HdotA = twopi*(np.inner(SSUniq.T[:3].T,A.T)+SSPhi)
     735    Gpm = sp.jn(M[:,np.newaxis,:]-1,HdotA)
     736    Gpp = sp.jn(M[:,np.newaxis,:]+1,HdotA)
     737    if Gpm.ndim > 3: #sum over multiple harmonics
     738        Gpm = np.sum(Gpm,axis=0)
     739        Gpp = np.sum(Gpp,axis=0)
     740    dGpdk = 0.5*(Gpm+Gpp)
     741    return np.real(dGpdk),np.imag(dGpdk)
    718742   
    719743def posFourier(tau,psin,pcos):

  • trunk/GSASIIphsGUI.py

    r1633 r1635  
    20812081            waveSizer.Add(waveHead)
    20822082            if len(waveBlk):
     2083                nFour = 0
    20832084                for iwave,wave in enumerate(waveBlk):
     2085                    if waveType == 'Fourier':
     2086                        nFour += 1
    20842087                    if not iwave:
    2085                         CSI = G2spc.GetSSfxuinel(waveType,xyz,uij,SGData,SSGData)
     2088                        CSI = G2spc.GetSSfxuinel(waveType,nFour,xyz,SGData,SSGData)
    20862089                    else:
    2087                         CSI = G2spc.GetSSfxuinel('Fourier',xyz,uij,SGData,SSGData)
     2090                        CSI = G2spc.GetSSfxuinel('Fourier',nFour,xyz,SGData,SSGData)
    20882091                    waveName = 'Fourier'
    20892092                    if Stype == 'Sfrac':
    … …  
    21072110                    waveSizer.Add(wx.StaticText(waveData,label=' %s  parameters: %s'%(waveName,str(names).rstrip(']').lstrip('[').replace("'",''))),0,WACV)
    21082111                    for ival,val in enumerate(wave[0]):
    2109                         if CSI[Stype][0][ival] == 0:
    2110                             waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_READONLY)
    2111                             waveVal.SetBackgroundColour(VERY_LIGHT_GREY)
    2112                         else:
     2112                        if any(CSI[Stype][0][ival]):
    21132113                            waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_PROCESS_ENTER)
    21142114                            waveVal.Bind(wx.EVT_TEXT_ENTER,OnWaveVal)
    21152115                            waveVal.Bind(wx.EVT_KILL_FOCUS,OnWaveVal)
    21162116                            Indx[waveVal.GetId()] = [iatm,Stype,iwave,ival]
     2117                        else:
     2118                            waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_READONLY)
     2119                            waveVal.SetBackgroundColour(VERY_LIGHT_GREY)
    21172120                        Waves.Add(waveVal,0,WACV)
    21182121                        if len(wave[0]) > 6 and ival == 5:

  • trunk/GSASIIplot.py

    r1633 r1635  
    222222        Toolbar.__init__(self,plotCanvas)
    223223        self.plotCanvas = plotCanvas
    224         POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6 # remove one button
     224        POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6 # remove one button, nos. start at 1!
    225225        self.DeleteToolByPos(POSITION_OF_CONFIGURE_SUBPLOTS_BTN)
    226226        self.parent = self.GetParent()

  • trunk/GSASIIspc.py

    r1630 r1635  
    2121import math
    2222import sys
     23import copy
    2324import os.path as ospath
    2425
    … …  
    2930npsind = lambda x: np.sin(x*np.pi/180.)
    3031npcosd = lambda x: np.cos(x*np.pi/180.)
    31 DEBUG = True
     32DEBUG = False
    3233   
    3334################################################################################
    … …  
    14221423    return CSuinel[indx[1]]
    14231424   
    1424 def GetSSfxuinel(XYZ,UIJ,SGData,SSGData):
    1425     CSI = {'Sfrac':[[1,2],[1.,1.]],'Spos':[[1,2,3, 4,5,6],[1.,1.,1., 1.,1.,1.]],    #sin & cos
    1426         'Sadp':[[1,2,3,4,5,6, 7,8,9,10,11,12],[1.,1.,1.,1.,1.,1., 1.,1.,1.,1.,1.,1.]],
    1427         'Smag':[[1,2,3, 4,5,6],[1.,1.,1., 1.,1.,1.]]}
    1428     deltx = np.ones((3,4))*.01
    1429     deltx[:3,:3] = np.eye((3))*.001
    1430     deltu = np.eye((6))*.0001
     1425def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False):
     1426   
     1427    def fracCrenel(tau,Toff,Twid):
     1428        Tau = (tau-Toff)%1.
     1429        A = np.where(Tau<Twid,1.,0.)
     1430        return A
     1431       
     1432    def fracFourier(tau,nH,fsin,fcos):
     1433        SA = np.sin(2.*nH*np.pi*tau)
     1434        CB = np.cos(2.*nH*np.pi*tau)
     1435        A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis]
     1436        B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis]
     1437        return A+B
     1438       
     1439    def posFourier(tau,nH,psin,pcos):
     1440        SA = np.sin(2*nH*np.pi*tau)
     1441        CB = np.cos(2*nH*np.pi*tau)
     1442        A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis]
     1443        B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis]
     1444        return A+B   
     1445
     1446    def posSawtooth(tau,Toff,slopes):
     1447        Tau = (tau-Toff[:,np.newaxis])%1.
     1448        A = slopes[:,:,np.newaxis]*Tau
     1449        return A
     1450   
     1451    def posZigZag(tau,Toff,slopes):
     1452        Tau = (tau-Toff[:,np.newaxis])%1.
     1453        A = np.where(Tau <= 0.5,slopes[:,:,np.newaxis]*Tau,slopes[:,:,np.newaxis]*(1.-Tau))
     1454        return A
     1455       
     1456    print 'super space group: ',SSGData['SSpGrp']
     1457    CSI = {'Sfrac':[[[1,0],[2,0]],[[1.,0.],[1.,0.]]],
     1458        'Spos':[[[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],
     1459            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],    #sin & cos
     1460        'Sadp':[[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
     1461            [7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],
     1462            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
     1463            [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],
     1464        'Smag':[[[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],
     1465            [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]],}
    14311466    xyz = np.array(XYZ)%1.
    14321467    xyzt = np.array(XYZ+[0,])%1.
    1433     uij = np.array(UIJ)
    1434     SGOps = SGData['SGOps']
    1435     SSGOps = SSGData['SSGOps']
    1436     ssop = SSGOps[0]
    1437     sop = SGOps[0]
    1438     for iop,Op in enumerate(SGOps):
     1468    SGOps = copy.deepcopy(SGData['SGOps'])
     1469    siteSym = SytSym(XYZ,SGData)[0].strip()
     1470    print 'siteSym: ',siteSym
     1471    if siteSym == '1':   #"1" site symmetry
     1472        if debug:
     1473            return CSI,None,None,None,None
     1474        else:
     1475            return CSI
     1476    elif siteSym == '-1':   #"-1" site symmetry
     1477        CSI['Sfrac'][0] = [[1,0],[0,0]]
     1478        CSI['Spos'][0] = [[1,0,0],[2,0,0],[3,0,0], [0,0,0],[0,0,0],[0,0,0]]
     1479        CSI['Sadp'][0] = [[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],
     1480        [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]]
     1481        if debug:
     1482            return CSI,None,None,None,None
     1483        else:
     1484            return CSI
     1485    SSGOps = copy.deepcopy(SSGData['SSGOps'])
     1486    #expand ops to include inversions if any
     1487    if SGData['SGInv']:
     1488        for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
     1489            SGOps.append([-op[0],-op[1]%1.])
     1490            SSGOps.append([-sop[0],-sop[1]%1.])
     1491    #build set of sym ops around special poasition       
     1492    SSop = []
     1493    Sop = []
     1494    for iop,Op in enumerate(SGOps):         
    14391495        nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
    1440         if SGData['SGInv'] and np.allclose(xyz,-nxyz%1.,1.e-6):
    1441             ssop = SSGOps[iop]
    1442             ssop = [-ssop[0],-ssop[1]%1.]
    1443             sop = [-Op[0],-Op[1]%1.]
    1444             break
    1445         elif np.allclose(xyz,nxyz,1.e-6) and iop:
    1446             ssop = SSGOps[iop]
    1447             sop = SGOps[iop]
    1448             break
    1449     siteSym = SytSym(XYZ,SGData)[0].strip().split('(')[0]
    1450     OpText =  MT2text(sop).replace(' ','')
    1451     SSOptext = SSMT2text(ssop).replace(' ','')
    1452     if siteSym == '1':   #"1" site symmetry
     1496        if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
     1497            SSop.append(SSGOps[iop])
     1498            Sop.append(SGOps[iop])
     1499    OpText =  [MT2text(s).replace(' ','') for s in Sop]         #debug?
     1500    SSOpText = [SSMT2text(ss).replace(' ','') for ss in SSop]   #debug?
     1501    print 'special pos super operators: ',SSOpText
     1502    #setup displacement arrays
     1503    tau = np.linspace(0,1,49,True)
     1504    delt2 = np.eye(2)*0.001
     1505    delt4 = np.eye(4)*0.001
     1506    delt6 = np.eye(6)*0.001
     1507    delt12 = np.eye(12)*0.0001
     1508    #make modulation arrays - one parameter at a time
     1509    #site fractions
     1510    CSI['Sfrac'] = [np.zeros((2),dtype='i'),np.ones(2)]
     1511    if 'Crenel' in waveType:
     1512        dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze()
     1513    else:
     1514        dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
     1515    dFT = np.zeros_like(dF)
     1516    #positions       
     1517    if 'Fourier' in waveType:
     1518        dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
     1519          #3x6x12 modulated position array (X,Spos,tau)& force positive
     1520        CSI['Spos'] = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
     1521    elif waveType == 'Sawtooth':
     1522        CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])]
     1523    elif waveType == 'ZigZag':
     1524        CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])]
     1525    #anisotropic thermal motion
     1526    dU = posFourier(tau,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
     1527    CSI['Sadp'] = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
     1528       
     1529    FSC = np.ones(2,dtype='i')
     1530    VFSC = np.ones(2)
     1531    XSC = np.ones(6,dtype='i')
     1532    USC = np.ones(12,dtype='i')
     1533    dFTP = []
     1534    dXTP = []
     1535    dUTP = []
     1536    for sop,ssop in zip(Sop,SSop):
     1537        fsc = np.ones(2,dtype='i')
     1538        xsc = np.ones(6,dtype='i')
     1539        ssopinv = nl.inv(ssop[0])
     1540        mst = ssopinv[3][:3]
     1541        epsinv = ssopinv[3][3]
     1542        tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1543        if waveType == 'Fourier':
     1544            dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
     1545        elif waveType == 'Sawtooth':
     1546            dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
     1547        elif waveType == 'ZigZag':
     1548            dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]           
     1549        dXT = np.inner(sop[0],dXT.T)
     1550        dXT = np.swapaxes(dXT,1,2)
     1551#        dXT = dXT[:,:,np.argsort(tauT)]
     1552        dXTP.append(dXT)
     1553        if waveType == 'Fourier':
     1554            for i in range(3):
     1555                if np.allclose(dX[i,i,:],dXT[i,i,:]):
     1556                    xsc[i] = 1
     1557                else:
     1558                    xsc[i] = 0
     1559                if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
     1560                    xsc[i+3] = 1
     1561                else:
     1562                    xsc[i+3] = 0
     1563            XSC &= xsc
     1564           
     1565        fsc = np.ones(2,dtype='i')
     1566        vfsc = np.ones(2)
     1567        if 'Crenel' in waveType:
     1568            dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze()
     1569        else:
     1570            dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
     1571        dFT = nl.det(sop[0])*dFT
     1572        dFT = dFT[:,np.argsort(tauT)]
     1573        dFTP.append(dFT)
     1574        for i in range(2):
     1575            if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
     1576                fsc[i] = 1
     1577                vfsc[i] = 1.0               
     1578            else:
     1579                fsc[i] = 0
     1580                vfsc[i] = 0.
     1581        FSC &= fsc
     1582        VFSC = vfsc
     1583           
     1584        usc = np.ones(12,dtype='i')
     1585        # make 12x12x4x4 with tau layers?
     1586        dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x12 array
     1587        dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3)    #convert dUT to 12x12x3x3
     1588        dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3)
     1589        dUT = np.array(U2Uij(dUijT))
     1590        dUT = dUT[:,:,np.argsort(tauT)]
     1591        dUTP.append(dUT)
     1592        for i in range(6):
     1593            if np.allclose(dU[i,i,:],dUT[i,i,:]):
     1594                usc[i] = 1
     1595            else:
     1596                usc[i] = 0
     1597            if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]):
     1598                usc[i+6] = 1
     1599            else:
     1600                usc[i+6] = 0
     1601        USC &= usc
     1602    n = -1
     1603    for i,U in enumerate(USC):
     1604        if U:
     1605            n += 1
     1606            CSI['Sadp'][0][i][0] = n+1
     1607            CSI['Sadp'][1][i][0] = 1.0
     1608    if waveType == 'Fourier':
     1609        n = -1
     1610        for i,X in enumerate(XSC):
     1611            if X:
     1612                n += 1
     1613                CSI['Spos'][0][i][0] = n+1
     1614                CSI['Spos'][1][i][0] = 1.0
     1615    n = -1
     1616    for i,[F,VF] in enumerate(zip(FSC,VFSC)):
     1617        if F:
     1618            n += 1
     1619            CSI['Sfrac'][0][i] = n+1
     1620            CSI['Sfrac'][1][i] = VF
     1621        else:
     1622            CSI['Sfrac'][0][i] = 0
     1623            CSI['Sfrac'][1][i] = 0.           
     1624    if debug:
     1625        return CSI,[tau,tauT],[dF,dFTP],[dX,dXTP],[dU,dUTP]
     1626    else:
    14531627        return CSI
    1454     elif siteSym == '-1':   #"-1" site symmetry
    1455         CSI['Spos'][0] = [1,2,3, 0,0,0]
    1456         CSI['Sadp'][0] = [0,0,0,0,0,0, 1,2,3,4,5,6]
    1457         return CSI
    1458 #    print siteSym,OpText,SSOptext   
    1459     UniqAx = {'a':'a','b':'b','c':'g'}
    1460     if SGData['SGLaue'] == '2/m':
    1461         if UniqAx[SGData['SGUniq']] in SSGData['modSymb']:   #e.g. (0b0)
    1462             if 's' in SSGData['SSpGrp'].split('(')[1]:
    1463                 if siteSym == 'm':
    1464                     CSI['Spos'][0] = [1,0,2, 0,3,0]
    1465                     CSI['Sadp'][0] = [1,2,3,0,4,0, 0,0,0,5,0,6]
    1466                 elif siteSym == '2/m':
    1467                     CSI['Spos'][0] = [1,0,2, 0,0,0]
    1468                     CSI['Sadp'][0] = [0,0,0,0,0,0, 0,0,0,1,0,2]
    1469                 elif siteSym == '2':
    1470                     CSI['Spos'][0] = [1,0,2, 3,0,4]               
    1471                     CSI['Sadp'][0] = [0,0,0,1,0,2, 0,0,0,3,0,4]
    1472             elif '1/2' in SSGData['modSymb']:    #e.g. (0b1/2)
    1473                 if siteSym == 'm':
    1474                     CSI['Spos'][0] = [0,0,0, 0,0,0]
    1475                     CSI['Sadp'][0] = [0,0,0,0,0,0, 0,0,0,0,0,0]
    1476                 elif siteSym == '2/m':
    1477                     CSI['Spos'][0] = [1,0,2, 0,0,0]
    1478                     CSI['Sadp'][0] = [0,0,0,0,0,0, 1,2,3,0,4,0]
    1479                 elif siteSym == '2':
    1480                     CSI['Spos'][0] = [1,0,2, 3,0,4]
    1481                     CSI['Sadp'][0] = [0,0,0,1,0,2, 0,0,0,3,0,4]
    1482             else:
    1483                 if siteSym == 'm':
    1484                     CSI['Spos'][0] = [0,1,0, 2,0,3]
    1485                     CSI['Sadp'][0] = [0,0,0,1,0,2, 3,4,5,0,6,0]
    1486                 elif siteSym == '2/m':
    1487                     CSI['Spos'][0] = [0,1,0, 0,0,0]
    1488                     CSI['Sadp'][0] = [0,0,0,0,0,0, 1,2,3,0,4,0]
    1489                 elif siteSym == '2':
    1490                     CSI['Spos'][0] = [0,1,0, 0,2,0]
    1491                     CSI['Sadp'][0] = [1,2,3,0,4,0, 5,6,7,0,8,0]
    1492            
    1493         else:   #e.g. (a0g)
    1494             if 's' in SSGData['SSpGrp'].split('(')[1]:
    1495                 if siteSym == 'm':
    1496                     CSI['Spos'][0] = [0,1,0, 0,2,0]
    1497                     CSI['Sadp'][0] = [0,0,0,1,0,2, 0,0,0,3,0,4]
    1498                 elif siteSym == '2/m':
    1499                     CSI['Spos'][0] = [0,1,0, 0,0,0]
    1500                     CSI['Sadp'][0] = [0,0,0,0,0,0, 0,0,0,1,0,2]
    1501                 elif siteSym == '2':
    1502                     CSI['Spos'][0] = [0,1,0, 2,0,3]
    1503                     CSI['Sadp'][0] = [1,2,3,0,4,0, 0,0,0,5,0,6]
    1504             else:
    1505                 if siteSym == 'm':
    1506                     CSI['Spos'][0] = [1,0,2, 3,0,4]
    1507                     CSI['Sadp'][0] = [1,2,3,0,4,0, 5,6,7,0,8,0]
    1508                 elif siteSym == '2/m':
    1509                     CSI['Spos'][0] = [1,0,2, 0,0,0]
    1510                     CSI['Sadp'][0] = [0,0,0,0,0,0, 1,2,3,0,4,0]
    1511                 elif siteSym == '2':
    1512                     CSI['Spos'][0] = [1,0,2, 0,3,0]
    1513                     CSI['Sadp'][0] = [0,0,0,1,0,2, 3,4,5,0,6,0]
    1514         Sx = CSI['Spos'][0]
    1515         Su = CSI['Sadp'][0]
    1516         if SGData['SGUniq'] == 'a':
    1517             CSI['Spos'][0] = [Sx[1],Sx[2],Sx[0], Sx[4],Sx[5],Sx[3]]
    1518             CSI['Sadp'][0] = [Su[1],Sx[2],Sx[0],Su[4],Su[5],Su[3], Su[7],Su[8],Su[6],Su[10],Su[11],Su[9]]
    1519         elif SGData['SGUniq'] == 'c':
    1520             CSI['Spos'][0] = [Sx[2],Sx[0],Sx[1], Sx[5],Sx[3],Sx[4]]
    1521             CSI['Sadp'][0] = [Su[2],Su[0],Su[1],Su[5],Su[3],Su[4], Su[8],Su[6],Su[7],Su[11],Su[9],Su[10]]
    1522 #        return CSI
    1523 #    elif SGData['SGLaue'] == 'mmm':
    1524 #    elif SGData['SGLaue'] in ['4/m','4/mmm']:
    1525 #    elif SGData['SGLaue'] in ['3','3m1','31m']:
    1526 #    elif SGData['SGLaue'] in ['6/m','6/mmm']:
    1527 #       
    1528     xsin = np.zeros(3,dtype='i')
    1529     xcos = np.zeros(3,dtype='i')
    1530     usin = np.zeros(6,dtype='i')
    1531     ucos = np.zeros(6,dtype='i')
    1532     csi = np.ones((6),dtype='i')*-1
    1533     for i,idelt in enumerate(deltx):
    1534 #        print 'idelt:',idelt
    1535         nxyzt = np.inner(ssop[0],(xyzt+idelt))+ssop[1]
    1536         nxyzt[3] -= ssop[1][3]
    1537 #        print 'nxyz',nxyzt
    1538         xsin[i] = np.allclose((xyzt-idelt),nxyzt,1.e-6)
    1539 #        print 'sin ',(xyzt-idelt),xsin[i]
    1540         xcos[i] = np.allclose((xyzt+idelt),nxyzt,1.e-6)
    1541 #        print 'cos ',(xyzt+idelt),xcos[i]
    1542     n = -1
    1543     for i,isin in enumerate(xsin):
    1544         if isin:
    1545             n += 1
    1546             csi[i] = n
    1547     for i,icos in enumerate(xcos):
    1548         if icos:
    1549             n += 1
    1550             csi[i+3] = n
    1551 #    print csi
    1552 #    print CSI['Spos'][0]
    1553 #    print xsin,xcos
    1554     for i,idelt in enumerate(deltu):
    1555         nuij = U2Uij(np.inner(sop[0],np.inner(np.abs(Uij2U(uij+idelt)),sop[0])))
    1556         usin[i] = np.equal(np.abs(uij-idelt),nuij)[i]
    1557         ucos[i] = np.equal(np.abs(uij+idelt),nuij)[i]
    1558 #    print CSI['Sadp'][0]
    1559 #    print usin,ucos
    1560     return CSI
    15611628   
    15621629def MustrainNames(SGData):

  • trunk/GSASIIstrIO.py

    r1630 r1635  
    10921092                    waveType = AtomSS['waveType']
    10931093                    phaseDict[pfx+'waveType:'+str(i)] = waveType   
    1094                     CSI = G2spc.GetSSfxuinel(at[cx:cx+3],at[cia+2:cia+8],SGData,SSGData)
    10951094                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
    10961095                        Waves = AtomSS[Stype]
    10971096                        uId,uCoef = CSI[Stype]
    10981097                        for iw,wave in enumerate(Waves):
     1098                            if not iw:
     1099                                CSI = G2spc.GetSSfxuinel(waveType,at[cx:cx+3],SGData,SSGData)
     1100                            else:
     1101                                CSI = G2spc.GetSSfxuinel('Fourier',at[cx:cx+3],SGData,SSGData)
    10991102                            stiw = str(i)+':'+str(iw)
    11001103                            if Stype == 'Spos':
    … …  
    14961499        cx,ct,cs,cia = General['AtomPtrs']
    14971500        print >>pFile,'\n Modulation waves'
    1498 #        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tzero','Xslope','Yslope','Zslope'],
    1499 #            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
    1500 #            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
    1501 #        print >>pFile,135*'-'
    1502 #        for i,at in enumerate(Atoms):
    1503 #            AtomSS = at[-1]['SS1']
    1504 #            for Stype in ['Sfrac','Spos','Sadp','Smag']:
    1505 #                Waves = AtomSS[Stype]
    1506 #                if len(Waves):
    1507 #                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
    1508 #                        %(at[ct-1],at[cs],Stype,AtomSS['waveType'])
    1509 #                    line = ''
    1510 #                    for item in names[Stype]:
    1511 #                        line += '%8s '%(item)
    1512 #                    print >>pFile,line
    1513 #                for wave in Waves:                   
    1514 #                    line = ''
    1515 #                    for item in wave[0]:
    1516 #                        line += '%8.4f '%(item)
    1517 #                    line += ' Refine? '+str(wave[1])
    1518 #                    print >>pFile,line
     1501        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tzero','Xslope','Yslope','Zslope'],
     1502            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
     1503            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
     1504        print >>pFile,135*'-'
     1505        for i,at in enumerate(Atoms):
     1506            AtomSS = at[-1]['SS1']
     1507            waveType = AtomSS['waveType']
     1508            for Stype in ['Sfrac','Spos','Sadp','Smag']:
     1509                Waves = AtomSS[Stype]
     1510                if len(Waves):
     1511                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
     1512                        %(at[ct-1],at[cs],Stype,waveType)
     1513                    line = ''
     1514                    for iw,wave in enumerate(Waves):
     1515                        stiw = ':'+str(i)+':'+str(iw)
     1516                        namstr = '  names :'
     1517                        valstr = '  values:'
     1518                        sigstr = '  esds  :'
     1519                        if Stype == 'Spos':
     1520                            nt = 6
     1521                            ot = 0
     1522                            if waveType in ['Sawtooth','ZigZag'] and not iw:
     1523                                nt = 4
     1524                                ot = 6
     1525                            for j in range(nt):
     1526                                name = names['Spos'][j+ot]
     1527                                namstr += '%12s'%(name)
     1528                                valstr += '%12.4f'%(wave[0][j])
     1529                                if name+stiw in wavesSig:
     1530                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
     1531                                else:
     1532                                    sigstr += 12*' '
     1533                        elif Stype == 'Sfrac':
     1534                            ot = 0
     1535                            if 'Crenel' in waveType and not iw:
     1536                                ot = 2
     1537                            for j in range(2):
     1538                                name = names['Sfrac'][j+ot]
     1539                                namstr += '%12s'%(names['Sfrac'][j+ot])
     1540                                valstr += '%12.4f'%(wave[0][j])
     1541                                if name+stiw in wavesSig:
     1542                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
     1543                                else:
     1544                                    sigstr += 12*' '
     1545                        elif Stype == 'Sadp':
     1546                            for j in range(12):
     1547                                name = names['Sadp'][j]
     1548                                namstr += '%10s'%(names['Sadp'][j])
     1549                                valstr += '%10.6f'%(wave[0][j])
     1550                                if name+stiw in wavesSig:
     1551                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
     1552                                else:
     1553                                    sigstr += 10*' '
     1554                        elif Stype == 'Smag':
     1555                            for j in range(6):
     1556                                name = names['Smag'][j]
     1557                                namstr += '%12s'%(names['Smag'][j])
     1558                                valstr += '%12.4f'%(wave[0][j])
     1559                                if name+stiw in wavesSig:
     1560                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
     1561                                else:
     1562                                    sigstr += 12*' '
     1563                               
     1564                    print >>pFile,namstr
     1565                    print >>pFile,valstr
     1566                    print >>pFile,sigstr
    15191567       
    15201568               
    … …  
    17841832                            if Stype == 'Spos':
    17851833                                if waveType in ['ZigZag','Sawtooth'] and not iw:
    1786                                     names = [pfx+'Tzero:'+stiw,pfx+'Xslope:'+stiw,pfx+'Yslope:'+stiw,pfx+'Zslope:'+stiw]
     1834                                    names = ['Tzero:'+stiw,'Xslope:'+stiw,'Yslope:'+stiw,'Zslope:'+stiw]
    17871835                                else:
    1788                                     names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
    1789                                         pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
     1836                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
     1837                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
    17901838                            elif Stype == 'Sadp':
    1791                                 names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
    1792                                     pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
    1793                                     pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
    1794                                     pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
     1839                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
     1840                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
     1841                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
     1842                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
    17951843                            elif Stype == 'Sfrac':
    17961844                                if 'Crenel' in waveType and not iw:
    1797                                     names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
     1845                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
    17981846                                else:
    1799                                     names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
     1847                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
    18001848                            elif Stype == 'Smag':
    1801                                 names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
    1802                                     pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
     1849                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
     1850                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
    18031851                            for iname,name in enumerate(names):
    1804                                 AtomSS[Stype][iw][0][iname] = parmDict[name]
    1805                                 if name in sigDict:
    1806                                     wavesSig[name] = sigDict[name]
     1852                                AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
     1853                                if pfx+name in sigDict:
     1854                                    wavesSig[name] = sigDict[pfx+name]
    18071855                   
    18081856            PrintAtomsAndSig(General,Atoms,atomsSig)

  • trunk/GSASIIstrMath.py

    r1630 r1635  
    711711        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
    712712        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
    713         fb = 0.
     713        fb = np.zeros_like(fa)
    714714        if not SGData['SGInv']:
    715715            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
    716         fa = fa*GfpuA.T-fb*GfpuB.T
    717         fb = fb*GfpuA.T+fa*GfpuB.T
    718         fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
    719         fbs = np.sum(np.sum(fb,axis=1),axis=1)
     716        fa = fa*GfpuA-fb*GfpuB
     717        fb = fb*GfpuA+fa*GfpuB
     718        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
     719        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
    720720           
    721721        fasq = fas**2
    … …  
    943943            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
    944944        H = np.array(refl[:4])
    945         if H[3]:
    946             continue
    947945        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
    948946        SQfactor = 8.0*SQ*np.pi**2
    … …  
    955953        Phi = np.inner(H[:3],SGT)
    956954        SSPhi = np.inner(H,SSGT)
    957         GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)       
     955        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
     956        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
     957        #need ModulationDerv here dGAdXsin, etc 
    958958        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
    959959        sinp = np.sin(phase)
    … …  
    962962        biso = -SQfactor*Uisodata
    963963        Tiso = np.where(biso<1.,np.exp(biso),1.0)
    964         HbH = -np.inner(H,np.inner(bij,H))
     964        HbH = -np.inner(H[:3],np.inner(bij,H[:3]))
    965965        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
    966966        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
    … …  
    971971        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
    972972        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
     973        GfpuA = np.swapaxes(GfpuA,1,2)
     974        GfpuB = np.swapaxes(GfpuB,1,2)
     975        fa = fa*GfpuA-fb*GfpuB
     976        fb = fb*GfpuA+fa*GfpuB
    973977       
    974978        fas = np.sum(np.sum(fa,axis=1),axis=1)
    … …  
    976980        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
    977981        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
     982        fax = fax*GfpuA-fbx*GfpuB
     983        fbx = fbx*GfpuA+fax*GfpuB
    978984        #sum below is over Uniq
    979985        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
    … …  
    988994        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
    989995        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     996        #need dFdXsin, etc for modulations...
    990997        if not SGData['SGInv']:
    991998            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
    … …  
    9991006            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
    10001007            dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1008        #need dFdXsin, etc for modulations...
    10011009        #loop over atoms - each dict entry is list of derivatives for all the reflections
    10021010    for i in range(len(Mdata)):     
    … …  
    10121020        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
    10131021        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1022        #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
    10141023    dFdvDict[pfx+'BabA'] = dFdbab.T[0]
    10151024    dFdvDict[pfx+'BabU'] = dFdbab.T[1]
Note: See TracChangeset for help on using the changeset viewer.