source: trunk/GSASIIstrMath.py @ 3268

Last change on this file since 3268 was 3268, checked in by vondreele, 4 years ago

begin work on incommensurate mag SF calcs.
fix problem with fade calc. for incommensurate site fractions

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