source: trunk/GSASIIstrMath.py @ 4535

Last change on this file since 4535 was 4535, checked in by vondreele, 3 years ago

fix X-label for pink powder patterns
revert incomm. mag. str. factor calcs to previous best one - still not correct tho.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 224.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2020-07-31 14:04:06 +0000 (Fri, 31 Jul 2020) $
8# $Author: vondreele $
9# $Revision: 4535 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4535 2020-07-31 14:04:06Z vondreele $
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: 4535 $")
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'] or 'NB' 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'] or 'NB' 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'] or 'NB' 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    Nops = len(SGMT)*(1+SGData['SGInv'])
1493    if SGData['SGGray']:
1494        Nops *= 2
1495    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1496    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1497    SSCen = SSGData['SSGCen']
1498    FFtables = calcControls['FFtables']
1499    BLtables = calcControls['BLtables']
1500    MFtables = calcControls['MFtables']
1501    Amat,Bmat = G2lat.Gmat2AB(G)
1502    Flack = 1.0
1503    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1504        Flack = 1.-2.*parmDict[phfx+'Flack']
1505    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1506        GetAtomFXU(pfx,calcControls,parmDict)
1507    if not Xdata.size:          #no atoms in phase!
1508        return
1509    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1510    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1511    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1512
1513    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms - Mmod matched drawing & Bilbao drawings
1514   
1515        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
1516        if SGData['SGGray']:
1517            mXYZ = np.hstack((mXYZ,mXYZ))
1518
1519        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sim parts sum matches drawing
1520       
1521        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1522            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1523            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1524                GSdata = np.hstack((GSdata,-GSdata))     
1525            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1526            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1527            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1528       
1529    FF = np.zeros(len(Tdata))
1530    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1531        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1532    elif 'X' in calcControls[hfx+'histType']:
1533        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1534        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1535    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1536    bij = Mast*Uij
1537    blkSize = 48       #no. of reflections in a block
1538    nRef = refDict['RefList'].shape[0]
1539    SQ = 1./(2.*refDict['RefList'].T[5])**2
1540    if 'N' in calcControls[hfx+'histType']:
1541        dat = G2el.getBLvalues(BLtables)
1542        refDict['FF']['El'] = list(dat.keys())
1543        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1544        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1545        for iel,El in enumerate(refDict['FF']['El']):
1546            if El in MFtables:
1547                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1548    else:
1549        dat = G2el.getFFvalues(FFtables,0.)
1550        refDict['FF']['El'] = list(dat.keys())
1551        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1552        for iel,El in enumerate(refDict['FF']['El']):
1553            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1554#reflection processing begins here - big arrays!
1555    iBeg = 0
1556    while iBeg < nRef:
1557        iFin = min(iBeg+blkSize,nRef)
1558        mRef= iFin-iBeg
1559        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1560        H = refl.T[:4]                          #array(blkSize,4)
1561        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1562        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1563        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1564        Uniq = np.inner(H.T,SSGMT)
1565        UniqP = np.inner(HP.T,SGMT)
1566        Phi = np.inner(H.T,SSGT)
1567        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1568            Uniq = np.hstack((Uniq,-Uniq))
1569            Phi = np.hstack((Phi,-Phi))
1570            UniqP = np.hstack((UniqP,-UniqP))
1571        if 'T' in calcControls[hfx+'histType']:
1572            if 'P' in calcControls[hfx+'histType']:
1573                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1574            else:
1575                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1576            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1577            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1578        Bab = 0.
1579        if phfx+'BabA' in parmDict:
1580            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1581        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1582        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1583        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1584        phase = np.hstack([phase for cen in SSCen])
1585        sinp = np.sin(phase)
1586        cosp = np.cos(phase)
1587        biso = -SQfactor*Uisodata[:,nxs]
1588        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1589        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1590        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1591        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1592
1593        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:           
1594           
1595            phasem = twopi*np.inner(mXYZ,HP.T).T    #2pi(Q.r)
1596            cosm = np.cos(phasem)                   #Nref,nops,natm
1597            sinm = np.sin(phasem)
1598            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1599            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2*Nops)     #Nref,Natm
1600            HM = np.inner(uAmat,HP.T)                    #put into cartesian space X||H,Z||H*L; uAmat better than uBmat
1601            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T    # normalize  HP  Nref,hkl=Unit vectors || Q
1602
1603            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1604                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1605                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]
1606#  calc mag. structure factors; Nref,Ntau,Nops,Natm,Mxyz                           
1607            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]+   
1608                H[3,i]*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1609                       
1610            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+   
1611                H[3,i]*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1612
1613            if not SGData['SGGray']:           
1614                fams += fam0[:,nxs,:,:,:]
1615                fbms += fbm0[:,nxs,:,:,:]
1616                               
1617#sum ops & atms                               
1618            fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Mxyz; sum ops & atoms
1619            fbsm = np.sum(np.sum(fbms,axis=-2),axis=-2)
1620#put into cartesian space
1621            facm = np.inner(fasm,uBmat.T)
1622            fbcm = np.inner(fbsm,uBmat.T)
1623#form e.F dot product
1624            eDotFa = np.sum(eM[:,nxs,:]*facm,axis=-1)    #Nref,Ntau       
1625            eDotFb = np.sum(eM[:,nxs,:]*fbcm,axis=-1)
1626#intensity Halpern & Johnson
1627            fass = np.sum((facm-eM[:,nxs,:]*eDotFa[:,:,nxs])**2,axis=-1)
1628            fbss = np.sum((fbcm-eM[:,nxs,:]*eDotFb[:,:,nxs])**2,axis=-1)
1629## #do integration           
1630            fas = np.sum(fass*glWt[nxs,:],axis=1)
1631            fbs = np.sum(fbss*glWt[nxs,:],axis=1)
1632           
1633            refl.T[10] = fas+fbs   #Sum(fams**2,Mxyz) Re + Im
1634            refl.T[11] = atan2d(fbs,fas)
1635
1636        else:
1637            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1638            if 'T' in calcControls[hfx+'histType']:
1639                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1640                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1641            else:
1642                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1643                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1644            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1645            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1646            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1647            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1648           
1649            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1650            refl.T[11] = atan2d(fbs[0],fas[0])  #use only tau=0; ignore f' & f"
1651        if 'P' not in calcControls[hfx+'histType']:
1652            refl.T[8] = np.copy(refl.T[10])               
1653        iBeg += blkSize
1654    return copy.deepcopy(refDict['RefList'])
1655
1656def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1657    '''
1658    Compute super structure factors for all h,k,l,m for phase - twins only
1659    puts the result, F^2, in each ref[8+im] in refList
1660    works on blocks of 32 reflections for speed
1661    input:
1662   
1663    :param dict refDict: where
1664        'RefList' list where each ref = h,k,l,m,it,d,...
1665        'FF' dict of form factors - filed in below
1666    :param np.array G:      reciprocal metric tensor
1667    :param str pfx:    phase id string
1668    :param dict SGData: space group info. dictionary output from SpcGroup
1669    :param dict calcControls:
1670    :param dict ParmDict:
1671
1672    '''
1673    phfx = pfx.split(':')[0]+hfx
1674    ast = np.sqrt(np.diag(G))
1675    Mast = twopisq*np.multiply.outer(ast,ast)   
1676    SGInv = SGData['SGInv']
1677    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1678    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1679    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1680    FFtables = calcControls['FFtables']
1681    BLtables = calcControls['BLtables']
1682    MFtables = calcControls['MFtables']
1683    Flack = 1.0
1684    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1685        Flack = 1.-2.*parmDict[phfx+'Flack']
1686    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1687    TwDict = refDict.get('TwDict',{})           
1688    if 'S' in calcControls[hfx+'histType']:
1689        NTL = calcControls[phfx+'NTL']
1690        NM = calcControls[phfx+'TwinNMN']+1
1691        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1692        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1693        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1694    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1695        GetAtomFXU(pfx,calcControls,parmDict)
1696    if not Xdata.size:          #no atoms in phase!
1697        return
1698    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1699    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1700    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1701    FF = np.zeros(len(Tdata))
1702    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1703        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1704    elif 'X' in calcControls[hfx+'histType']:
1705        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1706        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1707    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1708    bij = Mast*Uij
1709    blkSize = 32       #no. of reflections in a block
1710    nRef = refDict['RefList'].shape[0]
1711    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1712        SQ = 1./(2.*refDict['RefList'].T[5])**2
1713        if 'N' in calcControls[hfx+'histType']:
1714            dat = G2el.getBLvalues(BLtables)
1715            refDict['FF']['El'] = list(dat.keys())
1716            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1717            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1718            for iel,El in enumerate(refDict['FF']['El']):
1719                if El in MFtables:
1720                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1721        else:
1722            dat = G2el.getFFvalues(FFtables,0.)
1723            refDict['FF']['El'] = list(dat.keys())
1724            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1725            for iel,El in enumerate(refDict['FF']['El']):
1726                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1727#    time0 = time.time()
1728#reflection processing begins here - big arrays!
1729    iBeg = 0
1730    while iBeg < nRef:
1731        iFin = min(iBeg+blkSize,nRef)
1732        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1733        H = refl[:,:4]                          #array(blkSize,4)
1734        H3 = refl[:,:3]
1735        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1736        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1737        H3 = np.inner(H3,TwinLaw)       
1738        TwMask = np.any(HP,axis=-1)
1739        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1740            for ir in range(blkSize):
1741                iref = ir+iBeg
1742                if iref in TwDict:
1743                    for i in TwDict[iref]:
1744                        for n in range(NTL):
1745                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1746                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1747            TwMask = np.any(HP,axis=-1)
1748        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1749        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1750        Uniq = np.inner(H,SSGMT)
1751        Uniq3 = np.inner(H3,SGMT)
1752        UniqP = np.inner(HP,SGMT)
1753        Phi = np.inner(H,SSGT)
1754        if SGInv:   #if centro - expand HKL sets
1755            Uniq = np.hstack((Uniq,-Uniq))
1756            Uniq3 = np.hstack((Uniq3,-Uniq3))
1757            Phi = np.hstack((Phi,-Phi))
1758            UniqP = np.hstack((UniqP,-UniqP))
1759        if 'T' in calcControls[hfx+'histType']:
1760            if 'P' in calcControls[hfx+'histType']:
1761                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1762            else:
1763                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1764            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1765            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1766        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1767        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1768        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1769        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1770        sinp = np.sin(phase)
1771        cosp = np.cos(phase)
1772        biso = -SQfactor*Uisodata[:,nxs]
1773        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1774        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1775        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1776        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1777        if 'T' in calcControls[hfx+'histType']:
1778            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1779            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1780        else:
1781            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1782            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1783        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1784        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1785        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1786        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1787        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1788        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1789        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1790            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1791        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1792        iBeg += blkSize
1793#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1794
1795def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1796    '''
1797    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1798    Only Fourier component are done analytically here
1799    input:
1800   
1801    :param dict refDict: where
1802        'RefList' list where each ref = h,k,l,m,it,d,...
1803        'FF' dict of form factors - filled in below
1804    :param int im: = 1 (could be eliminated)
1805    :param np.array G:      reciprocal metric tensor
1806    :param str hfx:    histogram id string
1807    :param str pfx:    phase id string
1808    :param dict SGData: space group info. dictionary output from SpcGroup
1809    :param dict SSGData: super space group info.
1810    :param dict calcControls:
1811    :param dict ParmDict:
1812   
1813    :returns: dict dFdvDict: dictionary of derivatives
1814    '''
1815    phfx = pfx.split(':')[0]+hfx
1816    ast = np.sqrt(np.diag(G))
1817    Mast = twopisq*np.multiply.outer(ast,ast)
1818    SGInv = SGData['SGInv']
1819    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1820    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1821    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1822    FFtables = calcControls['FFtables']
1823    BLtables = calcControls['BLtables']
1824    nRef = len(refDict['RefList'])
1825    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1826        GetAtomFXU(pfx,calcControls,parmDict)
1827    if not Xdata.size:          #no atoms in phase!
1828        return {}
1829    mSize = len(Mdata)  #no. atoms
1830    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1831    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1832    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1833    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1834    FF = np.zeros(len(Tdata))
1835    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
1836        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1837    elif 'X' in calcControls[hfx+'histType']:
1838        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1839        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1840    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1841    bij = Mast*Uij
1842    if not len(refDict['FF']):
1843        if 'N' in calcControls[hfx+'histType']:
1844            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1845        else:
1846            dat = G2el.getFFvalues(FFtables,0.)       
1847        refDict['FF']['El'] = list(dat.keys())
1848        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1849    dFdvDict = {}
1850    dFdfr = np.zeros((nRef,mSize))
1851    dFdx = np.zeros((nRef,mSize,3))
1852    dFdui = np.zeros((nRef,mSize))
1853    dFdua = np.zeros((nRef,mSize,6))
1854    dFdbab = np.zeros((nRef,2))
1855    dFdfl = np.zeros((nRef))
1856    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1857    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1858    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1859    Flack = 1.0
1860    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1861        Flack = 1.-2.*parmDict[phfx+'Flack']
1862    time0 = time.time()
1863    nRef = len(refDict['RefList'])/100
1864    for iref,refl in enumerate(refDict['RefList']):
1865        if 'T' in calcControls[hfx+'histType']:
1866            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1867        H = np.array(refl[:4])
1868        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1869        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1870        SQfactor = 8.0*SQ*np.pi**2
1871        Bab = 0.0
1872        if phfx+'BabA' in parmDict:
1873            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1874            Bab = parmDict[phfx+'BabA']*dBabdA
1875        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1876        FF = refDict['FF']['FF'][iref].T[Tindx]
1877        Uniq = np.inner(H,SSGMT)
1878        Phi = np.inner(H,SSGT)
1879        UniqP = np.inner(HP,SGMT)
1880        if SGInv:   #if centro - expand HKL sets
1881            Uniq = np.vstack((Uniq,-Uniq))
1882            Phi = np.hstack((Phi,-Phi))
1883            UniqP = np.vstack((UniqP,-UniqP))
1884        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1885        sinp = np.sin(phase)
1886        cosp = np.cos(phase)
1887        occ = Mdata*Fdata/Uniq.shape[0]
1888        biso = -SQfactor*Uisodata[:,nxs]
1889        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1890        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1891        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1892        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1893        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1894        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1895        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1896        fotp = FPP*Tcorr            #ops x atoms
1897        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1898        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1899        # GfpuA is 2 x ops x atoms
1900        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1901        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1902        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1903        fag = fa*GfpuA[0]-fb*GfpuA[1]
1904        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1905       
1906        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1907        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1908        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1909        fbx = np.array([fot*cosp,-fotp*sinp])
1910        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1911        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1912        #sum below is over Uniq
1913        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1914        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1915        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1916        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1917        dfadui = np.sum(-SQfactor*fag,axis=1)
1918        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1919        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1920        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1921        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1922        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1923        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1924        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1925        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1926        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1927        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1928        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1929        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1930        if not SGData['SGInv']:   #Flack derivative
1931            dfadfl = np.sum(-FPP*Tcorr*sinp)
1932            dfbdfl = np.sum(FPP*Tcorr*cosp)
1933        else:
1934            dfadfl = 1.0
1935            dfbdfl = 1.0
1936        SA = fas[0]+fas[1]      #float = A+A'
1937        SB = fbs[0]+fbs[1]      #float = B+B'
1938        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1939            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1940            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1941                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1942            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1943                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1944            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1945                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1946            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1947                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1948            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1949                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1950            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1951                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1952            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1953                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1954        else:                       #OK, I think
1955            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1956            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1957            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1958            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1959            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1960                           
1961            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1962            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1963            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1964        if phfx+'BabA' in parmDict:
1965            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1966                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1967        #loop over atoms - each dict entry is list of derivatives for all the reflections
1968        if not iref%100 :
1969            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1970    for i in range(len(Mdata)):     #loop over atoms
1971        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1972        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1973        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1974        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1975        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1976        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1977        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1978        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1979        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1980        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1981        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1982        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1983            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1984            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1985        nx = 0
1986        if waveTypes[i] in ['Block','ZigZag']:
1987            nx = 1 
1988        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1989            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1990            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1991            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1992            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1993            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1994            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1995        for j in range(USSdata.shape[1]):       #loop over waves
1996            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1997            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1998            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1999            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2000            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2001            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2002            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2003            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2004            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2005            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2006            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2007            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2008           
2009    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
2010    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2011    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2012    return dFdvDict
2013
2014def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2015    '''
2016    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
2017    input:
2018   
2019    :param dict refDict: where
2020        'RefList' list where each ref = h,k,l,m,it,d,...
2021        'FF' dict of form factors - filled in below
2022    :param int im: = 1 (could be eliminated)
2023    :param np.array G:      reciprocal metric tensor
2024    :param str hfx:    histogram id string
2025    :param str pfx:    phase id string
2026    :param dict SGData: space group info. dictionary output from SpcGroup
2027    :param dict SSGData: super space group info.
2028    :param dict calcControls:
2029    :param dict ParmDict:
2030   
2031    :returns: dict dFdvDict: dictionary of derivatives
2032    '''
2033
2034    trefDict = copy.deepcopy(refDict)
2035    dM = 1.e-4
2036    dFdvDict = {}
2037    for parm in parmDict:
2038        if ':' not in parm:
2039            continue
2040        if parm.split(':')[2] in ['Tmin','Tmax','Xmax','Ymax','Zmax','Fzero','Fwid',
2041            'MXsin','MXcos','MYsin','MYcos','MZsin','MZcos','AMx','AMy','AMz',]:
2042            parmDict[parm] += dM
2043            prefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2044            parmDict[parm] -= 2*dM
2045            mrefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2046            parmDict[parm] += dM
2047            dFdvDict[parm] = (prefList[:,9+im]-mrefList[:,9+im])/(2.*dM)
2048    return dFdvDict
2049   
2050def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2051    'Needs a doc string'
2052    phfx = pfx.split(':')[0]+hfx
2053    ast = np.sqrt(np.diag(G))
2054    Mast = twopisq*np.multiply.outer(ast,ast)
2055    SGInv = SGData['SGInv']
2056    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2057    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2058    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2059    FFtables = calcControls['FFtables']
2060    BLtables = calcControls['BLtables']
2061    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2062    TwDict = refDict.get('TwDict',{})           
2063    if 'S' in calcControls[hfx+'histType']:
2064        NTL = calcControls[phfx+'NTL']
2065        NM = calcControls[phfx+'TwinNMN']+1
2066        TwinLaw = calcControls[phfx+'TwinLaw']
2067        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2068    nTwin = len(TwinLaw)       
2069    nRef = len(refDict['RefList'])
2070    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2071        GetAtomFXU(pfx,calcControls,parmDict)
2072    if not Xdata.size:          #no atoms in phase!
2073        return {}
2074    mSize = len(Mdata)  #no. atoms
2075    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2076    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
2077    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2078    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2079    FF = np.zeros(len(Tdata))
2080    if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
2081        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2082    elif 'X' in calcControls[hfx+'histType']:
2083        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2084        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2085    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2086    bij = Mast*Uij
2087    if not len(refDict['FF']):
2088        if 'N' in calcControls[hfx+'histType']:
2089            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2090        else:
2091            dat = G2el.getFFvalues(FFtables,0.)       
2092        refDict['FF']['El'] = list(dat.keys())
2093        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2094    dFdvDict = {}
2095    dFdfr = np.zeros((nRef,nTwin,mSize))
2096    dFdx = np.zeros((nRef,nTwin,mSize,3))
2097    dFdui = np.zeros((nRef,nTwin,mSize))
2098    dFdua = np.zeros((nRef,nTwin,mSize,6))
2099    dFdbab = np.zeros((nRef,nTwin,2))
2100    dFdtw = np.zeros((nRef,nTwin))
2101    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2102    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2103    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2104    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2105    Flack = 1.0
2106    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2107        Flack = 1.-2.*parmDict[phfx+'Flack']
2108    time0 = time.time()
2109    nRef = len(refDict['RefList'])/100
2110    for iref,refl in enumerate(refDict['RefList']):
2111        if 'T' in calcControls[hfx+'histType']:
2112            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2113        H = np.array(refl[:4])
2114        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2115        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2116        TwMask = np.any(H,axis=-1)
2117        if TwinLaw.shape[0] > 1 and TwDict:
2118            if iref in TwDict:
2119                for i in TwDict[iref]:
2120                    for n in range(NTL):
2121                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2122            TwMask = np.any(H,axis=-1)
2123        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2124        SQfactor = 8.0*SQ*np.pi**2
2125        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2126        Bab = parmDict[phfx+'BabA']*dBabdA
2127        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2128        FF = refDict['FF']['FF'][iref].T[Tindx]
2129        Uniq = np.inner(H,SSGMT)
2130        Phi = np.inner(H,SSGT)
2131        UniqP = np.inner(HP,SGMT)
2132        if SGInv:   #if centro - expand HKL sets
2133            Uniq = np.vstack((Uniq,-Uniq))
2134            Phi = np.hstack((Phi,-Phi))
2135            UniqP = np.vstack((UniqP,-UniqP))
2136        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2137        sinp = np.sin(phase)
2138        cosp = np.cos(phase)
2139        occ = Mdata*Fdata/Uniq.shape[0]
2140        biso = -SQfactor*Uisodata[:,nxs]
2141        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2142        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2143        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2144        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2145        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2146        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2147        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2148        fotp = FPP*Tcorr            #ops x atoms
2149        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2150        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2151        # GfpuA is 2 x ops x atoms
2152        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2153        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2154        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2155        fag = fa*GfpuA[0]-fb*GfpuA[1]
2156        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2157       
2158        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2159        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2160        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2161        fbx = np.array([fot*cosp,-fotp*sinp])
2162        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2163        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2164        #sum below is over Uniq
2165        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2166        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2167        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2168        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2169        dfadui = np.sum(-SQfactor*fag,axis=1)
2170        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2171        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2172        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2173        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2174        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2175        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2176        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2177        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2178        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2179        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2180        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2181        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2182        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2183        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2184#        GSASIIpath.IPyBreak()
2185        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2186        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2187        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2188        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)]
2189        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)]
2190        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)]
2191        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)]
2192        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2193
2194        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2195        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2196        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2197        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2198#            GSASIIpath.IPyBreak()
2199        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2200            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2201        #loop over atoms - each dict entry is list of derivatives for all the reflections
2202        if not iref%100 :
2203            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2204    for i in range(len(Mdata)):     #loop over atoms
2205        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2206        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2207        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2208        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2209        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2210        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2211        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2212        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2213        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2214        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2215        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2216        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2217            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2218            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2219        nx = 0
2220        if waveTypes[i] in ['Block','ZigZag']:
2221            nx = 1 
2222            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2223            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2224            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2225            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2226            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2227        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2228            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2229            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2230            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2231            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2232            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2233            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2234        for j in range(USSdata.shape[1]):       #loop over waves
2235            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2236            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2237            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2238            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2239            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2240            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2241            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2242            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2243            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2244            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2245            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2246            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2247           
2248#        GSASIIpath.IPyBreak()
2249    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2250    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2251    return dFdvDict
2252   
2253def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2254    ''' Single crystal extinction function; returns extinction & derivative
2255    '''
2256    extCor = 1.0
2257    dervDict = {}
2258    dervCor = 1.0
2259    if calcControls[phfx+'EType'] != 'None':
2260        SQ = 1/(4.*ref[4+im]**2)
2261        if 'C' in parmDict[hfx+'Type']:           
2262            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2263        else:   #'T'
2264            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2265        if 'SXC' in parmDict[hfx+'Type']:
2266            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2267            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2268            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2269            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2270        elif 'SNT' in parmDict[hfx+'Type']:
2271            AV = 1.e7/parmDict[pfx+'Vol']**2
2272            PL = SQ
2273            PLZ = AV*ref[9+im]*ref[12+im]**2
2274        elif 'SNC' in parmDict[hfx+'Type']:
2275            AV = 1.e7/parmDict[pfx+'Vol']**2
2276            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2277            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2278           
2279        if 'Primary' in calcControls[phfx+'EType']:
2280            PLZ *= 1.5
2281        else:
2282            if 'C' in parmDict[hfx+'Type']:
2283                PLZ *= calcControls[phfx+'Tbar']
2284            else: #'T'
2285                PLZ *= ref[13+im]      #t-bar
2286        if 'Primary' in calcControls[phfx+'EType']:
2287            PLZ *= 1.5
2288            PSIG = parmDict[phfx+'Ep']
2289        elif 'I & II' in calcControls[phfx+'EType']:
2290            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2291        elif 'Type II' in calcControls[phfx+'EType']:
2292            PSIG = parmDict[phfx+'Es']
2293        else:       # 'Secondary Type I'
2294            PSIG = parmDict[phfx+'Eg']/PL
2295           
2296        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2297        AL = 0.025+0.285*cos2T
2298        BG = 0.02-0.025*cos2T
2299        BL = 0.15-0.2*(0.75-cos2T)**2
2300        if cos2T < 0.:
2301            BL = -0.45*cos2T
2302        CG = 2.
2303        CL = 2.
2304        PF = PLZ*PSIG
2305       
2306        if 'Gaussian' in calcControls[phfx+'EApprox']:
2307            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2308            extCor = np.sqrt(PF4)
2309            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2310        else:
2311            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2312            extCor = np.sqrt(PF4)
2313            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2314
2315        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2316        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2317            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2318        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2319            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2320        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2321            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2322               
2323    return 1./extCor,dervDict,dervCor
2324   
2325def Dict2Values(parmdict, varylist):
2326    '''Use before call to leastsq to setup list of values for the parameters
2327    in parmdict, as selected by key in varylist'''
2328    return [parmdict[key] for key in varylist] 
2329   
2330def Values2Dict(parmdict, varylist, values):
2331    ''' Use after call to leastsq to update the parameter dictionary with
2332    values corresponding to keys in varylist'''
2333    parmdict.update(zip(varylist,values))
2334   
2335def GetNewCellParms(parmDict,varyList):
2336    '''Compute unit cell tensor terms from varied Aij and Dij values.
2337    Terms are included in the dict only if Aij or Dij is varied.
2338    '''
2339    newCellDict = {}
2340    Anames = ['A'+str(i) for i in range(6)]
2341    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2342    for item in varyList:
2343        keys = item.split(':')
2344        if keys[2] in Ddict:
2345            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2346            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2347            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2348    return newCellDict          # is e.g. {'0::D11':A0-D11}
2349   
2350def ApplyXYZshifts(parmDict,varyList):
2351    '''
2352    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2353   
2354    :param dict parmDict: parameter dictionary
2355    :param list varyList: list of variables (not used!)
2356    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2357
2358    '''
2359    newAtomDict = {}
2360    for item in parmDict:
2361        if 'dA' in item:
2362            parm = ''.join(item.split('d'))
2363            parmDict[parm] += parmDict[item]
2364            newAtomDict[item] = [parm,parmDict[parm]]
2365    return newAtomDict
2366   
2367def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2368    'Spherical harmonics texture'
2369    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2370    if 'T' in calcControls[hfx+'histType']:
2371        tth = parmDict[hfx+'2-theta']
2372    else:
2373        tth = refl[5+im]
2374    odfCor = 1.0
2375    H = refl[:3]
2376    cell = G2lat.Gmat2cell(g)
2377    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2378    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2379    phi,beta = G2lat.CrsAng(H,cell,SGData)
2380    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2381    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2382    for item in SHnames:
2383        L,M,N = eval(item.strip('C'))
2384        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2385        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2386        Lnorm = G2lat.Lnorm(L)
2387        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2388    return odfCor
2389   
2390def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2391    'Spherical harmonics texture derivatives'
2392    if 'T' in calcControls[hfx+'histType']:
2393        tth = parmDict[hfx+'2-theta']
2394    else:
2395        tth = refl[5+im]
2396    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2397    odfCor = 1.0
2398    dFdODF = {}
2399    dFdSA = [0,0,0]
2400    H = refl[:3]
2401    cell = G2lat.Gmat2cell(g)
2402    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2403    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2404    phi,beta = G2lat.CrsAng(H,cell,SGData)
2405    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2406    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2407    for item in SHnames:
2408        L,M,N = eval(item.strip('C'))
2409        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2410        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2411        Lnorm = G2lat.Lnorm(L)
2412        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2413        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2414        for i in range(3):
2415            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2416    return odfCor,dFdODF,dFdSA
2417   
2418def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2419    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2420    if 'T' in calcControls[hfx+'histType']:
2421        tth = parmDict[hfx+'2-theta']
2422    else:
2423        tth = refl[5+im]
2424    odfCor = 1.0
2425    H = refl[:3]
2426    cell = G2lat.Gmat2cell(g)
2427    Sangls = [0.,0.,0.]
2428    if 'Bragg' in calcControls[hfx+'instType']:
2429        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2430        IFCoup = True
2431    else:
2432        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2433        IFCoup = False
2434    phi,beta = G2lat.CrsAng(H,cell,SGData)
2435    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2436    SHnames = calcControls[phfx+'SHnames']
2437    for item in SHnames:
2438        L,N = eval(item.strip('C'))
2439        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2440        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2441        Lnorm = G2lat.Lnorm(L)
2442        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2443    return np.squeeze(odfCor)
2444   
2445def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2446    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2447    if 'T' in calcControls[hfx+'histType']:
2448        tth = parmDict[hfx+'2-theta']
2449    else:
2450        tth = refl[5+im]
2451    odfCor = 1.0
2452    dFdODF = {}
2453    H = refl[:3]
2454    cell = G2lat.Gmat2cell(g)
2455    Sangls = [0.,0.,0.]
2456    if 'Bragg' in calcControls[hfx+'instType']:
2457        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2458        IFCoup = True
2459    else:
2460        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2461        IFCoup = False
2462    phi,beta = G2lat.CrsAng(H,cell,SGData)
2463    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2464    SHnames = calcControls[phfx+'SHnames']
2465    for item in SHnames:
2466        L,N = eval(item.strip('C'))
2467        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2468        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2469        Lnorm = G2lat.Lnorm(L)
2470        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2471        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2472    return odfCor,dFdODF
2473   
2474def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2475    'March-Dollase preferred orientation correction'
2476    POcorr = 1.0
2477    MD = parmDict[phfx+'MD']
2478    if MD != 1.0:
2479        MDAxis = calcControls[phfx+'MDAxis']
2480        sumMD = 0
2481        for H in uniq:           
2482            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2483            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2484            sumMD += A**3
2485        POcorr = sumMD/len(uniq)
2486    return POcorr
2487   
2488def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2489    'Needs a doc string'
2490    POcorr = 1.0
2491    POderv = {}
2492    if calcControls[phfx+'poType'] == 'MD':
2493        MD = parmDict[phfx+'MD']
2494        MDAxis = calcControls[phfx+'MDAxis']
2495        sumMD = 0
2496        sumdMD = 0
2497        for H in uniq:           
2498            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2499            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2500            sumMD += A**3
2501            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2502        POcorr = sumMD/len(uniq)
2503        POderv[phfx+'MD'] = sumdMD/len(uniq)
2504    else:   #spherical harmonics
2505        if calcControls[phfx+'SHord']:
2506            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2507    return POcorr,POderv
2508   
2509def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2510    'Needs a doc string'
2511    if 'Debye' in calcControls[hfx+'instType']:
2512        if 'T' in calcControls[hfx+'histType']:
2513            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2514        else:
2515            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2516    else:
2517        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2518   
2519def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2520    'Needs a doc string'
2521    if 'Debye' in calcControls[hfx+'instType']:
2522        if 'T' in calcControls[hfx+'histType']:
2523            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2524        else:
2525            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2526    else:
2527        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2528       
2529def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2530    'Needs a doc string'
2531    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2532    pi2 = np.sqrt(2./np.pi)
2533    if 'T' in calcControls[hfx+'histType']:
2534        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2535        wave = refl[14+im]
2536    else:   #'C'W
2537        sth2 = sind(refl[5+im]/2.)**2
2538        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2539    c2th = 1.-2.0*sth2
2540    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2541    if 'X' in calcControls[hfx+'histType']:
2542        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2543    xfac = flv2*parmDict[phfx+'Extinction']
2544    exb = 1.0
2545    if xfac > -1.:
2546        exb = 1./np.sqrt(1.+xfac)
2547    exl = 1.0
2548    if 0 < xfac <= 1.:
2549        xn = np.array([xfac**(i+1) for i in range(6)])
2550        exl += np.sum(xn*coef)
2551    elif xfac > 1.:
2552        xfac2 = 1./np.sqrt(xfac)
2553        exl = pi2*(1.-0.125/xfac)*xfac2
2554    return exb*sth2+exl*(1.-sth2)
2555   
2556def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2557    'Needs a doc string'
2558    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2559    pi2 = np.sqrt(2./np.pi)
2560    if 'T' in calcControls[hfx+'histType']:
2561        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2562        wave = refl[14+im]
2563    else:   #'C'W
2564        sth2 = sind(refl[5+im]/2.)**2
2565        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2566    c2th = 1.-2.0*sth2
2567    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2568    if 'X' in calcControls[hfx+'histType']:
2569        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2570    xfac = flv2*parmDict[phfx+'Extinction']
2571    dbde = -500.*flv2
2572    if xfac > -1.:
2573        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2574    dlde = 0.
2575    if 0 < xfac <= 1.:
2576        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2577        dlde = np.sum(xn*coef)/xfac
2578    elif xfac > 1.:
2579        xfac2 = 1./np.sqrt(xfac)
2580        dlde = 0.5*flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2581       
2582    return dbde*sth2+dlde*(1.-sth2)
2583   
2584def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2585    'Needs a doc string'    #need powder extinction!
2586    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2587    if 'X' in parmDict[hfx+'Type']:
2588        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2589    POcorr = 1.0
2590    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2591        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2592    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2593        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2594    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2595        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2596    Icorr *= POcorr
2597    AbsCorr = 1.0
2598    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2599    Icorr *= AbsCorr
2600    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2601    Icorr *= ExtCorr
2602    return Icorr,POcorr,AbsCorr,ExtCorr
2603   
2604def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2605    'Needs a doc string'    #need powder extinction derivs!
2606    dIdsh = 1./parmDict[hfx+'Scale']
2607    dIdsp = 1./parmDict[phfx+'Scale']
2608    if 'X' in parmDict[hfx+'Type']:
2609        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2610        dIdPola /= pola
2611    else:       #'N'
2612        dIdPola = 0.0
2613    dFdODF = {}
2614    dFdSA = [0,0,0]
2615    dIdPO = {}
2616    if pfx+'SHorder' in parmDict:
2617        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2618        for iSH in dFdODF:
2619            dFdODF[iSH] /= odfCor
2620        for i in range(3):
2621            dFdSA[i] /= odfCor
2622    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2623        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2624        for iPO in dIdPO:
2625            dIdPO[iPO] /= POcorr
2626    if 'T' in parmDict[hfx+'Type']:
2627        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2628        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2629    else:
2630        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2631        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2632    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2633       
2634def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2635    'Needs a doc string'
2636    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:     #All checked & OK
2637        costh = cosd(refl[5+im]/2.)
2638        #crystallite size
2639        if calcControls[phfx+'SizeType'] == 'isotropic':
2640            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2641        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2642            H = np.array(refl[:3])
2643            P = np.array(calcControls[phfx+'SizeAxis'])
2644            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2645            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2646            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2647        else:           #ellipsoidal crystallites
2648            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2649            H = np.array(refl[:3])
2650            lenR = G2pwd.ellipseSize(H,Sij,GB)
2651            Sgam = 1.8*wave/(np.pi*costh*lenR)
2652        #microstrain               
2653        if calcControls[phfx+'MustrainType'] == 'isotropic':
2654            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2655        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2656            H = np.array(refl[:3])
2657            P = np.array(calcControls[phfx+'MustrainAxis'])
2658            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2659            Si = parmDict[phfx+'Mustrain;i']
2660            Sa = parmDict[phfx+'Mustrain;a']
2661            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2662        else:       #generalized - P.W. Stephens model
2663            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2664            Sum = 0
2665            for i,strm in enumerate(Strms):
2666                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2667            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2668    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2669        #crystallite size
2670        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2671            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2672        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2673            H = np.array(refl[:3])
2674            P = np.array(calcControls[phfx+'SizeAxis'])
2675            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2676            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2677            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2678        else:           #ellipsoidal crystallites   #OK
2679            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2680            H = np.array(refl[:3])
2681            lenR = G2pwd.ellipseSize(H,Sij,GB)
2682            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2683        #microstrain               
2684        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2685            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2686        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2687            H = np.array(refl[:3])
2688            P = np.array(calcControls[phfx+'MustrainAxis'])
2689            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2690            Si = parmDict[phfx+'Mustrain;i']
2691            Sa = parmDict[phfx+'Mustrain;a']
2692            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2693        else:       #generalized - P.W. Stephens model  OK
2694            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2695            Sum = 0
2696            for i,strm in enumerate(Strms):
2697                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2698            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2699           
2700    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2701    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2702    sig /= ateln2
2703    return sig,gam
2704       
2705def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2706    'Needs a doc string'
2707    gamDict = {}
2708    sigDict = {}
2709    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:         #All checked & OK
2710        costh = cosd(refl[5+im]/2.)
2711        tanth = tand(refl[5+im]/2.)
2712        #crystallite size derivatives
2713        if calcControls[phfx+'SizeType'] == 'isotropic':
2714            Sgam = 1.8*wave/(np.pi*costh*parmDict[phfx+'Size;i'])
2715            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh*parmDict[phfx+'Size;i']**2)
2716            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2717        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2718            H = np.array(refl[:3])
2719            P = np.array(calcControls[phfx+'SizeAxis'])
2720            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2721            Si = parmDict[phfx+'Size;i']
2722            Sa = parmDict[phfx+'Size;a']
2723            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2724            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2725            Sgam = gami*sqtrm
2726            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2727            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2728            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2729            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2730            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2731            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2732        else:           #ellipsoidal crystallites
2733            const = 1.8*wave/(np.pi*costh)
2734            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2735            H = np.array(refl[:3])
2736            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2737            Sgam = const/lenR
2738            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2739                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2740                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2741        gamDict[phfx+'Size;mx'] = Sgam
2742        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2743               
2744        #microstrain derivatives               
2745        if calcControls[phfx+'MustrainType'] == 'isotropic':
2746            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2747            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2748            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2749        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2750            H = np.array(refl[:3])
2751            P = np.array(calcControls[phfx+'MustrainAxis'])
2752            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2753            Si = parmDict[phfx+'Mustrain;i']
2754            Sa = parmDict[phfx+'Mustrain;a']
2755            gami = 0.018*Si*Sa*tanth/np.pi
2756            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2757            Mgam = gami/sqtrm
2758            dsi = -gami*Si*cosP**2/sqtrm**3
2759            dsa = -gami*Sa*sinP**2/sqtrm**3
2760            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2761            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2762            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2763            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2764        else:       #generalized - P.W. Stephens model
2765            const = 0.018*refl[4+im]**2*tanth/np.pi
2766            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2767            Sum = 0
2768            for i,strm in enumerate(Strms):
2769                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2770                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2771                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2772            Mgam = const*np.sqrt(Sum)
2773            for i in range(len(Strms)):
2774                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2775                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2776        gamDict[phfx+'Mustrain;mx'] = Mgam
2777        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2778    else:   #'T'OF - All checked & OK
2779        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2780            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2781            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2782            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2783        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2784            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2785            H = np.array(refl[:3])
2786            P = np.array(calcControls[phfx+'SizeAxis'])
2787            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2788            Si = parmDict[phfx+'Size;i']
2789            Sa = parmDict[phfx+'Size;a']
2790            gami = const/(Si*Sa)
2791            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2792            Sgam = gami*sqtrm
2793            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2794            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2795            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2796            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2797            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2798            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2799        else:           #OK  ellipsoidal crystallites
2800            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2801            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2802            H = np.array(refl[:3])
2803            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2804            Sgam = const/lenR
2805            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2806                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2807                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2808        gamDict[phfx+'Size;mx'] = Sgam  #OK
2809        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2810               
2811        #microstrain derivatives               
2812        if calcControls[phfx+'MustrainType'] == 'isotropic':
2813            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2814            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2815            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2816        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2817            H = np.array(refl[:3])
2818            P = np.array(calcControls[phfx+'MustrainAxis'])
2819            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2820            Si = parmDict[phfx+'Mustrain;i']
2821            Sa = parmDict[phfx+'Mustrain;a']
2822            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2823            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2824            Mgam = gami/sqtrm
2825            dsi = -gami*Si*cosP**2/sqtrm**3
2826            dsa = -gami*Sa*sinP**2/sqtrm**3
2827            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2828            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2829            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2830            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2831        else:       #generalized - P.W. Stephens model OK
2832            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2833            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2834            Sum = 0
2835            for i,strm in enumerate(Strms):
2836                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2837                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2838                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2839            Mgam = const*np.sqrt(Sum)
2840            for i in range(len(Strms)):
2841                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2842                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2843        gamDict[phfx+'Mustrain;mx'] = Mgam
2844        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2845       
2846    return sigDict,gamDict
2847       
2848def GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict):
2849    'Needs a doc string'
2850    if im:
2851        h,k,l,m = refl[:4]
2852        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2853        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2854    else:
2855        h,k,l = refl[:3]
2856        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2857    refl[4+im] = d
2858    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2859        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2860        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2861        if 'Bragg' in calcControls[hfx+'instType']:
2862            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2863                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2864        else:               #Debye-Scherrer - simple but maybe not right
2865            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+(parmDict[hfx+'DisplaceY']+parmDict[phfx+'LayerDisp'])*sind(pos))
2866    elif 'T' in calcControls[hfx+'histType']:
2867        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2868        #do I need sample position effects - maybe?
2869    return pos
2870
2871def GetReflPosDerv(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict):
2872    'Needs a doc string'
2873    dpr = 180./np.pi
2874    if im:
2875        h,k,l,m = refl[:4]
2876        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2877        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2878        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
2879    else:
2880        m = 0
2881        h,k,l = refl[:3]       
2882        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2883    dst = np.sqrt(dstsq)
2884    dsp = 1./dst
2885    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2886        pos = refl[5+im]-parmDict[hfx+'Zero']
2887        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2888        dpdw = const*dst
2889        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2890        dpdZ = 1.0
2891        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2892            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2893        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2894        if 'Bragg' in calcControls[hfx+'instType']:
2895            dpdSh = -4.*shft*cosd(pos/2.0)
2896            dpdTr = -shft*sind(pos)*100.0
2897            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2898        else:               #Debye-Scherrer - simple but maybe not right
2899            dpdXd = -shft*cosd(pos)
2900            dpdYd = -shft*sind(pos)
2901            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2902    elif 'T' in calcControls[hfx+'histType']:
2903        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2904        dpdZ = 1.0
2905        dpdDC = dsp
2906        dpdDA = dsp**2
2907        dpdDB = 1./dsp
2908        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2909            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2910        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2911           
2912def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2913    'Needs a doc string'
2914    laue = SGData['SGLaue']
2915    uniq = SGData['SGUniq']
2916    h,k,l = refl[:3]
2917    if laue in ['m3','m3m']:
2918        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2919            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2920    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2921        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2922    elif laue in ['3R','3mR']:
2923        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2924    elif laue in ['4/m','4/mmm']:
2925        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2926    elif laue in ['mmm']:
2927        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2928    elif laue in ['2/m']:
2929        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2930        if uniq == 'a':
2931            Dij += parmDict[phfx+'D23']*k*l
2932        elif uniq == 'b':
2933            Dij += parmDict[phfx+'D13']*h*l
2934        elif uniq == 'c':
2935            Dij += parmDict[phfx+'D12']*h*k
2936    else:
2937        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2938            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2939    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2940        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2941    else:
2942        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2943           
2944def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2945    'Needs a doc string'
2946    laue = SGData['SGLaue']
2947    uniq = SGData['SGUniq']
2948    h,k,l = refl[:3]
2949    if laue in ['m3','m3m']:
2950        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2951            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2952    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2953        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2954    elif laue in ['3R','3mR']:
2955        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2956    elif laue in ['4/m','4/mmm']:
2957        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2958    elif laue in ['mmm']:
2959        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2960    elif laue in ['2/m']:
2961        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2962        if uniq == 'a':
2963            dDijDict[phfx+'D23'] = k*l
2964        elif uniq == 'b':
2965            dDijDict[phfx+'D13'] = h*l
2966        elif uniq == 'c':
2967            dDijDict[phfx+'D12'] = h*k
2968    else:
2969        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2970            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2971    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
2972        for item in dDijDict:
2973            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2974    else:
2975        for item in dDijDict:
2976            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2977    return dDijDict
2978   
2979def GetDij(phfx,SGData,parmDict):
2980    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2981    return G2spc.HStrainVals(HSvals,SGData)
2982               
2983def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2984    '''Compute the observed structure factors for Powder histograms and store in reflection array
2985    Multiprocessing support added
2986    '''
2987    if GSASIIpath.GetConfigValue('Show_timing',False):
2988        starttime = time.time() #; print 'start GetFobsSq'
2989    histoList = list(Histograms.keys())
2990    histoList.sort()
2991    Ka2 = shl = lamRatio = kRatio = None
2992    for histogram in histoList:
2993        if 'PWDR' in histogram[:4]:
2994            Histogram = Histograms[histogram]
2995            hId = Histogram['hId']
2996            hfx = ':%d:'%(hId)
2997            Limits = calcControls[hfx+'Limits']
2998            if 'C' in calcControls[hfx+'histType']:
2999                shl = max(parmDict[hfx+'SH/L'],0.0005)
3000                Ka2 = False
3001                kRatio = 0.0
3002                if hfx+'Lam1' in list(parmDict.keys()):
3003                    Ka2 = True
3004                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3005                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
3006            x,y,w,yc,yb,yd = Histogram['Data']
3007            xMask = ma.getmaskarray(x)
3008            xB = np.searchsorted(x,Limits[0])
3009            xF = np.searchsorted(x,Limits[1])
3010            ymb = np.array(y-yb)
3011            ymb = np.where(ymb,ymb,1.0)
3012            ycmb = np.array(yc-yb)
3013            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3014            refLists = Histogram['Reflection Lists']
3015            for phase in refLists:
3016                if phase not in Phases:     #skips deleted or renamed phases silently!
3017                    continue
3018                Phase = Phases[phase]
3019                if histogram not in Phase['Histograms']:
3020                    continue
3021                im = 0
3022                if Phase['General'].get('Modulated',False):
3023                    im = 1
3024                pId = Phase['pId']
3025                phfx = '%d:%d:'%(pId,hId)
3026                refDict = refLists[phase]
3027                sumFo = 0.0
3028                sumdF = 0.0
3029                sumFosq = 0.0
3030                sumdFsq = 0.0
3031                sumInt = 0.0
3032                nExcl = 0
3033                # test to see if we are using multiprocessing below
3034                useMP,ncores = G2mp.InitMP()
3035                if len(refDict['RefList']) < 100: useMP = False       
3036                if useMP: # multiprocessing: create a set of initialized Python processes
3037                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3038                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3039                    profArgs = [[] for i in range(G2mp.ncores)]
3040                else:
3041                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3042                if 'C' in calcControls[hfx+'histType']:
3043                    # are we multiprocessing?
3044                    for iref,refl in enumerate(refDict['RefList']):
3045                        if useMP: 
3046                            profArgs[iref%G2mp.ncores].append((refl,iref))
3047                        else:
3048                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3049                            if type(icod) is tuple:
3050                                refl[8+im] = icod[0]
3051                                sumInt += icod[1]
3052                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3053                            elif icod == -1:
3054                                refl[3+im] *= -1
3055                                nExcl += 1
3056                            elif icod == -2:
3057                                break
3058                    if useMP:
3059                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3060                            sumInt += sInt
3061                            for refl8im,irefl in resList:
3062                                if refl8im is None:
3063                                    refDict['RefList'][irefl][3+im] *= -1
3064                                    nExcl += 1
3065                                else:
3066                                    refDict['RefList'][irefl][8+im] = refl8im
3067                                    if parmDict[phfx+'LeBail']:
3068                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3069                elif 'T' in calcControls[hfx+'histType']:
3070                    for iref,refl in enumerate(refDict['RefList']):
3071                        if useMP: 
3072                            profArgs[iref%G2mp.ncores].append((refl,iref))
3073                        else:
3074                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3075                            if type(icod) is tuple:
3076                                refl[8+im] = icod[0]
3077                                sumInt += icod[1]
3078                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3079                            elif icod == -1:
3080                                refl[3+im] *= -1
3081                                nExcl += 1
3082                            elif icod == -2:
3083                                break
3084                    if useMP:
3085                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3086                            sumInt += sInt
3087                            for refl8im,irefl in resList:
3088                                if refl8im is None:
3089                                    refDict['RefList'][irefl][3+im] *= -1
3090                                    nExcl += 1
3091                                else:
3092                                    refDict['RefList'][irefl][8+im] = refl8im
3093                                    if parmDict[phfx+'LeBail']:
3094                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3095                elif 'B' in calcControls[hfx+'histType']:
3096                    for iref,refl in enumerate(refDict['RefList']):
3097                        if useMP: 
3098                            profArgs[iref%G2mp.ncores].append((refl,iref))
3099                        else:
3100                            icod= G2mp.ComputeFobsSqPink(refl,iref)
3101                            if type(icod) is tuple:
3102                                refl[8+im] = icod[0]
3103                                sumInt += icod[1]
3104                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3105                            elif icod == -1:
3106                                refl[3+im] *= -1
3107                                nExcl += 1
3108                            elif icod == -2:
3109                                break
3110                    if useMP:
3111                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqPinkbatch,profArgs):
3112                            sumInt += sInt
3113                            for refl8im,irefl in resList:
3114                                if refl8im is None:
3115                                    refDict['RefList'][irefl][3+im] *= -1
3116                                    nExcl += 1
3117                                else:
3118                                    refDict['RefList'][irefl][8+im] = refl8im
3119                                    if parmDict[phfx+'LeBail']:
3120                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3121                if useMP: MPpool.terminate()
3122                sumFo = 0.0
3123                sumdF = 0.0
3124                sumFosq = 0.0
3125                sumdFsq = 0.0
3126                for iref,refl in enumerate(refDict['RefList']):
3127                    Fo = np.sqrt(np.abs(refl[8+im]))
3128                    Fc = np.sqrt(np.abs(refl[9]+im))
3129                    sumFo += Fo
3130                    sumFosq += refl[8+im]**2
3131                    sumdF += np.abs(Fo-Fc)
3132                    sumdFsq += (refl[8+im]-refl[9+im])**2
3133                if sumFo:
3134                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3135                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3136                else:
3137                    Histogram['Residuals'][phfx+'Rf'] = 100.
3138                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3139                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3140                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3141                Histogram['Residuals']['hId'] = hId
3142        elif 'HKLF' in histogram[:4]:
3143            Histogram = Histograms[histogram]
3144            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3145    if GSASIIpath.GetConfigValue('Show_timing',False):
3146        print ('GetFobsSq t=',time.time()-starttime)
3147               
3148def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup,histogram=None):
3149    'Computes the powder pattern for a histogram based on contributions from all used phases'
3150    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3151   
3152    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3153        U = parmDict[hfx+'U']
3154        V = parmDict[hfx+'V']
3155        W = parmDict[hfx+'W']
3156        X = parmDict[hfx+'X']
3157        Y = parmDict[hfx+'Y']
3158        Z = parmDict[hfx+'Z']
3159        tanPos = tand(refl[5+im]/2.0)
3160        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3161        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3162        sig = max(0.001,sig)
3163        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3164        gam = max(0.001,gam)
3165        return sig,gam
3166               
3167    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3168        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3169            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3170        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3171        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3172        sig += Ssig
3173        gam += Sgam
3174        return sig,gam
3175       
3176    def GetReflAlpBet(refl,im,hfx,parmDict):
3177        alp = parmDict[hfx+'alpha']/refl[4+im]
3178        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3179        return alp,bet
3180       
3181    def GetPinkReflAlpBet(refl,im,hfx,parmDict):
3182        tanPos = tand(refl[5+im]/2.0)
3183        alp = parmDict[hfx+'alpha-0']+parmDict[hfx+'alpha-0']*tanPos
3184        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']*tanPos
3185        return alp,bet
3186
3187    hId = Histogram['hId']
3188    hfx = ':%d:'%(hId)
3189    bakType = calcControls[hfx+'bakType']
3190    fixedBkg = {i:Histogram['Background'][1][i] for i in Histogram['Background'][1] if i.startswith("_")} 
3191    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixedBkg)
3192    yc = np.zeros_like(yb)
3193    cw = np.diff(ma.getdata(x))
3194    cw = np.append(cw,cw[-1])
3195       
3196    if 'C' in calcControls[hfx+'histType']:   
3197        shl = max(parmDict[hfx+'SH/L'],0.002)
3198        Ka2 = False
3199        if hfx+'Lam1' in (parmDict.keys()):
3200            wave = parmDict[hfx+'Lam1']
3201            Ka2 = True
3202            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3203            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3204        else:
3205            wave = parmDict[hfx+'Lam']
3206    elif 'B' in calcControls[hfx+'histType']:
3207            wave = parmDict[hfx+'Lam']
3208    else:
3209        shl = 0.
3210    for phase in Histogram['Reflection Lists']:
3211        refDict = Histogram['Reflection Lists'][phase]
3212        if phase not in Phases:     #skips deleted or renamed phases silently!
3213            continue
3214        Phase = Phases[phase]
3215        if histogram and not histogram in Phase['Histograms']:
3216            continue
3217        pId = Phase['pId']
3218        pfx = '%d::'%(pId)
3219        phfx = '%d:%d:'%(pId,hId)
3220        hfx = ':%d:'%(hId)
3221        SGData = Phase['General']['SGData']
3222        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3223        im = 0
3224        if Phase['General'].get('Modulated',False):
3225            SSGData = Phase['General']['SSGData']
3226            im = 1  #offset in SS reflection list
3227            #??
3228        Dij = GetDij(phfx,SGData,parmDict)
3229        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do someting if Dij << 0.
3230        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3231        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3232            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3233        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3234        Vst = np.sqrt(nl.det(G))    #V*
3235        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3236            if im:
3237                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3238            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3239                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3240            else:
3241                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3242        badPeak = False
3243        # test to see if we are using multiprocessing here
3244        useMP,ncores = G2mp.InitMP()
3245        if len(refDict['RefList']) < 100: useMP = False       
3246        if useMP: # multiprocessing: create a set of initialized Python processes
3247            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3248            profArgs = [[] for i in range(ncores)]
3249        if 'C' in calcControls[hfx+'histType']:
3250            for iref,refl in enumerate(refDict['RefList']):
3251                if im:
3252                    h,k,l,m = refl[:4]
3253                else:
3254                    h,k,l = refl[:3]
3255                Uniq = np.inner(refl[:3],SGMT)
3256                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
3257                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3258                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3259                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3260                refl[11+im] *= Vst*Lorenz
3261                 
3262                if Phase['General'].get('doPawley'):
3263                    try:
3264                        if im:
3265                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3266                        else:
3267                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3268                        refl[9+im] = parmDict[pInd]
3269                    except KeyError:
3270#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3271                        continue
3272                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3273                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3274                iFin = np.searchsorted(x,refl[5+im]+fmax)
3275                if not iBeg+iFin:       #peak below low limit - skip peak
3276                    continue
3277                elif not iBeg-iFin:     #peak above high limit - done
3278                    break
3279                elif iBeg > iFin:   #bad peak coeff - skip
3280                    badPeak = True
3281                    continue
3282                if useMP:
3283                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3284                else:
3285                    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
3286                if Ka2:
3287                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3288                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3289                    iBeg = np.searchsorted(x,pos2-fmin)
3290                    iFin = np.searchsorted(x,pos2+fmax)
3291                    if not iBeg+iFin:       #peak below low limit - skip peak
3292                        continue
3293                    elif not iBeg-iFin:     #peak above high limit - done
3294                        return yc,yb
3295                    elif iBeg > iFin:   #bad peak coeff - skip
3296                        continue
3297                    if useMP:
3298                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3299                    else:
3300                        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
3301        elif 'B' in calcControls[hfx+'histType']:
3302            for iref,refl in enumerate(refDict['RefList']):
3303                if im:
3304                    h,k,l,m = refl[:4]
3305                else:
3306                    h,k,l = refl[:3]
3307                Uniq = np.inner(refl[:3],SGMT)
3308                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
3309                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3310                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3311                refl[12+im:14+im] = GetPinkReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3312                refl[11+im],refl[14+im],refl[15+im],refl[16+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3313                refl[11+im] *= Vst*Lorenz
3314                 
3315                if Phase['General'].get('doPawley'):
3316                    try:
3317                        if im:
3318                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3319                        else:
3320                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3321                        refl[9+im] = parmDict[pInd]
3322                    except KeyError:
3323#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3324                        continue
3325                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.)
3326                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3327                iFin = np.searchsorted(x,refl[5+im]+fmax)
3