source: trunk/GSASIIstrMath.py @ 4213

Last change on this file since 4213 was 4213, checked in by toby, 2 years ago

multiple small changes to allow docs build under 3.x

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 218.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2019-12-20 05:04:03 +0000 (Fri, 20 Dec 2019) $
8# $Author: toby $
9# $Revision: 4213 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4213 2019-12-20 05:04:03Z toby $
12########### SVN repository information ###################
13from __future__ import division, print_function
14import time
15import copy
16import numpy as np
17import numpy.ma as ma
18import numpy.linalg as nl
19import scipy.stats as st
20import multiprocessing as mp
21import GSASIIpath
22GSASIIpath.SetVersionNumber("$Revision: 4213 $")
23import GSASIIElem as G2el
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import GSASIIpwd as G2pwd
27import GSASIImapvars as G2mv
28import GSASIImath as G2mth
29import GSASIIobj as G2obj
30import GSASIImpsubs as G2mp
31#G2mp.InitMP(False)  # This disables multiprocessing
32
33sind = lambda x: np.sin(x*np.pi/180.)
34cosd = lambda x: np.cos(x*np.pi/180.)
35tand = lambda x: np.tan(x*np.pi/180.)
36asind = lambda x: 180.*np.arcsin(x)/np.pi
37acosd = lambda x: 180.*np.arccos(x)/np.pi
38atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
39   
40try:  # fails on doc build
41    ateln2 = 8.0*np.log(2.0)
42    twopi = 2.0*np.pi
43    twopisq = 2.0*np.pi**2
44except TypeError:
45    pass
46nxs = np.newaxis
47
48################################################################################
49##### Rigid Body Models
50################################################################################
51       
52def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
53    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
54    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
55    '''
56    atxIds = ['Ax:','Ay:','Az:']
57    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
58    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
59    if not RBIds['Vector'] and not RBIds['Residue']:
60        return
61    VRBIds = RBIds['Vector']
62    RRBIds = RBIds['Residue']
63    if Update:
64        RBData = rigidbodyDict
65    else:
66        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
67    if RBIds['Vector']:                       # first update the vector magnitudes
68        VRBData = RBData['Vector']
69        for i,rbId in enumerate(VRBIds):
70            if VRBData[rbId]['useCount']:
71                for j in range(len(VRBData[rbId]['VectMag'])):
72                    name = '::RBV;'+str(j)+':'+str(i)
73                    VRBData[rbId]['VectMag'][j] = parmDict[name]
74    for phase in Phases:
75        Phase = Phases[phase]
76        General = Phase['General']
77        cx,ct,cs,cia = General['AtomPtrs']
78        cell = General['Cell'][1:7]
79        Amat,Bmat = G2lat.cell2AB(cell)
80        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
81        pfx = str(Phase['pId'])+'::'
82        if Update:
83            RBModels = Phase['RBModels']
84        else:
85            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
86        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
87            jrb = VRBIds.index(RBObj['RBId'])
88            rbsx = str(irb)+':'+str(jrb)
89            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
90                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
91            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
92                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
93            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
94            TLS = RBObj['ThermalMotion']
95            if 'T' in TLS[0]:
96                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
97                    TLS[1][i] = parmDict[pfx+pt+rbsx]
98            if 'L' in TLS[0]:
99                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
100                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
101            if 'S' in TLS[0]:
102                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
103                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
104            if 'U' in TLS[0]:
105                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
106            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
107            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
108            for i,x in enumerate(XYZ):
109                atId = RBObj['Ids'][i]
110                for j in [0,1,2]:
111                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
112                if UIJ[i][0] == 'A':
113                    for j in range(6):
114                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
115                elif UIJ[i][0] == 'I':
116                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
117           
118        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
119            jrb = RRBIds.index(RBObj['RBId'])
120            rbsx = str(irb)+':'+str(jrb)
121            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
122                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
123            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
124                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
125            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
126            TLS = RBObj['ThermalMotion']
127            if 'T' in TLS[0]:
128                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
129                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
130            if 'L' in TLS[0]:
131                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
132                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
133            if 'S' in TLS[0]:
134                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
135                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
136            if 'U' in TLS[0]:
137                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
138            for itors,tors in enumerate(RBObj['Torsions']):
139                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
140            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
141            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
142            for i,x in enumerate(XYZ):
143                atId = RBObj['Ids'][i]
144                for j in [0,1,2]:
145                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
146                if UIJ[i][0] == 'A':
147                    for j in range(6):
148                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
149                elif UIJ[i][0] == 'I':
150                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
151                   
152def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
153    'Needs a doc string'
154    atxIds = ['dAx:','dAy:','dAz:']
155    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
156    OIds = ['Oa:','Oi:','Oj:','Ok:']
157    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
158    if not RBIds['Vector'] and not RBIds['Residue']:
159        return
160    VRBIds = RBIds['Vector']
161    RRBIds = RBIds['Residue']
162    RBData = rigidbodyDict
163    for item in parmDict:
164        if 'RB' in item:
165            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
166    General = Phase['General']
167    cx,ct,cs,cia = General['AtomPtrs']
168    cell = General['Cell'][1:7]
169    Amat,Bmat = G2lat.cell2AB(cell)
170    rpd = np.pi/180.
171    rpd2 = rpd**2
172    g = nl.inv(np.inner(Bmat,Bmat))
173    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
174        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
175    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
176    pfx = str(Phase['pId'])+'::'
177    RBModels =  Phase['RBModels']
178    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
179        VModel = RBData['Vector'][RBObj['RBId']]
180        Q = RBObj['Orient'][0]
181        jrb = VRBIds.index(RBObj['RBId'])
182        rbsx = str(irb)+':'+str(jrb)
183        dXdv = []
184        for iv in range(len(VModel['VectMag'])):
185            dCdv = []
186            for vec in VModel['rbVect'][iv]:
187                dCdv.append(G2mth.prodQVQ(Q,vec))
188            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
189        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
190        for ia,atId in enumerate(RBObj['Ids']):
191            atNum = AtLookup[atId]
192            dx = 0.00001
193            for iv in range(len(VModel['VectMag'])):
194                for ix in [0,1,2]:
195                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
196            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
197                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
198            for iv in range(4):
199                Q[iv] -= dx
200                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
201                Q[iv] += 2.*dx
202                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
203                Q[iv] -= dx
204                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
205                for ix in [0,1,2]:
206                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
207            X = G2mth.prodQVQ(Q,Cart[ia])
208            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
209            dFdu = G2lat.U6toUij(dFdu.T)
210            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
211            dFdu = G2lat.UijtoU6(dFdu)
212            atNum = AtLookup[atId]
213            if 'T' in RBObj['ThermalMotion'][0]:
214                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
215                    dFdvDict[pfx+name+rbsx] += dFdu[i]
216            if 'L' in RBObj['ThermalMotion'][0]:
217                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
218                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
219                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
220                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
221                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
222                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
223                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
224                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
225                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
226            if 'S' in RBObj['ThermalMotion'][0]:
227                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
228                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
229                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
230                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
231                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
232                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
233                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
234                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
235            if 'U' in RBObj['ThermalMotion'][0]:
236                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
237
238
239    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
240        Q = RBObj['Orient'][0]
241        jrb = RRBIds.index(RBObj['RBId'])
242        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
243        rbsx = str(irb)+':'+str(jrb)
244        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
245        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
246            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
247            orId,pvId = torData[itors][:2]
248            pivotVec = Cart[orId]-Cart[pvId]
249            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
250            QB = G2mth.AVdeg2Q(0.001,pivotVec)
251            for ir in torData[itors][3]:
252                atNum = AtLookup[RBObj['Ids'][ir]]
253                rVec = Cart[ir]-Cart[pvId]
254                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
255                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
256                for ix in [0,1,2]:
257                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
258        for ia,atId in enumerate(RBObj['Ids']):
259            atNum = AtLookup[atId]
260            dx = 0.00001
261            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
262                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
263            for iv in range(4):
264                Q[iv] -= dx
265                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
266                Q[iv] += 2.*dx
267                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
268                Q[iv] -= dx
269                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
270                for ix in [0,1,2]:
271                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
272            X = G2mth.prodQVQ(Q,Cart[ia])
273            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
274            dFdu = G2lat.U6toUij(dFdu.T)
275            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
276            dFdu = G2lat.UijtoU6(dFdu)
277            atNum = AtLookup[atId]
278            if 'T' in RBObj['ThermalMotion'][0]:
279                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
280                    dFdvDict[pfx+name+rbsx] += dFdu[i]
281            if 'L' in RBObj['ThermalMotion'][0]:
282                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
283                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
284                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
285                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
286                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
287                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
288                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
289                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
290                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
291            if 'S' in RBObj['ThermalMotion'][0]:
292                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
293                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
294                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
295                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
296                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
297                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
298                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
299                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
300            if 'U' in RBObj['ThermalMotion'][0]:
301                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
302   
303################################################################################
304##### Penalty & restraint functions
305################################################################################
306
307def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
308    'Compute user-supplied and built-in restraint functions'
309    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
310    pNames = []
311    pVals = []
312    pWt = []
313    negWt = {}
314    pWsum = {}
315    pWnum = {}
316    for phase in Phases:
317        pId = Phases[phase]['pId']
318        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
319        General = Phases[phase]['General']
320        cx,ct,cs,cia = General['AtomPtrs']
321        textureData = General['SH Texture']
322        SGData = General['SGData']
323        Atoms = Phases[phase]['Atoms']
324        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
325        cell = General['Cell'][1:7]
326        Amat,Bmat = G2lat.cell2AB(cell)
327        if phase not in restraintDict:
328            continue
329        phaseRest = restraintDict[phase]
330        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
331            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
332            ['ChemComp','Sites'],['Texture','HKLs'],['General','General'],]
333        for name,rest in names:
334            pWsum[name] = 0.
335            pWnum[name] = 0
336            if name not in phaseRest:
337                continue
338            itemRest = phaseRest[name]
339            if itemRest[rest] and itemRest['Use']:
340                wt = itemRest.get('wtFactor',1.)
341                if name in ['Bond','Angle','Plane','Chiral']:
342                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
343                        pNames.append(str(pId)+':'+name+':'+str(i))
344                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
345                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
346                        if name == 'Bond':
347                            calc = G2mth.getRestDist(XYZ,Amat)
348                        elif name == 'Angle':
349                            calc = G2mth.getRestAngle(XYZ,Amat)
350                        elif name == 'Plane':
351                            calc = G2mth.getRestPlane(XYZ,Amat)
352                        elif name == 'Chiral':
353                            calc = G2mth.getRestChiral(XYZ,Amat)
354                        pVals.append(obs-calc)
355                        pWt.append(wt/esd**2)
356                        pWsum[name] += wt*((obs-calc)/esd)**2
357                        pWnum[name] += 1
358                elif name in ['Torsion','Rama']:
359                    coeffDict = itemRest['Coeff']
360                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
361                        pNames.append(str(pId)+':'+name+':'+str(i))
362                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
363                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
364                        if name == 'Torsion':
365                            tor = G2mth.getRestTorsion(XYZ,Amat)
366                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
367                        else:
368                            phi,psi = G2mth.getRestRama(XYZ,Amat)
369                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
370                        pVals.append(restr)
371                        pWt.append(wt/esd**2)
372                        pWsum[name] += wt*(restr/esd)**2
373                        pWnum[name] += 1
374                elif name == 'ChemComp':
375                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
376                        pNames.append(str(pId)+':'+name+':'+str(i))
377                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
378                        frac = np.array(G2mth.GetAtomFracByID(pId,parmDict,AtLookup,indx))
379                        calc = np.sum(mul*frac*factors)
380                        pVals.append(obs-calc)
381                        pWt.append(wt/esd**2)                   
382                        pWsum[name] += wt*((obs-calc)/esd)**2
383                        pWnum[name] += 1
384                elif name == 'Texture':
385                    SHkeys = list(textureData['SH Coeff'][1].keys())
386                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
387                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
388                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
389                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
390                        PH = np.array(hkl)
391                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
392                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
393                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
394                        Z1 = ma.masked_greater(Z,0.0)           #is this + or -?
395                        IndZ1 = np.array(ma.nonzero(Z1))
396                        for ind in IndZ1.T:
397                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
398                            pVals.append(Z1[ind[0]][ind[1]])
399                            pWt.append(wt/esd1**2)
400                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
401                            pWnum[name] += 1
402                        if ifesd2:
403                            Z2 = 1.-Z
404                            for ind in np.ndindex(grid,grid):
405                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
406                                pVals.append(Z2[ind[0]][ind[1]])
407                                pWt.append(wt/esd2**2)
408                                pWsum[name] += wt*(Z2/esd2)**2
409                                pWnum[name] += 1
410                elif name == 'General':
411                    for i,(eq,obs,esd) in enumerate(itemRest[rest]):
412                        pNames.append(str(pId)+':'+name+':'+str(i))
413                        calcobj = G2obj.ExpressionCalcObj(eq)
414                        calcobj.SetupCalc(parmDict)
415                        calc = calcobj.EvalExpression()
416                        pVals.append(obs-calc)
417                        pWt.append(wt/esd**2)                   
418                        pWsum[name] += wt*((obs-calc)/esd)**2
419                        pWnum[name] += 1
420       
421    for phase in Phases:
422        name = 'SH-Pref.Ori.'
423        pId = Phases[phase]['pId']
424        General = Phases[phase]['General']
425        SGData = General['SGData']
426        cell = General['Cell'][1:7]
427        pWsum[name] = 0.0
428        pWnum[name] = 0
429        for hist in Phases[phase]['Histograms']:
430            if not Phases[phase]['Histograms'][hist]['Use']:
431                continue
432            if hist in Histograms and 'PWDR' in hist:
433                hId = Histograms[hist]['hId']
434                phfx = '%d:%d:'%(pId,hId)
435                if calcControls[phfx+'poType'] == 'SH':
436                    toler = calcControls[phfx+'SHtoler']
437                    wt = 1./toler**2
438                    HKLs = np.array(calcControls[phfx+'SHhkl'])
439                    SHnames = calcControls[phfx+'SHnames']
440                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
441                    for i,PH in enumerate(HKLs):
442                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
443                        SH3Coef = {}
444                        for item in SHcof:
445                            L,N = eval(item.strip('C'))
446                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
447                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
448                        X = np.linspace(0,90.0,26)
449                        Y = ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)       #+ or -?
450                        IndY = ma.nonzero(Y)
451                        for ind in IndY[0]:
452                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
453                            pVals.append(Y[ind])
454                            pWt.append(wt)
455                            pWsum[name] += wt*(Y[ind])**2
456                            pWnum[name] += 1
457    pWsum['PWLref'] = 0.
458    pWnum['PWLref'] = 0
459    for item in varyList:
460        if 'PWLref' in item and parmDict[item] < 0.:
461            pId = int(item.split(':')[0])
462            if negWt[pId]:
463                pNames.append(item)
464                pVals.append(parmDict[item])
465                pWt.append(negWt[pId])
466                pWsum['PWLref'] += negWt[pId]*(parmDict[item])**2
467                pWnum['PWLref'] += 1
468    pVals = np.array(pVals)
469    pWt = np.array(pWt)         #should this be np.sqrt?
470    return pNames,pVals,pWt,pWsum,pWnum
471   
472def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
473    '''Compute derivatives on user-supplied and built-in restraint
474    (penalty) functions
475
476    where pNames is list of restraint labels
477
478    returns pDerv with partial derivatives by variable# in varList and
479       restraint# in pNames (pDerv[variable#][restraint#])
480    '''
481    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
482    pDerv = np.zeros((len(varyList),len(pVal)))
483    for pName in pNames: # loop over restraints
484        if 'General' == pName.split(':')[1]:
485            # initialize for General restraint(s) here
486            GeneralInit = True
487            parmDict0 = parmDict.copy()
488            # setup steps for each parameter
489            stepDict = {}
490            for parm in varyList:
491                stepDict[parm] = G2obj.getVarStep(parm,parmDict)
492            break
493    for phase in Phases:
494#        if phase not in restraintDict:
495#            continue
496        pId = Phases[phase]['pId']
497        General = Phases[phase]['General']
498        cx,ct,cs,cia = General['AtomPtrs']
499        SGData = General['SGData']
500        Atoms = Phases[phase]['Atoms']
501        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
502        cell = General['Cell'][1:7]
503        Amat,Bmat = G2lat.cell2AB(cell)
504        textureData = General['SH Texture']
505
506        SHkeys = list(textureData['SH Coeff'][1].keys())
507        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
508        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
509        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
510        sam = SamSym[textureData['Model']]
511        phaseRest = restraintDict.get(phase,{})
512        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
513            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
514            'ChemComp':'Sites','Texture':'HKLs'}
515        lasthkl = np.array([0,0,0])
516        for ip,pName in enumerate(pNames): # loop over restraints
517            pnames = pName.split(':')
518            if pId == int(pnames[0]):
519                name = pnames[1]
520                if 'PWL' in pName:
521                    pDerv[varyList.index(pName)][ip] += 1.
522                    continue
523                elif 'SH-' in pName:
524                    continue
525                Id = int(pnames[2]) 
526                itemRest = phaseRest[name]
527                if name in ['Bond','Angle','Plane','Chiral']:
528                    indx,ops,obs,esd = itemRest[names[name]][Id]
529                    dNames = []
530                    for ind in indx:
531                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
532                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
533                    if name == 'Bond':
534                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
535                    elif name == 'Angle':
536                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
537                    elif name == 'Plane':
538                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
539                    elif name == 'Chiral':
540                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
541                elif name in ['Torsion','Rama']:
542                    coffDict = itemRest['Coeff']
543                    indx,ops,cofName,esd = itemRest[names[name]][Id]
544                    dNames = []
545                    for ind in indx:
546                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
547                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
548                    if name == 'Torsion':
549                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
550                    else:
551                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
552                elif name == 'ChemComp':
553                    indx,factors,obs,esd = itemRest[names[name]][Id]
554                    dNames = []
555                    for ind in indx:
556                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
557                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
558                        deriv = mul*factors
559                elif 'Texture' in name:
560                    deriv = []
561                    dNames = []
562                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][Id]
563                    hkl = np.array(hkl)
564                    if np.any(lasthkl-hkl):
565                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
566                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
567                        lasthkl = copy.copy(hkl)                       
568                    if 'unit' in name:
569                        pass
570                    else:
571                        gam = float(pnames[3])
572                        psi = float(pnames[4])
573                        for SHname in ODFln:
574                            l,m,n = eval(SHname[1:])
575                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
576                            dNames += [str(pId)+'::'+SHname]
577                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
578                elif name == 'General':
579                    deriv = []
580                    dNames = []
581                    eq,obs,esd = itemRest[name][Id]
582                    calcobj = G2obj.ExpressionCalcObj(eq)
583                    parmlist = list(eq.assgnVars.values()) # parameters used in this expression
584                    for parm in parmlist: # expand list if any parms are determined by constraints
585                        if parm in G2mv.dependentVars:
586                            parmlist += G2mv.independentVars
587                            break
588                    for ind,var in enumerate(varyList):
589                        drv = 0
590                        if var in parmlist:
591                            step = stepDict.get(var,1e-5)
592                            calc = []
593                            # apply step to parameter
594                            oneparm = True
595                            for s in -step,2*step:
596                                parmDict[var] += s
597                                # extend shift if needed to other parameters
598                                if var in G2mv.independentVars:
599                                    G2mv.Dict2Map(parmDict,[])
600                                    oneparm = False
601                                elif var in G2mv.dependentVars:
602                                    G2mv.Map2Dict(parmDict,[])
603                                    oneparm = False
604                                if 'RB' in var:
605                                    ApplyRBModels(parmDict,Phases,rigidbodyDict)
606# test
607                                    oneparm = False
608                                calcobj.SetupCalc(parmDict)
609                                calc.append(calcobj.EvalExpression())
610                            drv = (calc[1]-calc[0])*.5/step
611                            # restore the dict
612                            if oneparm:
613                                parmDict[var] = parmDict0[var]
614                            else:
615                                parmDict = parmDict0.copy()
616                        else:
617                            drv = 0
618                        pDerv[ind][ip] = drv
619                # Add derivatives into matrix, if needed
620                for dName,drv in zip(dNames,deriv):
621                    try:
622                        ind = varyList.index(dName)
623                        pDerv[ind][ip] += drv
624                    except ValueError:
625                        pass
626       
627        lasthkl = np.array([0,0,0])
628        for ip,pName in enumerate(pNames):
629            deriv = []
630            dNames = []
631            pnames = pName.split(':')
632            if 'SH-' in pName and pId == int(pnames[0]):
633                hId = int(pnames[1])
634                phfx = '%d:%d:'%(pId,hId)
635                psi = float(pnames[4])
636                HKLs = calcControls[phfx+'SHhkl']
637                SHnames = calcControls[phfx+'SHnames']
638                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
639                hkl = np.array(HKLs[int(pnames[3])])     
640                if np.any(lasthkl-hkl):
641                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
642                    SH3Coef = {}
643                    for item in SHcof:
644                        L,N = eval(item.strip('C'))
645                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
646                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
647                    lasthkl = copy.copy(hkl)                       
648                for SHname in SHnames:
649                    l,n = eval(SHname[1:])
650                    SH3name = 'C%d,0,%d'%(l,n)
651                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
652                    dNames += [phfx+SHname]
653                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
654            for dName,drv in zip(dNames,deriv):
655                try:
656                    ind = varyList.index(dName)
657                    pDerv[ind][ip] += drv
658                except ValueError:
659                    pass
660    return pDerv
661
662################################################################################
663##### Function & derivative calculations
664################################################################################       
665                   
666def GetAtomFXU(pfx,calcControls,parmDict):
667    'Needs a doc string'
668    Natoms = calcControls['Natoms'][pfx]
669    Tdata = Natoms*[' ',]
670    Mdata = np.zeros(Natoms)
671    IAdata = Natoms*[' ',]
672    Fdata = np.zeros(Natoms)
673    Xdata = np.zeros((3,Natoms))
674    dXdata = np.zeros((3,Natoms))
675    Uisodata = np.zeros(Natoms)
676    Uijdata = np.zeros((6,Natoms))
677    Gdata = np.zeros((3,Natoms))
678    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
679        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
680        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
681        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
682        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5],
683        'AMx:':Gdata[0],'AMy:':Gdata[1],'AMz:':Gdata[2],}
684    for iatm in range(Natoms):
685        for key in keys:
686            parm = pfx+key+str(iatm)
687            if parm in parmDict:
688                keys[key][iatm] = parmDict[parm]
689    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
690    Gdata = np.where(Gdata,Gdata,1.e-8)         #avoid divide by zero in derivative calc.
691   
692    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata
693   
694def GetAtomSSFXU(pfx,calcControls,parmDict):
695    'Needs a doc string'
696    Natoms = calcControls['Natoms'][pfx]
697    maxSSwave = calcControls['maxSSwave'][pfx]
698    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
699        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
700    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
701    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
702    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
703    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
704    waveTypes = []
705    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
706        'Tmin:':XSSdata[0],'Tmax:':XSSdata[1],'Xmax:':XSSdata[2],'Ymax:':XSSdata[3],'Zmax:':XSSdata[4],
707        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
708        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
709        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
710        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
711    for iatm in range(Natoms):
712        wavetype = [parmDict.get(pfx+kind+'waveType:'+str(iatm),'') for kind in ['F','P','A','M']]
713        waveTypes.append(wavetype)
714        for key in keys:
715            for m in range(Nwave[key[0]]):
716                parm = pfx+key+str(iatm)+':%d'%(m)
717                if parm in parmDict:
718                    keys[key][m][iatm] = parmDict[parm]
719    return waveTypes,FSSdata,XSSdata,USSdata,MSSdata
720   
721def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
722    ''' Compute structure factors for all h,k,l for phase
723    puts the result, F^2, in each ref[8] in refList
724    operates on blocks of 100 reflections for speed
725    input:
726   
727    :param dict refDict: where
728        'RefList' list where each ref = h,k,l,it,d,...
729        'FF' dict of form factors - filed in below
730    :param np.array G:      reciprocal metric tensor
731    :param str pfx:    phase id string
732    :param dict SGData: space group info. dictionary output from SpcGroup
733    :param dict calcControls:
734    :param dict ParmDict:
735
736    '''       
737    phfx = pfx.split(':')[0]+hfx
738    ast = np.sqrt(np.diag(G))
739    Mast = twopisq*np.multiply.outer(ast,ast)
740    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
741    SGT = np.array([ops[1] for ops in SGData['SGOps']])
742    FFtables = calcControls['FFtables']
743    BLtables = calcControls['BLtables']
744    Amat,Bmat = G2lat.Gmat2AB(G)
745    Flack = 1.0
746    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
747        Flack = 1.-2.*parmDict[phfx+'Flack']
748    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
749    TwDict = refDict.get('TwDict',{})           
750    if 'S' in calcControls[hfx+'histType']:
751        NTL = calcControls[phfx+'NTL']
752        NM = calcControls[phfx+'TwinNMN']+1
753        TwinLaw = calcControls[phfx+'TwinLaw']
754        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
755        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
756    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
757        GetAtomFXU(pfx,calcControls,parmDict)
758    if not Xdata.size:          #no atoms in phase!
759        return
760    if 'NC' in calcControls[hfx+'histType']:
761        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
762    elif 'X' in calcControls[hfx+'histType']:
763        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
764        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
765    Uij = np.array(G2lat.U6toUij(Uijdata))
766    bij = Mast*Uij.T
767    blkSize = 100       #no. of reflections in a block - size seems optimal
768    nRef = refDict['RefList'].shape[0]
769    SQ = 1./(2.*refDict['RefList'].T[4])**2
770    if 'N' in calcControls[hfx+'histType']:
771        dat = G2el.getBLvalues(BLtables)
772        refDict['FF']['El'] = list(dat.keys())
773        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
774    else:       #'X'
775        dat = G2el.getFFvalues(FFtables,0.)
776        refDict['FF']['El'] = list(dat.keys())
777        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
778        for iel,El in enumerate(refDict['FF']['El']):
779            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
780#reflection processing begins here - big arrays!
781    iBeg = 0
782    while iBeg < nRef:
783        iFin = min(iBeg+blkSize,nRef)
784        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
785        H = refl.T[:3]                          #array(blkSize,3)
786        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
787        TwMask = np.any(H,axis=-1)
788        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
789            for ir in range(blkSize):
790                iref = ir+iBeg
791                if iref in TwDict:
792                    for i in TwDict[iref]:
793                        for n in range(NTL):
794                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
795            TwMask = np.any(H,axis=-1)
796        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
797        SQfactor = 4.0*SQ*twopisq               #ditto prev.
798        if 'T' in calcControls[hfx+'histType']:
799            if 'P' in calcControls[hfx+'histType']:
800                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
801            else:
802                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
803            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
804            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
805        Uniq = np.inner(H,SGMT)
806        Phi = np.inner(H,SGT)
807        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
808        sinp = np.sin(phase)
809        cosp = np.cos(phase)
810        biso = -SQfactor*Uisodata[:,nxs]
811        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
812        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
813        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
814        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
815        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
816        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
817        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
818        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
819            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
820            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
821        else:
822            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
823            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
824        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
825        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
826        if SGData['SGInv']: #centrosymmetric; B=0
827            fbs[0] *= 0.
828            fas[1] *= 0.
829        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
830            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
831            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
832        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
833            if len(TwinLaw) > 1:
834                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
835                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
836                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
837                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
838            else:   # checked correct!!
839                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
840                refl.T[7] = np.copy(refl.T[9])               
841                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
842#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
843        iBeg += blkSize
844#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
845   
846def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
847    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
848    faster than StructureFactorDerv - correct for powders/nontwins!!
849    input:
850   
851    :param dict refDict: where
852        'RefList' list where each ref = h,k,l,it,d,...
853        'FF' dict of form factors - filled in below
854    :param np.array G:      reciprocal metric tensor
855    :param str hfx:    histogram id string
856    :param str pfx:    phase id string
857    :param dict SGData: space group info. dictionary output from SpcGroup
858    :param dict calcControls:
859    :param dict parmDict:
860   
861    :returns: dict dFdvDict: dictionary of derivatives
862    '''
863    phfx = pfx.split(':')[0]+hfx
864    ast = np.sqrt(np.diag(G))
865    Mast = twopisq*np.multiply.outer(ast,ast)
866    SGMT = np.array([ops[0] for ops in SGData['SGOps']])
867    SGT = np.array([ops[1] for ops in SGData['SGOps']])
868    FFtables = calcControls['FFtables']
869    BLtables = calcControls['BLtables']
870    Amat,Bmat = G2lat.Gmat2AB(G)
871    nRef = len(refDict['RefList'])
872    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
873        GetAtomFXU(pfx,calcControls,parmDict)
874    if not Xdata.size:          #no atoms in phase!
875        return {}
876    mSize = len(Mdata)
877    FF = np.zeros(len(Tdata))
878    if 'NC' in calcControls[hfx+'histType']:
879        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
880    elif 'X' in calcControls[hfx+'histType']:
881        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
882        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
883    Uij = np.array(G2lat.U6toUij(Uijdata))
884    bij = Mast*Uij.T
885    dFdvDict = {}
886    dFdfr = np.zeros((nRef,mSize))
887    dFdx = np.zeros((nRef,mSize,3))
888    dFdui = np.zeros((nRef,mSize))
889    dFdua = np.zeros((nRef,mSize,6))
890    dFdbab = np.zeros((nRef,2))
891    dFdfl = np.zeros((nRef))
892    Flack = 1.0
893    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
894        Flack = 1.-2.*parmDict[phfx+'Flack']
895#    time0 = time.time()
896#reflection processing begins here - big arrays!
897    iBeg = 0
898    blkSize = 32       #no. of reflections in a block - optimized for speed
899    while iBeg < nRef:
900        iFin = min(iBeg+blkSize,nRef)
901        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
902        H = refl.T[:3].T
903        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
904        SQfactor = 8.0*SQ*np.pi**2
905        if 'T' in calcControls[hfx+'histType']:
906            if 'P' in calcControls[hfx+'histType']:
907                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
908            else:
909                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
910            FP = np.repeat(FP.T,len(SGT),axis=0)
911            FPP = np.repeat(FPP.T,len(SGT),axis=0)
912        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
913        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
914        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
915        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
916        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
917        Phi = np.inner(H,SGT)
918        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
919        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
920        cosp = np.cos(phase)
921        occ = Mdata*Fdata/len(SGT)
922        biso = -SQfactor*Uisodata[:,nxs]
923        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
924        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
925        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
926        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
927        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
928        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
929        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
930        if len(FPP.shape) > 1:
931            fotp = np.reshape(FPP,cosp.shape)*Tcorr
932        else:
933            fotp = FPP*Tcorr     
934        if 'T' in calcControls[hfx+'histType']:
935            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
936            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
937        else:
938            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
939            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
940        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
941        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
942        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,nEqv,nAtoms)
943        fbx = np.array([fot*cosp,-fotp*sinp])
944        #sum below is over Uniq
945        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
946        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
947        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
948        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
949        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
950        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
951        if not SGData['SGInv']:
952            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
953            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
954            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
955            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
956            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
957            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
958            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
959        else:
960            dfbdfr = np.zeros_like(dfadfr)
961            dfbdx = np.zeros_like(dfadx)
962            dfbdui = np.zeros_like(dfadui)
963            dfbdua = np.zeros_like(dfadua)
964            dfbdba = np.zeros_like(dfadba)
965            dfadfl = 0.0
966            dfbdfl = 0.0
967        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
968        SA = fas[0]+fas[1]
969        SB = fbs[0]+fbs[1]
970        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
971            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
972            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
973            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
974            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
975        else:
976            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
977            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
978            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
979            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
980            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
981        dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \
982                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
983        iBeg += blkSize
984#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
985        #loop over atoms - each dict entry is list of derivatives for all the reflections
986    for i in range(len(Mdata)):
987        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
988        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
989        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
990        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
991        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
992        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
993        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
994        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
995        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
996        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
997        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
998    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
999    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1000    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1001    return dFdvDict
1002   
1003def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1004    ''' Compute neutron magnetic structure factors for all h,k,l for phase
1005    puts the result, F^2, in each ref[8] in refList
1006    operates on blocks of 100 reflections for speed
1007    input:
1008   
1009    :param dict refDict: where
1010        'RefList' list where each ref = h,k,l,it,d,...
1011        'FF' dict of form factors - filed in below
1012    :param np.array G:      reciprocal metric tensor
1013    :param str pfx:    phase id string
1014    :param dict SGData: space group info. dictionary output from SpcGroup
1015    :param dict calcControls:
1016    :param dict ParmDict:
1017       
1018    :returns: copy of new refList - used in calculating numerical derivatives
1019
1020    '''       
1021    g = nl.inv(G)
1022    ast = np.sqrt(np.diag(G))
1023    ainv = np.sqrt(np.diag(g))
1024    GS = G/np.outer(ast,ast)
1025    Ginv = g/np.outer(ainv,ainv)
1026    uAmat = G2lat.Gmat2AB(GS)[0]
1027    Mast = twopisq*np.multiply.outer(ast,ast)
1028    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1029    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1030    Ncen = len(SGData['SGCen'])
1031    Nops = len(SGMT)*Ncen
1032    if not SGData['SGFixed']:
1033        Nops *= (1+SGData['SGInv'])
1034    MFtables = calcControls['MFtables']
1035    Bmat = G2lat.Gmat2AB(G)[1]
1036    TwinLaw = np.ones(1)
1037#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
1038#    TwDict = refDict.get('TwDict',{})           
1039#    if 'S' in calcControls[hfx+'histType']:
1040#        NTL = calcControls[phfx+'NTL']
1041#        NM = calcControls[phfx+'TwinNMN']+1
1042#        TwinLaw = calcControls[phfx+'TwinLaw']
1043#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1044#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1045    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1046        GetAtomFXU(pfx,calcControls,parmDict)
1047    if not Xdata.size:          #no atoms in phase!
1048        return
1049    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
1050    Gdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
1051    if SGData['SGInv'] and not SGData['SGFixed']:
1052        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1053    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering--> [Mxyz,nops,natms]
1054    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
1055    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1056    Kdata = np.inner(Gdata.T,uAmat).T
1057    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1058    Kdata /= Kmean     #Cartesian unit vectors
1059    Uij = np.array(G2lat.U6toUij(Uijdata))
1060    bij = Mast*Uij.T
1061    blkSize = 100       #no. of reflections in a block - size seems optimal
1062    nRef = refDict['RefList'].shape[0]
1063    SQ = 1./(2.*refDict['RefList'].T[4])**2
1064    refDict['FF']['El'] = list(MFtables.keys())
1065    refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
1066    for iel,El in enumerate(refDict['FF']['El']):
1067        refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1068#reflection processing begins here - big arrays!
1069    iBeg = 0
1070    while iBeg < nRef:
1071        iFin = min(iBeg+blkSize,nRef)
1072        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1073        H = refl.T[:3].T                          #array(blkSize,3)
1074#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
1075#        TwMask = np.any(H,axis=-1)
1076#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1077#            for ir in range(blkSize):
1078#                iref = ir+iBeg
1079#                if iref in TwDict:
1080#                    for i in TwDict[iref]:
1081#                        for n in range(NTL):
1082#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1083#            TwMask = np.any(H,axis=-1)
1084        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
1085        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1086        Uniq = np.inner(H,SGMT)
1087        Phi = np.inner(H,SGT)
1088        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1089        biso = -SQfactor*Uisodata[:,nxs]
1090        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
1091        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1092        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1093        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1094        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1095        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1096        if SGData['SGInv']:
1097            if not SGData['SGFixed']:
1098                mphase = np.hstack((phase,-phase))  #OK
1099            else:
1100                mphase = phase
1101        else:
1102            mphase = phase                    #
1103        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1104        mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
1105        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1106        cosm = np.cos(mphase)                               #ditto
1107        HM = np.inner(Bmat,H)                             #put into cartesian space
1108        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Kdata = MAGS & HM = UVEC in magstrfc.for both OK
1109        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
1110        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1111        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1112        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1113        fams = np.sum(np.sum(fam,axis=-1),axis=-1)     #Mxyz,Nref  Sum(sum(fam,atoms),ops)
1114        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)     #ditto
1115        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)   #Sum(fams**2,Mxyz) Re + Im
1116        refl.T[7] = np.copy(refl.T[9])               
1117        refl.T[10] = atan2d(fbms[0],fams[0]) #- what is phase for mag refl?
1118#        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
1119#            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
1120#            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1121#        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
1122#            if len(TwinLaw) > 1:
1123#                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
1124#                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
1125#                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
1126#                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
1127#            else:   # checked correct!!
1128#                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
1129#                refl.T[7] = np.copy(refl.T[9])               
1130#                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1131##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
1132        iBeg += blkSize
1133#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
1134    return copy.deepcopy(refDict['RefList'])
1135
1136def MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1137    '''Compute magnetic structure factor derivatives numerically - for powders/nontwins only
1138    input:
1139   
1140    :param dict refDict: where
1141        'RefList' list where each ref = h,k,l,it,d,...
1142        'FF' dict of form factors - filled in below
1143    :param np.array G:      reciprocal metric tensor
1144    :param str hfx:    histogram id string
1145    :param str pfx:    phase id string
1146    :param dict SGData: space group info. dictionary output from SpcGroup
1147    :param dict calcControls:
1148    :param dict parmDict:
1149   
1150    :returns: dict dFdvDict: dictionary of magnetic derivatives
1151    '''
1152   
1153    trefDict = copy.deepcopy(refDict)
1154    dM = 1.e-6
1155    dFdvDict = {}
1156    for parm in parmDict:
1157        if 'AM' in parm:
1158            parmDict[parm] += dM
1159            prefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
1160            parmDict[parm] -= 2*dM
1161            mrefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
1162            parmDict[parm] += dM
1163            dFdvDict[parm] = (prefList[:,9]-mrefList[:,9])/(2.*dM)
1164    return dFdvDict
1165           
1166def MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1167    '''Compute nonmagnetic structure factor derivatives on blocks of reflections in magnetic structures - for powders/nontwins only
1168    input:
1169   
1170    :param dict refDict: where
1171        'RefList' list where each ref = h,k,l,it,d,...
1172        'FF' dict of form factors - filled in below
1173    :param np.array G:      reciprocal metric tensor
1174    :param str hfx:    histogram id string
1175    :param str pfx:    phase id string
1176    :param dict SGData: space group info. dictionary output from SpcGroup
1177    :param dict calcControls:
1178    :param dict parmDict:
1179   
1180    :returns: dict dFdvDict: dictionary of derivatives
1181    '''
1182   
1183    g = nl.inv(G)
1184    ast = np.sqrt(np.diag(G))
1185    ainv = np.sqrt(np.diag(g))
1186    GS = G/np.outer(ast,ast)
1187    Ginv = g/np.outer(ainv,ainv)
1188    uAmat = G2lat.Gmat2AB(GS)[0]
1189    Mast = twopisq*np.multiply.outer(ast,ast)
1190    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1191    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1192    Ncen = len(SGData['SGCen'])
1193    Nops = len(SGMT)*Ncen
1194    if not SGData['SGFixed']:
1195        Nops *= (1+SGData['SGInv'])
1196    Bmat = G2lat.Gmat2AB(G)[1]
1197    nRef = len(refDict['RefList'])
1198    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1199        GetAtomFXU(pfx,calcControls,parmDict)
1200    if not Xdata.size:          #no atoms in phase!
1201        return {}
1202    mSize = len(Mdata)
1203    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
1204    Gones = np.ones_like(Gdata)
1205    Gdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
1206    Gones = np.inner(Gones.T,SGMT).T
1207    if SGData['SGInv'] and not SGData['SGFixed']:
1208        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1209        Gones = np.hstack((Gones,-Gones))       #inversion if any
1210    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
1211    Gones = np.hstack([Gones for icen in range(Ncen)])        #dup over cell centering
1212    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
1213    Gones = SGData['MagMom'][nxs,:,nxs]*Gones   #flip vectors according to spin flip
1214    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1215    Kdata = np.inner(Gdata.T,uAmat).T     #Cartesian unit vectors
1216    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1217    Kdata /= Kmean
1218    Uij = np.array(G2lat.U6toUij(Uijdata))
1219    bij = Mast*Uij.T
1220    dFdvDict = {}
1221    dFdfr = np.zeros((nRef,mSize))
1222    dFdx = np.zeros((nRef,mSize,3))
1223    dFdui = np.zeros((nRef,mSize))
1224    dFdua = np.zeros((nRef,mSize,6))
1225    time0 = time.time()
1226#reflection processing begins here - big arrays!
1227    iBeg = 0
1228    blkSize = 5       #no. of reflections in a block - optimized for speed
1229    while iBeg < nRef:
1230        iFin = min(iBeg+blkSize,nRef)
1231        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1232        H = refl.T[:3].T
1233        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1234        SQfactor = 8.0*SQ*np.pi**2
1235        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
1236        Phi = np.inner(H,SGT)
1237        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1238        occ = Mdata*Fdata/Nops
1239        biso = -SQfactor*Uisodata[:,nxs]
1240        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
1241        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1242        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1243        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1244        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))
1245        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1246        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1247        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1248        if SGData['SGInv']:
1249            if not SGData['SGFixed']:
1250                mphase = np.hstack((phase,-phase))  #OK
1251                Uniq = np.hstack((Uniq,-Uniq))      #Nref,Nops,hkl
1252                Hij = np.hstack((Hij,Hij))
1253            else:
1254                mphase = phase
1255        else:
1256            mphase = phase                    #
1257        Hij = np.concatenate(np.array([Hij for cen in SGData['SGCen']]),axis=1)
1258        Uniq = np.hstack([Uniq for cen in SGData['SGCen']])
1259        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1260        mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
1261        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1262        cosm = np.cos(mphase)                               #ditto
1263        HM = np.inner(Bmat.T,H)                             #put into cartesian space
1264        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit cartesian vector for H
1265        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
1266        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1267       
1268        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
1269        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]
1270        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                      #Mxyz,Nref
1271        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)
1272        famx = -Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*sinm[nxs,:,:,:]   #Mxyz,Nref,Nops,Natom
1273        fbmx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*cosm[nxs,:,:,:]
1274        #sums below are over Nops - real part
1275        dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK
1276        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK
1277        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
1278        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK
1279        # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
1280        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
1281        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2)
1282        dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
1283        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
1284        #accumulate derivatives   
1285        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
1286        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
1287        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
1288        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok
1289        iBeg += blkSize
1290    print (' %d derivative time %.4f\r'%(nRef,time.time()-time0))
1291        #loop over atoms - each dict entry is list of derivatives for all the reflections
1292    for i in range(len(Mdata)):
1293        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1294        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1295        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1296        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1297        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1298        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1299        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1300        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1301        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1302        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1303        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1304    return dFdvDict
1305       
1306def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1307    '''Compute structure factor derivatives on blocks of reflections - for twins only
1308    faster than StructureFactorDervTw
1309    input:
1310   
1311    :param dict refDict: where
1312        'RefList' list where each ref = h,k,l,it,d,...
1313        'FF' dict of form factors - filled in below
1314    :param np.array G:      reciprocal metric tensor
1315    :param str hfx:    histogram id string
1316    :param str pfx:    phase id string
1317    :param dict SGData: space group info. dictionary output from SpcGroup
1318    :param dict calcControls:
1319    :param dict parmDict:
1320   
1321    :returns: dict dFdvDict: dictionary of derivatives
1322    '''
1323    phfx = pfx.split(':')[0]+hfx
1324    ast = np.sqrt(np.diag(G))
1325    Mast = twopisq*np.multiply.outer(ast,ast)
1326    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1327    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1328    FFtables = calcControls['FFtables']
1329    BLtables = calcControls['BLtables']
1330    TwDict = refDict.get('TwDict',{})           
1331    NTL = calcControls[phfx+'NTL']
1332    NM = calcControls[phfx+'TwinNMN']+1
1333    TwinLaw = calcControls[phfx+'TwinLaw']
1334    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1335    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1336    nTwin = len(TwinLaw)       
1337    nRef = len(refDict['RefList'])
1338    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1339        GetAtomFXU(pfx,calcControls,parmDict)
1340    if not Xdata.size:          #no atoms in phase!
1341        return {}
1342    mSize = len(Mdata)
1343    FF = np.zeros(len(Tdata))
1344    if 'NC' in calcControls[hfx+'histType']:
1345        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1346    elif 'X' in calcControls[hfx+'histType']:
1347        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1348        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1349    Uij = np.array(G2lat.U6toUij(Uijdata))
1350    bij = Mast*Uij.T
1351    dFdvDict = {}
1352    dFdfr = np.zeros((nRef,nTwin,mSize))
1353    dFdx = np.zeros((nRef,nTwin,mSize,3))
1354    dFdui = np.zeros((nRef,nTwin,mSize))
1355    dFdua = np.zeros((nRef,nTwin,mSize,6))
1356    dFdbab = np.zeros((nRef,nTwin,2))
1357    dFdtw = np.zeros((nRef,nTwin))
1358    time0 = time.time()
1359#reflection processing begins here - big arrays!
1360    iBeg = 0
1361    blkSize = 16       #no. of reflections in a block - optimized for speed
1362    while iBeg < nRef:
1363        iFin = min(iBeg+blkSize,nRef)
1364        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1365        H = refl.T[:3]
1366        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
1367        TwMask = np.any(H,axis=-1)
1368        for ir in range(blkSize):
1369            iref = ir+iBeg
1370            if iref in TwDict:
1371                for i in TwDict[iref]:
1372                    for n in range(NTL):
1373                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1374        TwMask = np.any(H,axis=-1)
1375        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1376        SQfactor = 8.0*SQ*np.pi**2
1377        if 'T' in calcControls[hfx+'histType']:
1378            if 'P' in calcControls[hfx+'histType']:
1379                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1380            else:
1381                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1382            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1383            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1384        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1385        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
1386        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1387        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1388        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
1389        Phi = np.inner(H,SGT)
1390        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1391        sinp = np.sin(phase)
1392        cosp = np.cos(phase)
1393        occ = Mdata*Fdata/len(SGT)
1394        biso = -SQfactor*Uisodata[:,nxs]
1395        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1396        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1397        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1398        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,nTwin,len(SGT),6))
1399        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1400        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
1401        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1402        fotp = FPP*Tcorr       
1403        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1404            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
1405            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
1406        else:
1407            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
1408            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
1409        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1410        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1411        if SGData['SGInv']: #centrosymmetric; B=0
1412            fbs[0] *= 0.
1413            fas[1] *= 0.
1414        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
1415        fbx = np.array([fot*cosp,-fotp*sinp])
1416        #sum below is over Uniq
1417        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
1418        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1419        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
1420        dfadx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fax[:,:,:,:,:,nxs],axis=-3),axis=0) # nRef x nTwin x nAtoms x xyz; sum on ops & A,A'
1421        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 
1422        if not SGData['SGInv']:
1423            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
1424            dfadba /= 2.
1425#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1426            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
1427            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
1428            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
1429        else:
1430            dfbdfr = np.zeros_like(dfadfr)
1431            dfbdx = np.zeros_like(dfadx)
1432            dfbdui = np.zeros_like(dfadui)
1433            dfbdua = np.zeros_like(dfadua)
1434#            dfbdba = np.zeros_like(dfadba)
1435        SA = fas[0]+fas[1]
1436        SB = fbs[0]+fbs[1]
1437        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1438        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1439        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1440        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1441        if SGData['SGInv']: #centrosymmetric; B=0
1442            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1443        else:               
1444            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1445#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1446#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1447        iBeg += blkSize
1448    print (' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0))
1449    #loop over atoms - each dict entry is list of derivatives for all the reflections
1450    for i in range(len(Mdata)):     #these all OK
1451        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1452        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1453        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1454        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1455        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1456        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1457        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1458        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1459        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1460        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1461        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1462    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1463    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1464    for i in range(nTwin):
1465        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1466    return dFdvDict
1467   
1468def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1469    '''
1470    Compute super structure factors for all h,k,l,m for phase - no twins
1471    puts the result, F^2, in each ref[9] in refList
1472    works on blocks of 32 reflections for speed
1473    input:
1474   
1475    :param dict refDict: where
1476        'RefList' list where each ref = h,k,l,m,it,d,...
1477        'FF' dict of form factors - filed in below
1478    :param np.array G:      reciprocal metric tensor
1479    :param str pfx:    phase id string
1480    :param dict SGData: space group info. dictionary output from SpcGroup
1481    :param dict calcControls:
1482    :param dict ParmDict:
1483
1484    '''
1485    phfx = pfx.split(':')[0]+hfx
1486    ast = np.sqrt(np.diag(G))
1487    GS = G/np.outer(ast,ast)
1488    uAmat,uBmat = G2lat.Gmat2AB(GS)
1489    Mast = twopisq*np.multiply.outer(ast,ast)   
1490    SGInv = SGData['SGInv']
1491    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1492    Ncen = len(SGData['SGCen'])
1493    Nops = len(SGMT)*(1+SGData['SGInv'])
1494    if SGData['SGGray']:
1495        Nops *= 2
1496    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1497    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1498    SSCen = SSGData['SSGCen']
1499    FFtables = calcControls['FFtables']
1500    BLtables = calcControls['BLtables']
1501    MFtables = calcControls['MFtables']
1502    Amat,Bmat = G2lat.Gmat2AB(G)
1503    Flack = 1.0
1504    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1505        Flack = 1.-2.*parmDict[phfx+'Flack']
1506    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1507        GetAtomFXU(pfx,calcControls,parmDict)
1508    if not Xdata.size:          #no atoms in phase!
1509        return
1510    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1511    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1512    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1513
1514    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms
1515       
1516        mXYZ = np.array([[XYZ[0] for XYZ in list(G2spc.GenAtom(xyz,SGData,All=True,Move=True))] for xyz in (Xdata+dXdata).T]) #Natn,Nop,xyz
1517        if SGData['SGGray']:
1518            mXYZ = np.hstack((mXYZ,mXYZ))
1519        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sim parts sum matches drawing
1520        Mmod = MmodA+MmodB
1521       
1522        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1523            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1524            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1525                GSdata = np.hstack((GSdata,-GSdata))     
1526            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1527            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1528            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1529            Mmod += GSdata[nxs,:,:,:]
1530           
1531#        Kdata = np.inner(Mmod,uAmat)    #Ntau,Nop,Natm,Mxyz
1532#        Mag = np.sqrt(np.sum(Kdata**2,axis=-1))
1533#        Kdata /= Mag[:,:,:,nxs]     #Cartesian unit vectors
1534#        Kdata = np.nan_to_num(Kdata)    #Ntau,Nops,Natm,Mxyz       
1535
1536    FF = np.zeros(len(Tdata))
1537    if 'NC' in calcControls[hfx+'histType']:
1538        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1539    elif 'X' in calcControls[hfx+'histType']:
1540        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1541        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1542    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1543    bij = Mast*Uij
1544    blkSize = 48       #no. of reflections in a block
1545    nRef = refDict['RefList'].shape[0]
1546    SQ = 1./(2.*refDict['RefList'].T[5])**2
1547    if 'N' in calcControls[hfx+'histType']:
1548        dat = G2el.getBLvalues(BLtables)
1549        refDict['FF']['El'] = list(dat.keys())
1550        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1551        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1552        for iel,El in enumerate(refDict['FF']['El']):
1553            if El in MFtables:
1554                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1555    else:
1556        dat = G2el.getFFvalues(FFtables,0.)
1557        refDict['FF']['El'] = list(dat.keys())
1558        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1559        for iel,El in enumerate(refDict['FF']['El']):
1560            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1561#    time0 = time.time()
1562#reflection processing begins here - big arrays!
1563    iBeg = 0
1564    while iBeg < nRef:
1565        iFin = min(iBeg+blkSize,nRef)
1566        mRef= iFin-iBeg
1567        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1568        H = refl.T[:4]                          #array(blkSize,4)
1569        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1570        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1571        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1572        Uniq = np.inner(H.T,SSGMT)
1573        UniqP = np.inner(HP.T,SGMT)
1574        Phi = np.inner(H.T,SSGT)
1575        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1576            Uniq = np.hstack((Uniq,-Uniq))
1577            Phi = np.hstack((Phi,-Phi))
1578            UniqP = np.hstack((UniqP,-UniqP))
1579        if 'T' in calcControls[hfx+'histType']:
1580            if 'P' in calcControls[hfx+'histType']:
1581                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1582            else:
1583                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1584            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1585            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1586        Bab = 0.
1587        if phfx+'BabA' in parmDict:
1588            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1589        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1590        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1591        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1592        phase = np.hstack([phase for cen in SSCen])
1593        sinp = np.sin(phase)
1594        cosp = np.cos(phase)
1595        biso = -SQfactor*Uisodata[:,nxs]
1596        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1597        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1598        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1599        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1600
1601        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: mag math here??           
1602           
1603            phasem = twopi*np.inner(HP.T,mXYZ)    #2pi(Q.r)
1604            phasem = np.swapaxes(phasem,1,2)            #Nref,Nops,Natm
1605            cosm = np.cos(phasem)
1606            sinm = np.sin(phasem)
1607            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1608            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2*Nops)     #Nref,Natm
1609                     
1610            HM = np.inner(uBmat,HP.T)                            #put into cartesian space X||H,Z||H*L
1611            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T               # normalize  HP  Nref,hkl=Unit vectors || Q
1612#for fixed moments --> m=0 reflections
1613            fam0 = 0.
1614            fbm0 = 0.
1615            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1616                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1617                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]   
1618#for modulated moments --> m != 0 reflections
1619                       
1620            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]-    \
1621                np.sign(H[3,i])*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1622                       
1623            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+    \
1624                np.sign(H[3,i])*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1625           
1626            if not SGData['SGGray']:
1627                fams += fam0[:,nxs,:,:,:]
1628                fbms += fbm0[:,nxs,:,:,:]
1629            else:
1630                fams *= Ncen
1631                fbms *= Ncen
1632               
1633# do sum on ops, atms 1st                       
1634            fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Ntau,Mxyz; sum ops & atoms
1635            fbsm = np.sum(np.sum(fbms,axis=-2),axis=-2)           
1636#put into cartesian space
1637            facm = np.inner(fasm,uAmat.T)
1638            fbcm = np.inner(fbsm,uAmat.T)
1639#form e.F dot product
1640            eDotFa = np.sum(eM[:,nxs,:]*facm,axis=-1)    #Nref,Ntau       
1641            eDotFb = np.sum(eM[:,nxs,:]*fbcm,axis=-1)
1642#intensity
1643            fass = np.sum(facm**2,axis=-1)-eDotFa**2
1644            fbss = np.sum(fbcm**2,axis=-1)-eDotFb**2
1645#do integration
1646           
1647            fas = np.sum(glWt*fass,axis=1)
1648            fbs = np.sum(glWt*fbss,axis=1)
1649           
1650            if SGData['SGInv']:
1651                fbs *= 4.
1652                fas = 0.
1653           
1654            refl.T[10] = fas+fbs   #Sum(fams**2,Mxyz) Re + Im
1655            refl.T[11] = atan2d(fbs,fas)
1656
1657        else:
1658            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1659            if 'T' in calcControls[hfx+'histType']:
1660                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1661                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1662            else:
1663                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1664                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1665            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1666            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1667            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1668            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1669           
1670            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1671            refl.T[11] = atan2d(fbs[0],fas[0])  #use only tau=0; ignore f' & f"
1672        if 'P' not in calcControls[hfx+'histType']:
1673            refl.T[8] = np.copy(refl.T[10])               
1674        iBeg += blkSize
1675#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1676    return copy.deepcopy(refDict['RefList'])
1677
1678def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1679    '''
1680    Compute super structure factors for all h,k,l,m for phase - twins only
1681    puts the result, F^2, in each ref[8+im] in refList
1682    works on blocks of 32 reflections for speed
1683    input:
1684   
1685    :param dict refDict: where
1686        'RefList' list where each ref = h,k,l,m,it,d,...
1687        'FF' dict of form factors - filed in below
1688    :param np.array G:      reciprocal metric tensor
1689    :param str pfx:    phase id string
1690    :param dict SGData: space group info. dictionary output from SpcGroup
1691    :param dict calcControls:
1692    :param dict ParmDict:
1693
1694    '''
1695    phfx = pfx.split(':')[0]+hfx
1696    ast = np.sqrt(np.diag(G))
1697    Mast = twopisq*np.multiply.outer(ast,ast)   
1698    SGInv = SGData['SGInv']
1699    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1700    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1701    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1702    FFtables = calcControls['FFtables']
1703    BLtables = calcControls['BLtables']
1704    MFtables = calcControls['MFtables']
1705    Flack = 1.0
1706    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1707        Flack = 1.-2.*parmDict[phfx+'Flack']
1708    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1709    TwDict = refDict.get('TwDict',{})           
1710    if 'S' in calcControls[hfx+'histType']:
1711        NTL = calcControls[phfx+'NTL']
1712        NM = calcControls[phfx+'TwinNMN']+1
1713        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1714        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1715        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1716    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1717        GetAtomFXU(pfx,calcControls,parmDict)
1718    if not Xdata.size:          #no atoms in phase!
1719        return
1720    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1721    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1722    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1723    FF = np.zeros(len(Tdata))
1724    if 'NC' in calcControls[hfx+'histType']:
1725        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1726    elif 'X' in calcControls[hfx+'histType']:
1727        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1728        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1729    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1730    bij = Mast*Uij
1731    blkSize = 32       #no. of reflections in a block
1732    nRef = refDict['RefList'].shape[0]
1733    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1734        SQ = 1./(2.*refDict['RefList'].T[5])**2
1735        if 'N' in calcControls[hfx+'histType']:
1736            dat = G2el.getBLvalues(BLtables)
1737            refDict['FF']['El'] = list(dat.keys())
1738            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1739            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1740            for iel,El in enumerate(refDict['FF']['El']):
1741                if El in MFtables:
1742                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1743        else:
1744            dat = G2el.getFFvalues(FFtables,0.)
1745            refDict['FF']['El'] = list(dat.keys())
1746            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1747            for iel,El in enumerate(refDict['FF']['El']):
1748                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1749#    time0 = time.time()
1750#reflection processing begins here - big arrays!
1751    iBeg = 0
1752    while iBeg < nRef:
1753        iFin = min(iBeg+blkSize,nRef)
1754        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1755        H = refl[:,:4]                          #array(blkSize,4)
1756        H3 = refl[:,:3]
1757        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1758        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1759        H3 = np.inner(H3,TwinLaw)       
1760        TwMask = np.any(HP,axis=-1)
1761        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1762            for ir in range(blkSize):
1763                iref = ir+iBeg
1764                if iref in TwDict:
1765                    for i in TwDict[iref]:
1766                        for n in range(NTL):
1767                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1768                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1769            TwMask = np.any(HP,axis=-1)
1770        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1771        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1772        Uniq = np.inner(H,SSGMT)
1773        Uniq3 = np.inner(H3,SGMT)
1774        UniqP = np.inner(HP,SGMT)
1775        Phi = np.inner(H,SSGT)
1776        if SGInv:   #if centro - expand HKL sets
1777            Uniq = np.hstack((Uniq,-Uniq))
1778            Uniq3 = np.hstack((Uniq3,-Uniq3))
1779            Phi = np.hstack((Phi,-Phi))
1780            UniqP = np.hstack((UniqP,-UniqP))
1781        if 'T' in calcControls[hfx+'histType']:
1782            if 'P' in calcControls[hfx+'histType']:
1783                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1784            else:
1785                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1786            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1787            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1788        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1789        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1790        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1791        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1792        sinp = np.sin(phase)
1793        cosp = np.cos(phase)
1794        biso = -SQfactor*Uisodata[:,nxs]
1795        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1796        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1797        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1798        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1799        if 'T' in calcControls[hfx+'histType']:
1800            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1801            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1802        else:
1803            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1804            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1805        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1806        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1807        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1808        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1809        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1810        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1811        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1812            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1813        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1814        iBeg += blkSize
1815#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1816
1817def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1818    '''
1819    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1820    Only Fourier component are done analytically here
1821    input:
1822   
1823    :param dict refDict: where
1824        'RefList' list where each ref = h,k,l,m,it,d,...
1825        'FF' dict of form factors - filled in below
1826    :param int im: = 1 (could be eliminated)
1827    :param np.array G:      reciprocal metric tensor
1828    :param str hfx:    histogram id string
1829    :param str pfx:    phase id string
1830    :param dict SGData: space group info. dictionary output from SpcGroup
1831    :param dict SSGData: super space group info.
1832    :param dict calcControls:
1833    :param dict ParmDict:
1834   
1835    :returns: dict dFdvDict: dictionary of derivatives
1836    '''
1837    phfx = pfx.split(':')[0]+hfx
1838    ast = np.sqrt(np.diag(G))
1839    Mast = twopisq*np.multiply.outer(ast,ast)
1840    SGInv = SGData['SGInv']
1841    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1842    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1843    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1844    FFtables = calcControls['FFtables']
1845    BLtables = calcControls['BLtables']
1846    nRef = len(refDict['RefList'])
1847    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1848        GetAtomFXU(pfx,calcControls,parmDict)
1849    if not Xdata.size:          #no atoms in phase!
1850        return {}
1851    mSize = len(Mdata)  #no. atoms
1852    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1853    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1854    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1855    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1856    FF = np.zeros(len(Tdata))
1857    if 'NC' in calcControls[hfx+'histType']:
1858        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1859    elif 'X' in calcControls[hfx+'histType']:
1860        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1861        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1862    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1863    bij = Mast*Uij
1864    if not len(refDict['FF']):
1865        if 'N' in calcControls[hfx+'histType']:
1866            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1867        else:
1868            dat = G2el.getFFvalues(FFtables,0.)       
1869        refDict['FF']['El'] = list(dat.keys())
1870        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1871    dFdvDict = {}
1872    dFdfr = np.zeros((nRef,mSize))
1873    dFdx = np.zeros((nRef,mSize,3))
1874    dFdui = np.zeros((nRef,mSize))
1875    dFdua = np.zeros((nRef,mSize,6))
1876    dFdbab = np.zeros((nRef,2))
1877    dFdfl = np.zeros((nRef))
1878    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1879    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1880    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1881    Flack = 1.0
1882    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1883        Flack = 1.-2.*parmDict[phfx+'Flack']
1884    time0 = time.time()
1885    nRef = len(refDict['RefList'])/100
1886    for iref,refl in enumerate(refDict['RefList']):
1887        if 'T' in calcControls[hfx+'histType']:
1888            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1889        H = np.array(refl[:4])
1890        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1891        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1892        SQfactor = 8.0*SQ*np.pi**2
1893        Bab = 0.0
1894        if phfx+'BabA' in parmDict:
1895            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1896            Bab = parmDict[phfx+'BabA']*dBabdA
1897        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1898        FF = refDict['FF']['FF'][iref].T[Tindx]
1899        Uniq = np.inner(H,SSGMT)
1900        Phi = np.inner(H,SSGT)
1901        UniqP = np.inner(HP,SGMT)
1902        if SGInv:   #if centro - expand HKL sets
1903            Uniq = np.vstack((Uniq,-Uniq))
1904            Phi = np.hstack((Phi,-Phi))
1905            UniqP = np.vstack((UniqP,-UniqP))
1906        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1907        sinp = np.sin(phase)
1908        cosp = np.cos(phase)
1909        occ = Mdata*Fdata/Uniq.shape[0]
1910        biso = -SQfactor*Uisodata[:,nxs]
1911        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1912        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1913        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1914        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1915        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1916        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1917        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1918        fotp = FPP*Tcorr            #ops x atoms
1919        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1920        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1921        # GfpuA is 2 x ops x atoms
1922        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1923        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1924        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1925        fag = fa*GfpuA[0]-fb*GfpuA[1]
1926        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1927       
1928        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1929        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1930        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1931        fbx = np.array([fot*cosp,-fotp*sinp])
1932        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1933        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1934        #sum below is over Uniq
1935        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1936        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1937        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1938        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1939        dfadui = np.sum(-SQfactor*fag,axis=1)
1940        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1941        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1942        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1943        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1944        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1945        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1946        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1947        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1948        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1949        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1950        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1951        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1952        if not SGData['SGInv']:   #Flack derivative
1953            dfadfl = np.sum(-FPP*Tcorr*sinp)
1954            dfbdfl = np.sum(FPP*Tcorr*cosp)
1955        else:
1956            dfadfl = 1.0
1957            dfbdfl = 1.0
1958        SA = fas[0]+fas[1]      #float = A+A'
1959        SB = fbs[0]+fbs[1]      #float = B+B'
1960        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1961            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1962            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1963                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1964            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1965                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1966            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1967                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1968            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1969                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1970            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1971                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1972            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1973                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1974            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1975                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1976        else:                       #OK, I think
1977            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1978            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1979            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1980            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1981            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1982                           
1983            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1984            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1985            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1986        if phfx+'BabA' in parmDict:
1987            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1988                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1989        #loop over atoms - each dict entry is list of derivatives for all the reflections
1990        if not iref%100 :
1991            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1992    for i in range(len(Mdata)):     #loop over atoms
1993        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1994        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1995        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1996        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1997        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1998        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1999        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2000        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2001        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2002        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2003        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2004        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2005            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2006            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2007        nx = 0
2008        if waveTypes[i] in ['Block','ZigZag']:
2009            nx = 1 
2010        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2011            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2012            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2013            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2014            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2015            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2016            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2017        for j in range(USSdata.shape[1]):       #loop over waves
2018            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2019            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2020            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2021            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2022            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2023            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2024            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2025            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2026            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2027            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2028            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2029            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2030           
2031    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
2032    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2033    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2034    return dFdvDict
2035
2036def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2037    '''
2038    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
2039    input:
2040   
2041    :param dict refDict: where
2042        'RefList' list where each ref = h,k,l,m,it,d,...
2043        'FF' dict of form factors - filled in below
2044    :param int im: = 1 (could be eliminated)
2045    :param np.array G:      reciprocal metric tensor
2046    :param str hfx:    histogram id string
2047    :param str pfx:    phase id string
2048    :param dict SGData: space group info. dictionary output from SpcGroup
2049    :param dict SSGData: super space group info.
2050    :param dict calcControls:
2051    :param dict ParmDict:
2052   
2053    :returns: dict dFdvDict: dictionary of derivatives
2054    '''
2055
2056    trefDict = copy.deepcopy(refDict)
2057    dM = 1.e-4
2058    dFdvDict = {}
2059    for parm in parmDict:
2060        if ':' not in parm:
2061            continue
2062        if parm.split(':')[2] in ['Tmin','Tmax','Xmax','Ymax','Zmax','Fzero','Fwid',
2063            'MXsin','MXcos','MYsin','MYcos','MZsin','MZcos','AMx','AMy','AMz',]:
2064            parmDict[parm] += dM
2065            prefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2066            parmDict[parm] -= 2*dM
2067            mrefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2068            parmDict[parm] += dM
2069            dFdvDict[parm] = (prefList[:,9+im]-mrefList[:,9+im])/(2.*dM)
2070    return dFdvDict
2071   
2072def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2073    'Needs a doc string'
2074    phfx = pfx.split(':')[0]+hfx
2075    ast = np.sqrt(np.diag(G))
2076    Mast = twopisq*np.multiply.outer(ast,ast)
2077    SGInv = SGData['SGInv']
2078    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2079    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2080    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2081    FFtables = calcControls['FFtables']
2082    BLtables = calcControls['BLtables']
2083    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2084    TwDict = refDict.get('TwDict',{})           
2085    if 'S' in calcControls[hfx+'histType']:
2086        NTL = calcControls[phfx+'NTL']
2087        NM = calcControls[phfx+'TwinNMN']+1
2088        TwinLaw = calcControls[phfx+'TwinLaw']
2089        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2090    nTwin = len(TwinLaw)       
2091    nRef = len(refDict['RefList'])
2092    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2093        GetAtomFXU(pfx,calcControls,parmDict)
2094    if not Xdata.size:          #no atoms in phase!
2095        return {}
2096    mSize = len(Mdata)  #no. atoms
2097    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2098    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
2099    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2100    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2101    FF = np.zeros(len(Tdata))
2102    if 'NC' in calcControls[hfx+'histType']:
2103        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2104    elif 'X' in calcControls[hfx+'histType']:
2105        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2106        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2107    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2108    bij = Mast*Uij
2109    if not len(refDict['FF']):
2110        if 'N' in calcControls[hfx+'histType']:
2111            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2112        else:
2113            dat = G2el.getFFvalues(FFtables,0.)       
2114        refDict['FF']['El'] = list(dat.keys())
2115        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2116    dFdvDict = {}
2117    dFdfr = np.zeros((nRef,nTwin,mSize))
2118    dFdx = np.zeros((nRef,nTwin,mSize,3))
2119    dFdui = np.zeros((nRef,nTwin,mSize))
2120    dFdua = np.zeros((nRef,nTwin,mSize,6))
2121    dFdbab = np.zeros((nRef,nTwin,2))
2122    dFdtw = np.zeros((nRef,nTwin))
2123    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2124    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2125    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2126    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2127    Flack = 1.0
2128    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2129        Flack = 1.-2.*parmDict[phfx+'Flack']
2130    time0 = time.time()
2131    nRef = len(refDict['RefList'])/100
2132    for iref,refl in enumerate(refDict['RefList']):
2133        if 'T' in calcControls[hfx+'histType']:
2134            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2135        H = np.array(refl[:4])
2136        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2137        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2138        TwMask = np.any(H,axis=-1)
2139        if TwinLaw.shape[0] > 1 and TwDict:
2140            if iref in TwDict:
2141                for i in TwDict[iref]:
2142                    for n in range(NTL):
2143                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2144            TwMask = np.any(H,axis=-1)
2145        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2146        SQfactor = 8.0*SQ*np.pi**2
2147        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2148        Bab = parmDict[phfx+'BabA']*dBabdA
2149        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2150        FF = refDict['FF']['FF'][iref].T[Tindx]
2151        Uniq = np.inner(H,SSGMT)
2152        Phi = np.inner(H,SSGT)
2153        UniqP = np.inner(HP,SGMT)
2154        if SGInv:   #if centro - expand HKL sets
2155            Uniq = np.vstack((Uniq,-Uniq))
2156            Phi = np.hstack((Phi,-Phi))
2157            UniqP = np.vstack((UniqP,-UniqP))
2158        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2159        sinp = np.sin(phase)
2160        cosp = np.cos(phase)
2161        occ = Mdata*Fdata/Uniq.shape[0]
2162        biso = -SQfactor*Uisodata[:,nxs]
2163        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2164        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2165        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2166        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2167        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2168        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2169        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2170        fotp = FPP*Tcorr            #ops x atoms
2171        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2172        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2173        # GfpuA is 2 x ops x atoms
2174        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2175        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2176        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2177        fag = fa*GfpuA[0]-fb*GfpuA[1]
2178        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2179       
2180        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2181        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2182        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2183        fbx = np.array([fot*cosp,-fotp*sinp])
2184        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2185        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2186        #sum below is over Uniq
2187        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2188        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2189        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2190        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2191        dfadui = np.sum(-SQfactor*fag,axis=1)
2192        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2193        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2194        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2195        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2196        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2197        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2198        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2199        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2200        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2201        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2202        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2203        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2204        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2205        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2206#        GSASIIpath.IPyBreak()
2207        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2208        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2209        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2210        dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0][it]+SA[it]*dfadfr[1][it]+SB[it]*dfbdfr[0][it]+SB[it]*dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
2211        dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it][0]+SA[it]*dfadx[it][1]+SB[it]*dfbdx[it][0]+SB[it]*dfbdx[it][1]) for it in range(nTwin)]
2212        dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[it][0]+SA[it]*dfadui[it][1]+SB[it]*dfbdui[it][0]+SB[it]*dfbdui[it][1]) for it in range(nTwin)]
2213        dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it][0]+SA[it]*dfadua[it][1]+SB[it]*dfbdua[it][0]+SB[it]*dfbdua[it][1]) for it in range(nTwin)]
2214        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2215
2216        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2217        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2218        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2219        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2220#            GSASIIpath.IPyBreak()
2221        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2222            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2223        #loop over atoms - each dict entry is list of derivatives for all the reflections
2224        if not iref%100 :
2225            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2226    for i in range(len(Mdata)):     #loop over atoms
2227        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2228        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2229        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2230        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2231        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2232        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2233        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2234        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2235        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2236        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2237        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2238        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2239            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2240            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2241        nx = 0
2242        if waveTypes[i] in ['Block','ZigZag']:
2243            nx = 1 
2244            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2245            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2246            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2247            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2248            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2249        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2250            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2251            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2252            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2253            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2254            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2255            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2256        for j in range(USSdata.shape[1]):       #loop over waves
2257            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2258            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2259            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2260            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2261            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2262            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2263            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2264            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2265            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2266            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2267            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2268            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2269           
2270#        GSASIIpath.IPyBreak()
2271    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2272    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2273    return dFdvDict
2274   
2275def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2276    ''' Single crystal extinction function; returns extinction & derivative
2277    '''
2278    extCor = 1.0
2279    dervDict = {}
2280    dervCor = 1.0
2281    if calcControls[phfx+'EType'] != 'None':
2282        SQ = 1/(4.*ref[4+im]**2)
2283        if 'C' in parmDict[hfx+'Type']:           
2284            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2285        else:   #'T'
2286            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2287        if 'SXC' in parmDict[hfx+'Type']:
2288            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2289            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2290            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2291            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2292        elif 'SNT' in parmDict[hfx+'Type']:
2293            AV = 1.e7/parmDict[pfx+'Vol']**2
2294            PL = SQ
2295            PLZ = AV*ref[9+im]*ref[12+im]**2
2296        elif 'SNC' in parmDict[hfx+'Type']:
2297            AV = 1.e7/parmDict[pfx+'Vol']**2
2298            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2299            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2300           
2301        if 'Primary' in calcControls[phfx+'EType']:
2302            PLZ *= 1.5
2303        else:
2304            if 'C' in parmDict[hfx+'Type']:
2305                PLZ *= calcControls[phfx+'Tbar']
2306            else: #'T'
2307                PLZ *= ref[13+im]      #t-bar
2308        if 'Primary' in calcControls[phfx+'EType']:
2309            PLZ *= 1.5
2310            PSIG = parmDict[phfx+'Ep']
2311        elif 'I & II' in calcControls[phfx+'EType']:
2312            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2313        elif 'Type II' in calcControls[phfx+'EType']:
2314            PSIG = parmDict[phfx+'Es']
2315        else:       # 'Secondary Type I'
2316            PSIG = parmDict[phfx+'Eg']/PL
2317           
2318        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2319        AL = 0.025+0.285*cos2T
2320        BG = 0.02-0.025*cos2T
2321        BL = 0.15-0.2*(0.75-cos2T)**2
2322        if cos2T < 0.:
2323            BL = -0.45*cos2T
2324        CG = 2.
2325        CL = 2.
2326        PF = PLZ*PSIG
2327       
2328        if 'Gaussian' in calcControls[phfx+'EApprox']:
2329            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2330            extCor = np.sqrt(PF4)
2331            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2332        else:
2333            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2334            extCor = np.sqrt(PF4)
2335            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2336
2337        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2338        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2339            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2340        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2341            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2342        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2343            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2344               
2345    return 1./extCor,dervDict,dervCor
2346   
2347def Dict2Values(parmdict, varylist):
2348    '''Use before call to leastsq to setup list of values for the parameters
2349    in parmdict, as selected by key in varylist'''
2350    return [parmdict[key] for key in varylist] 
2351   
2352def Values2Dict(parmdict, varylist, values):
2353    ''' Use after call to leastsq to update the parameter dictionary with
2354    values corresponding to keys in varylist'''
2355    parmdict.update(zip(varylist,values))
2356   
2357def GetNewCellParms(parmDict,varyList):
2358    '''Compute unit cell tensor terms from varied Aij and Dij values.
2359    Terms are included in the dict only if Aij or Dij is varied.
2360    '''
2361    newCellDict = {}
2362    Anames = ['A'+str(i) for i in range(6)]
2363    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2364    for item in varyList:
2365        keys = item.split(':')
2366        if keys[2] in Ddict:
2367            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2368            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2369            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2370    return newCellDict          # is e.g. {'0::D11':A0-D11}
2371   
2372def ApplyXYZshifts(parmDict,varyList):
2373    '''
2374    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2375   
2376    :param dict parmDict: parameter dictionary
2377    :param list varyList: list of variables (not used!)
2378    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2379
2380    '''
2381    newAtomDict = {}
2382    for item in parmDict:
2383        if 'dA' in item:
2384            parm = ''.join(item.split('d'))
2385            parmDict[parm] += parmDict[item]
2386            newAtomDict[item] = [parm,parmDict[parm]]
2387    return newAtomDict
2388   
2389def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2390    'Spherical harmonics texture'
2391    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2392    if 'T' in calcControls[hfx+'histType']:
2393        tth = parmDict[hfx+'2-theta']
2394    else:
2395        tth = refl[5+im]
2396    odfCor = 1.0
2397    H = refl[:3]
2398    cell = G2lat.Gmat2cell(g)
2399    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2400    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2401    phi,beta = G2lat.CrsAng(H,cell,SGData)
2402    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2403    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2404    for item in SHnames:
2405        L,M,N = eval(item.strip('C'))
2406        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2407        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2408        Lnorm = G2lat.Lnorm(L)
2409        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2410    return odfCor
2411   
2412def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2413    'Spherical harmonics texture derivatives'
2414    if 'T' in calcControls[hfx+'histType']:
2415        tth = parmDict[hfx+'2-theta']
2416    else:
2417        tth = refl[5+im]
2418    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2419    odfCor = 1.0
2420    dFdODF = {}
2421    dFdSA = [0,0,0]
2422    H = refl[:3]
2423    cell = G2lat.Gmat2cell(g)
2424    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2425    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2426    phi,beta = G2lat.CrsAng(H,cell,SGData)
2427    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2428    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2429    for item in SHnames:
2430        L,M,N = eval(item.strip('C'))
2431        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2432        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2433        Lnorm = G2lat.Lnorm(L)
2434        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2435        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2436        for i in range(3):
2437            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2438    return odfCor,dFdODF,dFdSA
2439   
2440def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2441    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2442    if 'T' in calcControls[hfx+'histType']:
2443        tth = parmDict[hfx+'2-theta']
2444    else:
2445        tth = refl[5+im]
2446    odfCor = 1.0
2447    H = refl[:3]
2448    cell = G2lat.Gmat2cell(g)
2449    Sangls = [0.,0.,0.]
2450    if 'Bragg' in calcControls[hfx+'instType']:
2451        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2452        IFCoup = True
2453    else:
2454        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2455        IFCoup = False
2456    phi,beta = G2lat.CrsAng(H,cell,SGData)
2457    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2458    SHnames = calcControls[phfx+'SHnames']
2459    for item in SHnames:
2460        L,N = eval(item.strip('C'))
2461        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2462        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2463        Lnorm = G2lat.Lnorm(L)
2464        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2465    return np.squeeze(odfCor)
2466   
2467def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2468    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2469    if 'T' in calcControls[hfx+'histType']:
2470        tth = parmDict[hfx+'2-theta']
2471    else:
2472        tth = refl[5+im]
2473    odfCor = 1.0
2474    dFdODF = {}
2475    H = refl[:3]
2476    cell = G2lat.Gmat2cell(g)
2477    Sangls = [0.,0.,0.]
2478    if 'Bragg' in calcControls[hfx+'instType']:
2479        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2480        IFCoup = True
2481    else:
2482        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2483        IFCoup = False
2484    phi,beta = G2lat.CrsAng(H,cell,SGData)
2485    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2486    SHnames = calcControls[phfx+'SHnames']
2487    for item in SHnames:
2488        L,N = eval(item.strip('C'))
2489        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2490        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2491        Lnorm = G2lat.Lnorm(L)
2492        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2493        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2494    return odfCor,dFdODF
2495   
2496def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2497    'March-Dollase preferred orientation correction'
2498    POcorr = 1.0
2499    MD = parmDict[phfx+'MD']
2500    if MD != 1.0:
2501        MDAxis = calcControls[phfx+'MDAxis']
2502        sumMD = 0
2503        for H in uniq:           
2504            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2505            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2506            sumMD += A**3
2507        POcorr = sumMD/len(uniq)
2508    return POcorr
2509   
2510def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2511    'Needs a doc string'
2512    POcorr = 1.0
2513    POderv = {}
2514    if calcControls[phfx+'poType'] == 'MD':
2515        MD = parmDict[phfx+'MD']
2516        MDAxis = calcControls[phfx+'MDAxis']
2517        sumMD = 0
2518        sumdMD = 0
2519        for H in uniq:           
2520            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2521            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2522            sumMD += A**3
2523            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2524        POcorr = sumMD/len(uniq)
2525        POderv[phfx+'MD'] = sumdMD/len(uniq)
2526    else:   #spherical harmonics
2527        if calcControls[phfx+'SHord']:
2528            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2529    return POcorr,POderv
2530   
2531def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2532    'Needs a doc string'
2533    if 'Debye' in calcControls[hfx+'instType']:
2534        if 'T' in calcControls[hfx+'histType']:
2535            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2536        else:
2537            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2538    else:
2539        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2540   
2541def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2542    'Needs a doc string'
2543    if 'Debye' in calcControls[hfx+'instType']:
2544        if 'T' in calcControls[hfx+'histType']:
2545            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2546        else:
2547            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2548    else:
2549        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2550       
2551def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2552    'Needs a doc string'
2553    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2554    pi2 = np.sqrt(2./np.pi)
2555    if 'T' in calcControls[hfx+'histType']:
2556        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2557        wave = refl[14+im]
2558    else:   #'C'W
2559        sth2 = sind(refl[5+im]/2.)**2
2560        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2561    c2th = 1.-2.0*sth2
2562    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2563    if 'X' in calcControls[hfx+'histType']:
2564        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2565    xfac = flv2*parmDict[phfx+'Extinction']
2566    exb = 1.0
2567    if xfac > -1.:
2568        exb = 1./np.sqrt(1.+xfac)
2569    exl = 1.0
2570    if 0 < xfac <= 1.:
2571        xn = np.array([xfac**(i+1) for i in range(6)])
2572        exl += np.sum(xn*coef)
2573    elif xfac > 1.:
2574        xfac2 = 1./np.sqrt(xfac)
2575        exl = pi2*(1.-0.125/xfac)*xfac2
2576    return exb*sth2+exl*(1.-sth2)
2577   
2578def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2579    'Needs a doc string'
2580    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2581    pi2 = np.sqrt(2./np.pi)
2582    if 'T' in calcControls[hfx+'histType']:
2583        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2584        wave = refl[14+im]
2585    else:   #'C'W
2586        sth2 = sind(refl[5+im]/2.)**2
2587        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2588    c2th = 1.-2.0*sth2
2589    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2590    if 'X' in calcControls[hfx+'histType']:
2591        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2592    xfac = flv2*parmDict[phfx+'Extinction']
2593    dbde = -500.*flv2
2594    if xfac > -1.:
2595        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2596    dlde = 0.
2597    if 0 < xfac <= 1.:
2598        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2599        dlde = np.sum(xn*coef)/xfac
2600    elif xfac > 1.:
2601        xfac2 = 1./np.sqrt(xfac)
2602        dlde = 0.5*flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2603       
2604    return dbde*sth2+dlde*(1.-sth2)
2605   
2606def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2607    'Needs a doc string'    #need powder extinction!
2608    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2609    if 'X' in parmDict[hfx+'Type']:
2610        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2611    POcorr = 1.0
2612    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2613        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2614    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2615        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2616    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2617        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2618    Icorr *= POcorr
2619    AbsCorr = 1.0
2620    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2621    Icorr *= AbsCorr
2622    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2623    Icorr *= ExtCorr
2624    return Icorr,POcorr,AbsCorr,ExtCorr
2625   
2626def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2627    'Needs a doc string'    #need powder extinction derivs!
2628    dIdsh = 1./parmDict[hfx+'Scale']
2629    dIdsp = 1./parmDict[phfx+'Scale']
2630    if 'X' in parmDict[hfx+'Type']:
2631        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2632        dIdPola /= pola
2633    else:       #'N'
2634        dIdPola = 0.0
2635    dFdODF = {}
2636    dFdSA = [0,0,0]
2637    dIdPO = {}
2638    if pfx+'SHorder' in parmDict:
2639        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2640        for iSH in dFdODF:
2641            dFdODF[iSH] /= odfCor
2642        for i in range(3):
2643            dFdSA[i] /= odfCor
2644    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2645        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2646        for iPO in dIdPO:
2647            dIdPO[iPO] /= POcorr
2648    if 'T' in parmDict[hfx+'Type']:
2649        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2650        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2651    else:
2652        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2653        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2654    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2655       
2656def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2657    'Needs a doc string'
2658    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2659        costh = cosd(refl[5+im]/2.)
2660        #crystallite size
2661        if calcControls[phfx+'SizeType'] == 'isotropic':
2662            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2663        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2664            H = np.array(refl[:3])
2665            P = np.array(calcControls[phfx+'SizeAxis'])
2666            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2667            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2668            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2669        else:           #ellipsoidal crystallites
2670            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2671            H = np.array(refl[:3])
2672            lenR = G2pwd.ellipseSize(H,Sij,GB)
2673            Sgam = 1.8*wave/(np.pi*costh*lenR)
2674        #microstrain               
2675        if calcControls[phfx+'MustrainType'] == 'isotropic':
2676            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2677        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2678            H = np.array(refl[:3])
2679            P = np.array(calcControls[phfx+'MustrainAxis'])
2680            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2681            Si = parmDict[phfx+'Mustrain;i']
2682            Sa = parmDict[phfx+'Mustrain;a']
2683            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2684        else:       #generalized - P.W. Stephens model
2685            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2686            Sum = 0
2687            for i,strm in enumerate(Strms):
2688                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2689            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2690    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2691        #crystallite size
2692        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2693            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2694        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2695            H = np.array(refl[:3])
2696            P = np.array(calcControls[phfx+'SizeAxis'])
2697            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2698            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2699            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2700        else:           #ellipsoidal crystallites   #OK
2701            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2702            H = np.array(refl[:3])
2703            lenR = G2pwd.ellipseSize(H,Sij,GB)
2704            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2705        #microstrain               
2706        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2707            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2708        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2709            H = np.array(refl[:3])
2710            P = np.array(calcControls[phfx+'MustrainAxis'])
2711            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2712            Si = parmDict[phfx+'Mustrain;i']
2713            Sa = parmDict[phfx+'Mustrain;a']
2714            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2715        else:       #generalized - P.W. Stephens model  OK
2716            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2717            Sum = 0
2718            for i,strm in enumerate(Strms):
2719                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2720            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2721           
2722    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2723    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2724    sig /= ateln2
2725    return sig,gam
2726       
2727def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2728    'Needs a doc string'
2729    gamDict = {}
2730    sigDict = {}
2731    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2732        costh = cosd(refl[5+im]/2.)
2733        tanth = tand(refl[5+im]/2.)
2734        #crystallite size derivatives
2735        if calcControls[phfx+'SizeType'] == 'isotropic':
2736            Sgam = 1.8*wave/(np.pi*costh*parmDict[phfx+'Size;i'])
2737            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh*parmDict[phfx+'Size;i']**2)
2738            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2739        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2740            H = np.array(refl[:3])
2741            P = np.array(calcControls[phfx+'SizeAxis'])
2742            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2743            Si = parmDict[phfx+'Size;i']
2744            Sa = parmDict[phfx+'Size;a']
2745            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2746            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2747            Sgam = gami*sqtrm
2748            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2749            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2750            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2751            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2752            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2753            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2754        else:           #ellipsoidal crystallites
2755            const = 1.8*wave/(np.pi*costh)
2756            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2757            H = np.array(refl[:3])
2758            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2759            Sgam = const/lenR
2760            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2761                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2762                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2763        gamDict[phfx+'Size;mx'] = Sgam
2764        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2765               
2766        #microstrain derivatives               
2767        if calcControls[phfx+'MustrainType'] == 'isotropic':
2768            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2769            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2770            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2771        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2772            H = np.array(refl[:3])
2773            P = np.array(calcControls[phfx+'MustrainAxis'])
2774            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2775            Si = parmDict[phfx+'Mustrain;i']
2776            Sa = parmDict[phfx+'Mustrain;a']
2777            gami = 0.018*Si*Sa*tanth/np.pi
2778            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2779            Mgam = gami/sqtrm
2780            dsi = -gami*Si*cosP**2/sqtrm**3
2781            dsa = -gami*Sa*sinP**2/sqtrm**3
2782            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2783            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2784            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2785            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2786        else:       #generalized - P.W. Stephens model
2787            const = 0.018*refl[4+im]**2*tanth/np.pi
2788            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2789            Sum = 0
2790            for i,strm in enumerate(Strms):
2791                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2792                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2793                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2794            Mgam = const*np.sqrt(Sum)
2795            for i in range(len(Strms)):
2796                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2797                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2798        gamDict[phfx+'Mustrain;mx'] = Mgam
2799        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2800    else:   #'T'OF - All checked & OK
2801        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2802            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2803            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2804            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2805        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2806            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2807            H = np.array(refl[:3])
2808            P = np.array(calcControls[phfx+'SizeAxis'])
2809            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2810            Si = parmDict[phfx+'Size;i']
2811            Sa = parmDict[phfx+'Size;a']
2812            gami = const/(Si*Sa)
2813            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2814            Sgam = gami*sqtrm
2815            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2816            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2817            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2818            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2819            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2820            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2821        else:           #OK  ellipsoidal crystallites
2822            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2823            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2824            H = np.array(refl[:3])
2825            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2826            Sgam = const/lenR
2827            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2828                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2829                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2830        gamDict[phfx+'Size;mx'] = Sgam  #OK
2831        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2832               
2833        #microstrain derivatives               
2834        if calcControls[phfx+'MustrainType'] == 'isotropic':
2835            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2836            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2837            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2838        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2839            H = np.array(refl[:3])
2840            P = np.array(calcControls[phfx+'MustrainAxis'])
2841            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2842            Si = parmDict[phfx+'Mustrain;i']
2843            Sa = parmDict[phfx+'Mustrain;a']
2844            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2845            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2846            Mgam = gami/sqtrm
2847            dsi = -gami*Si*cosP**2/sqtrm**3
2848            dsa = -gami*Sa*sinP**2/sqtrm**3
2849            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2850            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2851            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2852            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2853        else:       #generalized - P.W. Stephens model OK
2854            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2855            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2856            Sum = 0
2857            for i,strm in enumerate(Strms):
2858                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2859                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2860                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2861            Mgam = const*np.sqrt(Sum)
2862            for i in range(len(Strms)):
2863                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2864                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2865        gamDict[phfx+'Mustrain;mx'] = Mgam
2866        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2867       
2868    return sigDict,gamDict
2869       
2870def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2871    'Needs a doc string'
2872    if im:
2873        h,k,l,m = refl[:4]
2874        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2875        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2876    else:
2877        h,k,l = refl[:3]
2878        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2879    refl[4+im] = d
2880    if 'C' in calcControls[hfx+'histType']:
2881        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2882        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2883        if 'Bragg' in calcControls[hfx+'instType']:
2884            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2885                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2886        else:               #Debye-Scherrer - simple but maybe not right
2887            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2888    elif 'T' in calcControls[hfx+'histType']:
2889        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2890        #do I need sample position effects - maybe?
2891    return pos
2892
2893def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2894    'Needs a doc string'
2895    dpr = 180./np.pi
2896    if im:
2897        h,k,l,m = refl[:4]
2898        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2899        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2900        h,k,l = [h+m*vec[0],k+m*vec[1],l+m*vec[2]]          #do proj of hklm to hkl so dPdA & dPdV come out right
2901    else:
2902        m = 0
2903        h,k,l = refl[:3]       
2904        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2905    dst = np.sqrt(dstsq)
2906    dsp = 1./dst
2907    if 'C' in calcControls[hfx+'histType']:
2908        pos = refl[5+im]-parmDict[hfx+'Zero']
2909        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2910        dpdw = const*dst
2911        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2912        dpdZ = 1.0
2913        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2914            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2915        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2916        if 'Bragg' in calcControls[hfx+'instType']:
2917            dpdSh = -4.*shft*cosd(pos/2.0)
2918            dpdTr = -shft*sind(pos)*100.0
2919            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2920        else:               #Debye-Scherrer - simple but maybe not right
2921            dpdXd = -shft*cosd(pos)
2922            dpdYd = -shft*sind(pos)
2923            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2924    elif 'T' in calcControls[hfx+'histType']:
2925        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2926        dpdZ = 1.0
2927        dpdDC = dsp
2928        dpdDA = dsp**2
2929        dpdDB = 1./dsp
2930        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2931            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2932        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2933           
2934def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2935    'Needs a doc string'
2936    laue = SGData['SGLaue']
2937    uniq = SGData['SGUniq']
2938    h,k,l = refl[:3]
2939    if laue in ['m3','m3m']:
2940        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2941            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2942    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2943        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2944    elif laue in ['3R','3mR']:
2945        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2946    elif laue in ['4/m','4/mmm']:
2947        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2948    elif laue in ['mmm']:
2949        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2950    elif laue in ['2/m']:
2951        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2952        if uniq == 'a':
2953            Dij += parmDict[phfx+'D23']*k*l
2954        elif uniq == 'b':
2955            Dij += parmDict[phfx+'D13']*h*l
2956        elif uniq == 'c':
2957            Dij += parmDict[phfx+'D12']*h*k
2958    else:
2959        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2960            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2961    if 'C' in calcControls[hfx+'histType']:
2962        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2963    else:
2964        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2965           
2966def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2967    'Needs a doc string'
2968    laue = SGData['SGLaue']
2969    uniq = SGData['SGUniq']
2970    h,k,l = refl[:3]
2971    if laue in ['m3','m3m']:
2972        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2973            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2974    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2975        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2976    elif laue in ['3R','3mR']:
2977        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2978    elif laue in ['4/m','4/mmm']:
2979        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2980    elif laue in ['mmm']:
2981        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2982    elif laue in ['2/m']:
2983        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2984        if uniq == 'a':
2985            dDijDict[phfx+'D23'] = k*l
2986        elif uniq == 'b':
2987            dDijDict[phfx+'D13'] = h*l
2988        elif uniq == 'c':
2989            dDijDict[phfx+'D12'] = h*k
2990    else:
2991        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2992            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2993    if 'C' in calcControls[hfx+'histType']:
2994        for item in dDijDict:
2995            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2996    else:
2997        for item in dDijDict:
2998            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2999    return dDijDict
3000   
3001def GetDij(phfx,SGData,parmDict):
3002    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
3003    return G2spc.HStrainVals(HSvals,SGData)
3004               
3005def GetFobsSq(Histograms,Phases,parmDict,calcControls):
3006    '''Compute the observed structure factors for Powder histograms and store in reflection array
3007    Multiprocessing support added
3008    '''
3009    if GSASIIpath.GetConfigValue('Show_timing',False):
3010        starttime = time.time() #; print 'start GetFobsSq'
3011    histoList = list(Histograms.keys())
3012    histoList.sort()
3013    Ka2 = shl = lamRatio = kRatio = None
3014    for histogram in histoList:
3015        if 'PWDR' in histogram[:4]:
3016            Histogram = Histograms[histogram]
3017            hId = Histogram['hId']
3018            hfx = ':%d:'%(hId)
3019            Limits = calcControls[hfx+'Limits']
3020            if 'C' in calcControls[hfx+'histType']:
3021                shl = max(parmDict[hfx+'SH/L'],0.0005)
3022                Ka2 = False
3023                kRatio = 0.0
3024                if hfx+'Lam1' in list(parmDict.keys()):
3025                    Ka2 = True
3026                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3027                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
3028            x,y,w,yc,yb,yd = Histogram['Data']
3029            xMask = ma.getmaskarray(x)
3030            xB = np.searchsorted(x,Limits[0])
3031            xF = np.searchsorted(x,Limits[1])
3032            ymb = np.array(y-yb)
3033            ymb = np.where(ymb,ymb,1.0)
3034            ycmb = np.array(yc-yb)
3035            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3036            refLists = Histogram['Reflection Lists']
3037            for phase in refLists:
3038                if phase not in Phases:     #skips deleted or renamed phases silently!
3039                    continue
3040                Phase = Phases[phase]
3041                if histogram not in Phase['Histograms']:
3042                    continue
3043                im = 0
3044                if Phase['General'].get('Modulated',False):
3045                    im = 1
3046                pId = Phase['pId']
3047                phfx = '%d:%d:'%(pId,hId)
3048                refDict = refLists[phase]
3049                sumFo = 0.0
3050                sumdF = 0.0
3051                sumFosq = 0.0
3052                sumdFsq = 0.0
3053                sumInt = 0.0
3054                nExcl = 0
3055                # test to see if we are using multiprocessing below
3056                useMP,ncores = G2mp.InitMP()
3057                if len(refDict['RefList']) < 100: useMP = False       
3058                if useMP: # multiprocessing: create a set of initialized Python processes
3059                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3060                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3061                    profArgs = [[] for i in range(G2mp.ncores)]
3062                else:
3063                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3064                if 'C' in calcControls[hfx+'histType']:
3065                    # are we multiprocessing?
3066                    for iref,refl in enumerate(refDict['RefList']):
3067                        if useMP: 
3068                            profArgs[iref%G2mp.ncores].append((refl,iref))
3069                        else:
3070                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3071                            if type(icod) is tuple:
3072                                refl[8+im] = icod[0]
3073                                sumInt += icod[1]
3074                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3075                            elif icod == -1:
3076                                refl[3+im] *= -1
3077                                nExcl += 1
3078                            elif icod == -2:
3079                                break
3080                    if useMP:
3081                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3082                            sumInt += sInt
3083                            for refl8im,irefl in resList:
3084                                if refl8im is None:
3085                                    refDict['RefList'][irefl][3+im] *= -1
3086                                    nExcl += 1
3087                                else:
3088                                    refDict['RefList'][irefl][8+im] = refl8im
3089                                    if parmDict[phfx+'LeBail']:
3090                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3091                elif 'T' in calcControls[hfx+'histType']:
3092                    for iref,refl in enumerate(refDict['RefList']):
3093                        if useMP: 
3094                            profArgs[iref%G2mp.ncores].append((refl,iref))
3095                        else:
3096                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3097                            if type(icod) is tuple:
3098                                refl[8+im] = icod[0]
3099                                sumInt += icod[1]
3100                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3101                            elif icod == -1:
3102                                refl[3+im] *= -1
3103                                nExcl += 1
3104                            elif icod == -2:
3105                                break
3106                    if useMP:
3107                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3108                            sumInt += sInt
3109                            for refl8im,irefl in resList:
3110                                if refl8im is None:
3111                                    refDict['RefList'][irefl][3+im] *= -1
3112                                    nExcl += 1
3113                                else:
3114                                    refDict['RefList'][irefl][8+im] = refl8im
3115                                    if parmDict[phfx+'LeBail']:
3116                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3117                if useMP: MPpool.terminate()
3118                sumFo = 0.0
3119                sumdF = 0.0
3120                sumFosq = 0.0
3121                sumdFsq = 0.0
3122                for iref,refl in enumerate(refDict['RefList']):
3123                    Fo = np.sqrt(np.abs(refl[8+im]))
3124                    Fc = np.sqrt(np.abs(refl[9]+im))
3125                    sumFo += Fo
3126                    sumFosq += refl[8+im]**2
3127                    sumdF += np.abs(Fo-Fc)
3128                    sumdFsq += (refl[8+im]-refl[9+im])**2
3129                if sumFo:
3130                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3131                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3132                else:
3133                    Histogram['Residuals'][phfx+'Rf'] = 100.
3134                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3135                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3136                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3137                Histogram['Residuals']['hId'] = hId
3138        elif 'HKLF' in histogram[:4]:
3139            Histogram = Histograms[histogram]
3140            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3141    if GSASIIpath.GetConfigValue('Show_timing',False):
3142        print ('GetFobsSq t=',time.time()-starttime)
3143               
3144def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup,histogram=None):
3145    'Computes the powder pattern for a histogram based on contributions from all used phases'
3146    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3147   
3148    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3149        U = parmDict[hfx+'U']
3150        V = parmDict[hfx+'V']
3151        W = parmDict[hfx+'W']
3152        X = parmDict[hfx+'X']
3153        Y = parmDict[hfx+'Y']
3154        Z = parmDict[hfx+'Z']
3155        tanPos = tand(refl[5+im]/2.0)
3156        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3157        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3158        sig = max(0.001,sig)
3159        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3160        gam = max(0.001,gam)
3161        return sig,gam
3162               
3163    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3164        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3165            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3166        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3167        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3168        sig += Ssig
3169        gam += Sgam
3170        return sig,gam
3171       
3172    def GetReflAlpBet(refl,im,hfx,parmDict):
3173        alp = parmDict[hfx+'alpha']/refl[4+im]
3174        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3175        return alp,bet
3176       
3177    hId = Histogram['hId']
3178    hfx = ':%d:'%(hId)
3179    bakType = calcControls[hfx+'bakType']
3180    fixedBkg = {i:Histogram['Background'][1][i] for i in Histogram['Background'][1] if i.startswith("_")} 
3181    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixedBkg)
3182    yc = np.zeros_like(yb)
3183    cw = np.diff(ma.getdata(x))
3184    cw = np.append(cw,cw[-1])
3185       
3186    if 'C' in calcControls[hfx+'histType']:   
3187        shl = max(parmDict[hfx+'SH/L'],0.002)
3188        Ka2 = False
3189        if hfx+'Lam1' in (parmDict.keys()):
3190            wave = parmDict[hfx+'Lam1']
3191            Ka2 = True
3192            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3193            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3194        else:
3195            wave = parmDict[hfx+'Lam']
3196    else:
3197        shl = 0.
3198    for phase in Histogram['Reflection Lists']:
3199        refDict = Histogram['Reflection Lists'][phase]
3200        if phase not in Phases:     #skips deleted or renamed phases silently!
3201            continue
3202        Phase = Phases[phase]
3203        if histogram and not histogram in Phase['Histograms']:
3204            continue
3205        pId = Phase['pId']
3206        pfx = '%d::'%(pId)
3207        phfx = '%d:%d:'%(pId,hId)
3208        hfx = ':%d:'%(hId)
3209        SGData = Phase['General']['SGData']
3210        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3211        im = 0
3212        if Phase['General'].get('Modulated',False):
3213            SSGData = Phase['General']['SSGData']
3214            im = 1  #offset in SS reflection list
3215            #??
3216        Dij = GetDij(phfx,SGData,parmDict)
3217        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do someting if Dij << 0.
3218        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3219        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3220            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3221        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3222        Vst = np.sqrt(nl.det(G))    #V*
3223        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3224            if im:
3225                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3226            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3227                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3228            else:
3229                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3230        badPeak = False
3231        # test to see if we are using multiprocessing here
3232        useMP,ncores = G2mp.InitMP()
3233        if len(refDict['RefList']) < 100: useMP = False       
3234        if useMP: # multiprocessing: create a set of initialized Python processes
3235            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3236            profArgs = [[] for i in range(ncores)]
3237        if 'C' in calcControls[hfx+'histType']:
3238            for iref,refl in enumerate(refDict['RefList']):
3239                if im:
3240                    h,k,l,m = refl[:4]
3241                else:
3242                    h,k,l = refl[:3]
3243                Uniq = np.inner(refl[:3],SGMT)
3244                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
3245                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3246                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3247                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3248                refl[11+im] *= Vst*Lorenz
3249                 
3250                if Phase['General'].get('doPawley'):
3251                    try:
3252                        if im:
3253                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3254                        else:
3255                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3256                        refl[9+im] = parmDict[pInd]
3257                    except KeyError:
3258#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3259                        continue
3260                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3261                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3262                iFin = np.searchsorted(x,refl[5+im]+fmax)
3263                if not iBeg+iFin:       #peak below low limit - skip peak
3264                    continue
3265                elif not iBeg-iFin:     #peak above high limit - done
3266                    break
3267                elif iBeg > iFin:   #bad peak coeff - skip
3268                    badPeak = True
3269                    continue
3270                if useMP:
3271                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3272                else:
3273                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
3274                if Ka2:
3275                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3276                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3277                    iBeg = np.searchsorted(x,pos2-fmin)
3278                    iFin = np.searchsorted(x,pos2+fmax)
3279                    if not iBeg+iFin:       #peak below low limit - skip peak
3280                        continue
3281                    elif not iBeg-iFin:     #peak above high limit - done
3282                        return yc,yb
3283                    elif iBeg > iFin:   #bad peak coeff - skip
3284                        continue
3285                    if useMP:
3286                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3287                    else:
3288                        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
3289        elif 'T' in calcControls[hfx+'histType']:
3290            for iref,refl in enumerate(refDict['RefList']):
3291                if im:
3292                    h,k,l,m = refl[:4]
3293                else:
3294                    h,k,l = refl[:3]
3295                Uniq = np.inner(refl[:3],SGMT)
3296                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
3297                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
3298#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
3299                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3300                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3301                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3302                refl[11+im] *= Vst*Lorenz
3303                if Phase['General'].get('doPawley'):
3304                    try:
3305                        if im:
3306                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3307                        else:
3308                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3309                        refl[9+im] = parmDict[pInd]
3310                    except KeyError:
3311#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3312                        continue
3313                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
3314                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3315                iFin = np.searchsorted(x,refl[5+im]+fmax)
3316                if not iBeg+iFin:       #peak below low limit - skip peak
3317                    continue
3318                elif not iBeg-iFin:     #peak above high limit - done
3319                    break
3320                elif iBeg > iFin:   #bad peak coeff - skip
3321                    badPeak = True
3322                    continue
3323                if useMP:
3324                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin))
3325                else:
3326                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
3327#        print 'profile calc time: %.3fs'%(time.time()-time0)
3328        if useMP and 'C' in calcControls[hfx+'histType']:
3329            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
3330                yc += y
3331            MPpool.terminate()
3332        elif useMP:
3333            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
3334                yc += y
3335            MPpool.terminate()
3336    if badPeak:
3337        print ('ouch #4 bad profile coefficients yield negative peak width; some reflections skipped')
3338    if GSASIIpath.GetConfigValue('Show_timing',False):
3339        print ('getPowderProfile t=%.3f'%(time.time()-starttime))
3340    return yc,yb
3341   
3342def getPowderProfileDervMP(args):
3343    '''Computes the derivatives of the computed powder pattern with respect to all
3344    refined parameters.
3345    Multiprocessing version.
3346    '''
3347    import pytexture as ptx
3348    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics for each processor
3349    parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars = args[:9]
3350    prc,tprc,histogram = 0,1,None
3351    if len(args) >= 10: prc=args[9]
3352    if len(args) >= 11: tprc=args[10]
3353    if len(args) >= 12: his