source: trunk/GSASIIstrMath.py @ 3297

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

progress on magnetism - BNS system ok for monoclinic, orthorhombic, tetragonal, hexagonal & trigonal
a few exceptions remain (cubic, rhombohedral, etc.)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 240.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2018-02-27 22:25:17 +0000 (Tue, 27 Feb 2018) $
8# $Author: vondreele $
9# $Revision: 3297 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 3297 2018-02-27 22:25:17Z 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: 3297 $")
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        Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
692        Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
693        if SGData['SGInv'] and not SGData['SGFixed']:
694            Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
695        Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
696        Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
697        Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
698        if SGData['SGInv'] and not SGData['SGFixed']:
699            Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
700    if 'NC' in calcControls[hfx+'histType']:
701        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
702    elif 'X' in calcControls[hfx+'histType']:
703        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
704        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
705    Uij = np.array(G2lat.U6toUij(Uijdata))
706    bij = Mast*Uij.T
707    blkSize = 100       #no. of reflections in a block - size seems optimal
708    nRef = refDict['RefList'].shape[0]
709    SQ = 1./(2.*refDict['RefList'].T[4])**2
710    if 'N' in calcControls[hfx+'histType']:
711        dat = G2el.getBLvalues(BLtables)
712        refDict['FF']['El'] = list(dat.keys())
713        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
714        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
715        for iel,El in enumerate(refDict['FF']['El']):
716            if El in MFtables:
717                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
718    else:       #'X'
719        dat = G2el.getFFvalues(FFtables,0.)
720        refDict['FF']['El'] = list(dat.keys())
721        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
722        for iel,El in enumerate(refDict['FF']['El']):
723            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
724#reflection processing begins here - big arrays!
725    iBeg = 0
726    while iBeg < nRef:
727        iFin = min(iBeg+blkSize,nRef)
728        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
729        H = refl.T[:3]                          #array(blkSize,3)
730        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
731        TwMask = np.any(H,axis=-1)
732        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
733            for ir in range(blkSize):
734                iref = ir+iBeg
735                if iref in TwDict:
736                    for i in TwDict[iref]:
737                        for n in range(NTL):
738                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
739            TwMask = np.any(H,axis=-1)
740        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
741        SQfactor = 4.0*SQ*twopisq               #ditto prev.
742        if 'T' in calcControls[hfx+'histType']:
743            if 'P' in calcControls[hfx+'histType']:
744                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
745            else:
746                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
747            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
748            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
749        Uniq = np.inner(H,SGMT)
750        Phi = np.inner(H,SGT)
751        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
752        sinp = np.sin(phase)
753        cosp = np.cos(phase)
754        biso = -SQfactor*Uisodata[:,nxs]
755        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
756        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
757        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
758        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
759        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
760        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
761        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: math here??
762            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
763            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
764            if SGData['SGInv'] and not SGData['SGFixed']:
765                mphase = np.hstack((phase,-phase))
766            else:
767                mphase = phase
768            mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
769            mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
770            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
771            cosm = np.cos(mphase)                               #ditto
772            HM = np.inner(Bmat.T,H)                             #put into cartesian space
773            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
774            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
775            Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
776            fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
777            fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
778            fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
779            fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
780            refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
781            refl.T[7] = np.copy(refl.T[9])               
782            refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
783        else:
784            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
785            if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
786                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
787                fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
788            else:
789                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
790                fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
791            fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
792            fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
793            if SGData['SGInv']: #centrosymmetric; B=0
794                fbs[0] *= 0.
795                fas[1] *= 0.
796            if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
797                refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
798                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
799                if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
800                    refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
801            else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
802                if len(TwinLaw) > 1:
803                    refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
804                    refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
805                        np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
806                    refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
807                else:   # checked correct!!
808                    refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
809                    refl.T[7] = np.copy(refl.T[9])               
810                    refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
811#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
812        iBeg += blkSize
813#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
814   
815def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
816    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
817    faster than StructureFactorDerv - correct for powders/nontwins!!
818    input:
819   
820    :param dict refDict: where
821        'RefList' list where each ref = h,k,l,it,d,...
822        'FF' dict of form factors - filled in below
823    :param np.array G:      reciprocal metric tensor
824    :param str hfx:    histogram id string
825    :param str pfx:    phase id string
826    :param dict SGData: space group info. dictionary output from SpcGroup
827    :param dict calcControls:
828    :param dict parmDict:
829   
830    :returns: dict dFdvDict: dictionary of derivatives
831    '''
832    phfx = pfx.split(':')[0]+hfx
833    ast = np.sqrt(np.diag(G))
834    Mast = twopisq*np.multiply.outer(ast,ast)
835    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
836    SGT = np.array([ops[1] for ops in SGData['SGOps']])
837    FFtables = calcControls['FFtables']
838    BLtables = calcControls['BLtables']
839    Amat,Bmat = G2lat.Gmat2AB(G)
840    nRef = len(refDict['RefList'])
841    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
842        GetAtomFXU(pfx,calcControls,parmDict)
843    if not Xdata.size:          #no atoms in phase!
844        return {}
845    mSize = len(Mdata)
846    FF = np.zeros(len(Tdata))
847    if 'NC' in calcControls[hfx+'histType']:
848        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
849    elif 'X' in calcControls[hfx+'histType']:
850        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
851        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
852    Uij = np.array(G2lat.U6toUij(Uijdata))
853    bij = Mast*Uij.T
854    dFdvDict = {}
855    dFdfr = np.zeros((nRef,mSize))
856    dFdx = np.zeros((nRef,mSize,3))
857    dFdui = np.zeros((nRef,mSize))
858    dFdua = np.zeros((nRef,mSize,6))
859    dFdbab = np.zeros((nRef,2))
860    dFdfl = np.zeros((nRef))
861    Flack = 1.0
862    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
863        Flack = 1.-2.*parmDict[phfx+'Flack']
864    time0 = time.time()
865#reflection processing begins here - big arrays!
866    iBeg = 0
867    blkSize = 32       #no. of reflections in a block - optimized for speed
868    while iBeg < nRef:
869        iFin = min(iBeg+blkSize,nRef)
870        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
871        H = refl.T[:3].T
872        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
873        SQfactor = 8.0*SQ*np.pi**2
874        if 'T' in calcControls[hfx+'histType']:
875            if 'P' in calcControls[hfx+'histType']:
876                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
877            else:
878                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
879            FP = np.repeat(FP.T,len(SGT),axis=0)
880            FPP = np.repeat(FPP.T,len(SGT),axis=0)
881        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
882        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
883        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
884        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
885        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
886        Phi = np.inner(H,SGT)
887        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
888        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
889        cosp = np.cos(phase)
890        occ = Mdata*Fdata/len(SGT)
891        biso = -SQfactor*Uisodata[:,nxs]
892        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
893        HbH = np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
894        Tuij = np.where(HbH<1.,np.exp(-HbH),1.0).T
895        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
896        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
897        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
898        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
899        if len(FPP.shape) > 1:
900            fotp = np.reshape(FPP,cosp.shape)*Tcorr
901        else:
902            fotp = FPP*Tcorr     
903#            GSASIIpath.IPyBreak()
904        if 'T' in calcControls[hfx+'histType']:
905            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
906            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
907        else:
908            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
909            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
910        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
911        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
912        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,nEqv,nAtoms)
913        fbx = np.array([fot*cosp,-fotp*sinp])
914        #sum below is over Uniq
915        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
916        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
917        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
918        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
919        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
920        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
921        if not SGData['SGInv']:
922            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
923            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
924            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
925            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
926            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
927            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
928            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
929        else:
930            dfbdfr = np.zeros_like(dfadfr)
931            dfbdx = np.zeros_like(dfadx)
932            dfbdui = np.zeros_like(dfadui)
933            dfbdua = np.zeros_like(dfadua)
934            dfbdba = np.zeros_like(dfadba)
935            dfadfl = 0.0
936            dfbdfl = 0.0
937        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
938        SA = fas[0]+fas[1]
939        SB = fbs[0]+fbs[1]
940        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
941            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
942            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
943            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
944            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
945        else:
946            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
947            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
948            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
949            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
950            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
951        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)])+ \
952                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
953#        GSASIIpath.IPyBreak()
954        iBeg += blkSize
955#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
956        #loop over atoms - each dict entry is list of derivatives for all the reflections
957    for i in range(len(Mdata)):
958        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
959        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
960        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
961        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
962        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
963        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
964        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
965        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
966        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
967        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
968        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
969    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
970    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
971    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
972    return dFdvDict
973   
974def StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
975    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
976    input:
977   
978    :param dict refDict: where
979        'RefList' list where each ref = h,k,l,it,d,...
980        'FF' dict of form factors - filled in below
981    :param np.array G:      reciprocal metric tensor
982    :param str hfx:    histogram id string
983    :param str pfx:    phase id string
984    :param dict SGData: space group info. dictionary output from SpcGroup
985    :param dict calcControls:
986    :param dict parmDict:
987   
988    :returns: dict dFdvDict: dictionary of derivatives
989    '''
990    #TODO: fix mag math - moments parallel to crystal axes
991    ast = np.sqrt(np.diag(G))
992    Mast = twopisq*np.multiply.outer(ast,ast)
993    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
994    SGT = np.array([ops[1] for ops in SGData['SGOps']])
995    Ncen = len(SGData['SGCen'])
996    if SGData['SGFixed']:
997        Nops = len(SGMT)*Ncen
998    else:
999        Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
1000    Amat,Bmat = G2lat.Gmat2AB(G)
1001    nRef = len(refDict['RefList'])
1002    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1003        GetAtomFXU(pfx,calcControls,parmDict)
1004    if not Xdata.size:          #no atoms in phase!
1005        return {}
1006    mSize = len(Mdata)
1007    Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
1008#    Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
1009    dGdM = np.repeat(Gdata[:,nxs,:],Nops,axis=1)
1010#    Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
1011    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
1012    if SGData['SGInv'] and not SGData['SGFixed']:
1013        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1014    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
1015    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
1016#    Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps
1017#    Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last - Mxyz,Nops,Natms
1018    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1019    dGdm = (1.-Gdata**2)                        #1/Mag removed - canceled out in dqmx=sum(dqdm*dGdm)
1020    dFdMx = np.zeros((nRef,mSize,3))
1021    Uij = np.array(G2lat.U6toUij(Uijdata))
1022    bij = Mast*Uij.T
1023    dFdvDict = {}
1024    dFdfr = np.zeros((nRef,mSize))
1025    dFdx = np.zeros((nRef,mSize,3))
1026    dFdMx = np.zeros((3,nRef,mSize))
1027    dFdui = np.zeros((nRef,mSize))
1028    dFdua = np.zeros((nRef,mSize,6))
1029    time0 = time.time()
1030#reflection processing begins here - big arrays!
1031    iBeg = 0
1032    blkSize = 32       #no. of reflections in a block - optimized for speed
1033    while iBeg < nRef:
1034        iFin = min(iBeg+blkSize,nRef)
1035        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1036        H = refl.T[:3].T
1037        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1038        SQfactor = 8.0*SQ*np.pi**2
1039        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
1040        Phi = np.inner(H,SGT)
1041        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1042        occ = Mdata*Fdata/Nops
1043        biso = -SQfactor*Uisodata[:,nxs]
1044        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
1045        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1046        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1047        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1048        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))
1049        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1050        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1051        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1052        if SGData['SGInv'] and not SGData['SGFixed']:
1053            mphase = np.hstack((phase,-phase))
1054            Uniq = np.hstack((Uniq,-Uniq))      #Nref,Nops,hkl
1055            Hij = np.hstack((Hij,Hij))
1056        else:
1057            mphase = phase
1058        Hij = np.concatenate(np.array([Hij for cen in SGData['SGCen']]),axis=1)
1059        Uniq = np.hstack([Uniq for cen in SGData['SGCen']])
1060        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1061        mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
1062        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1063        cosm = np.cos(mphase)                               #ditto
1064        HM = np.inner(Bmat.T,H)                             #put into cartesian space
1065        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit vector for H
1066        eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
1067        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1068        NQ = np.where(np.abs(Q)>0.,1./np.abs(Q),0.)     #this sort of works esp for 1 axis moments
1069#        NQ2 = np.where(np.abs(Q)>0.,1./np.sqrt(np.sum(Q**2,axis=0)),0.)
1070        dqdm = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T   #Mxyz,Mxyz,Nref (3x3 matrix)
1071        dqmx = dqdm[:,:,:,nxs,nxs]*dGdm[:,nxs,nxs,:,:]
1072        dqmx2 = np.sum(dqmx,axis=1)   #matrix * vector = vector
1073#        dqmx1 = np.swapaxes(np.swapaxes(np.inner(dqdm.T,dGdm.T),0,1),2,3)
1074        dmx = NQ*Q*dGdM[:,nxs,:,:]-Q*dqmx2                                #*Mag canceled out of dqmx term
1075       
1076        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
1077        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]
1078        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                      #Mxyz,Nref
1079        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)
1080        famx = -Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*sinm[nxs,:,:,:]   #Mxyz,Nref,Nops,Natom
1081        fbmx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*cosm[nxs,:,:,:]
1082        #sums below are over Nops - real part
1083        dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK
1084        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK
1085        dfadmx = np.sum(dmx*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:],axis=2)
1086        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
1087        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)    #deriv OK? not U12 & U23 in sarc
1088        # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
1089        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
1090        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2)
1091        dfbdmx = np.sum(dmx*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:],axis=2)
1092        dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
1093        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
1094        #accumulate derivatives   
1095        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
1096        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
1097        dFdMx[:,iBeg:iFin,:] = 2.*(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx)                       #problems
1098        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
1099        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #problems U12 & U23 in sarc
1100        iBeg += blkSize
1101    print (' %d derivative time %.4f\r'%(nRef,time.time()-time0))
1102        #loop over atoms - each dict entry is list of derivatives for all the reflections
1103    for i in range(len(Mdata)):
1104        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1105        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1106        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1107        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1108        dFdvDict[pfx+'AMx:'+str(i)] = dFdMx[0,:,i]
1109        dFdvDict[pfx+'AMy:'+str(i)] = dFdMx[1,:,i]
1110        dFdvDict[pfx+'AMz:'+str(i)] = dFdMx[2,:,i]
1111        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1112        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1113        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1114        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1115        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
1116        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
1117        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
1118    return dFdvDict
1119       
1120def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1121    '''Compute structure factor derivatives on blocks of reflections - for twins only
1122    faster than StructureFactorDervTw
1123    input:
1124   
1125    :param dict refDict: where
1126        'RefList' list where each ref = h,k,l,it,d,...
1127        'FF' dict of form factors - filled in below
1128    :param np.array G:      reciprocal metric tensor
1129    :param str hfx:    histogram id string
1130    :param str pfx:    phase id string
1131    :param dict SGData: space group info. dictionary output from SpcGroup
1132    :param dict calcControls:
1133    :param dict parmDict:
1134   
1135    :returns: dict dFdvDict: dictionary of derivatives
1136    '''
1137    phfx = pfx.split(':')[0]+hfx
1138    ast = np.sqrt(np.diag(G))
1139    Mast = twopisq*np.multiply.outer(ast,ast)
1140    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1141    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1142    FFtables = calcControls['FFtables']
1143    BLtables = calcControls['BLtables']
1144    TwDict = refDict.get('TwDict',{})           
1145    NTL = calcControls[phfx+'NTL']
1146    NM = calcControls[phfx+'TwinNMN']+1
1147    TwinLaw = calcControls[phfx+'TwinLaw']
1148    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1149    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1150    nTwin = len(TwinLaw)       
1151    nRef = len(refDict['RefList'])
1152    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1153        GetAtomFXU(pfx,calcControls,parmDict)
1154    if not Xdata.size:          #no atoms in phase!
1155        return {}
1156    mSize = len(Mdata)
1157    FF = np.zeros(len(Tdata))
1158    if 'NC' in calcControls[hfx+'histType']:
1159        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1160    elif 'X' in calcControls[hfx+'histType']:
1161        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1162        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1163    Uij = np.array(G2lat.U6toUij(Uijdata))
1164    bij = Mast*Uij.T
1165    dFdvDict = {}
1166    dFdfr = np.zeros((nRef,nTwin,mSize))
1167    dFdx = np.zeros((nRef,nTwin,mSize,3))
1168    dFdui = np.zeros((nRef,nTwin,mSize))
1169    dFdua = np.zeros((nRef,nTwin,mSize,6))
1170    dFdbab = np.zeros((nRef,nTwin,2))
1171    dFdtw = np.zeros((nRef,nTwin))
1172    time0 = time.time()
1173#reflection processing begins here - big arrays!
1174    iBeg = 0
1175    blkSize = 16       #no. of reflections in a block - optimized for speed
1176    while iBeg < nRef:
1177        iFin = min(iBeg+blkSize,nRef)
1178        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1179        H = refl.T[:3]
1180        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
1181        TwMask = np.any(H,axis=-1)
1182        for ir in range(blkSize):
1183            iref = ir+iBeg
1184            if iref in TwDict:
1185                for i in TwDict[iref]:
1186                    for n in range(NTL):
1187                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1188        TwMask = np.any(H,axis=-1)
1189        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1190        SQfactor = 8.0*SQ*np.pi**2
1191        if 'T' in calcControls[hfx+'histType']:
1192            if 'P' in calcControls[hfx+'histType']:
1193                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1194            else:
1195                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1196            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1197            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1198        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1199        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
1200        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1201        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1202        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
1203        Phi = np.inner(H,SGT)
1204        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1205        sinp = np.sin(phase)
1206        cosp = np.cos(phase)
1207        occ = Mdata*Fdata/len(SGT)
1208        biso = -SQfactor*Uisodata[:,nxs]
1209        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1210        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1211        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1212        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,nTwin,len(SGT),6))
1213        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1214        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
1215        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1216        fotp = FPP*Tcorr       
1217        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1218            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
1219            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
1220        else:
1221            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
1222            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
1223        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1224        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1225        if SGData['SGInv']: #centrosymmetric; B=0
1226            fbs[0] *= 0.
1227            fas[1] *= 0.
1228        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
1229        fbx = np.array([fot*cosp,-fotp*sinp])
1230        #sum below is over Uniq
1231        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
1232        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1233        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
1234        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'
1235        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 
1236        if not SGData['SGInv']:
1237            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
1238            dfadba /= 2.
1239#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1240            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
1241            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
1242            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
1243        else:
1244            dfbdfr = np.zeros_like(dfadfr)
1245            dfbdx = np.zeros_like(dfadx)
1246            dfbdui = np.zeros_like(dfadui)
1247            dfbdua = np.zeros_like(dfadua)
1248#            dfbdba = np.zeros_like(dfadba)
1249        SA = fas[0]+fas[1]
1250        SB = fbs[0]+fbs[1]
1251#        GSASIIpath.IPyBreak()
1252        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1253        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1254        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1255        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1256        if SGData['SGInv']: #centrosymmetric; B=0
1257            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1258        else:               
1259            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1260#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1261#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1262        iBeg += blkSize
1263#        GSASIIpath.IPyBreak()
1264    print (' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0))
1265    #loop over atoms - each dict entry is list of derivatives for all the reflections
1266    for i in range(len(Mdata)):     #these all OK
1267        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1268        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1269        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1270        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1271        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1272        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1273        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1274        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1275        dFdvDict[pfx+'AU12:'+str(i)] = np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1276        dFdvDict[pfx+'AU13:'+str(i)] = np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1277        dFdvDict[pfx+'AU23:'+str(i)] = np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1278    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1279    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1280    for i in range(nTwin):
1281        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1282    return dFdvDict
1283   
1284def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1285    '''
1286    Compute super structure factors for all h,k,l,m for phase - no twins
1287    puts the result, F^2, in each ref[9] in refList
1288    works on blocks of 32 reflections for speed
1289    input:
1290   
1291    :param dict refDict: where
1292        'RefList' list where each ref = h,k,l,m,it,d,...
1293        'FF' dict of form factors - filed in below
1294    :param np.array G:      reciprocal metric tensor
1295    :param str pfx:    phase id string
1296    :param dict SGData: space group info. dictionary output from SpcGroup
1297    :param dict calcControls:
1298    :param dict ParmDict:
1299
1300    '''
1301    phfx = pfx.split(':')[0]+hfx
1302    ast = np.sqrt(np.diag(G))
1303    Mast = twopisq*np.multiply.outer(ast,ast)   
1304    SGInv = SGData['SGInv']
1305    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1306    Ncen = len(SGData['SGCen'])
1307    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
1308    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1309    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1310    FFtables = calcControls['FFtables']
1311    BLtables = calcControls['BLtables']
1312    MFtables = calcControls['MFtables']
1313    Amat,Bmat = G2lat.Gmat2AB(G)
1314    Flack = 1.0
1315    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1316        Flack = 1.-2.*parmDict[phfx+'Flack']
1317    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1318        GetAtomFXU(pfx,calcControls,parmDict)
1319    if not Xdata.size:          #no atoms in phase!
1320        return
1321
1322    if parmDict[pfx+'isMag']:       #TODO: fix the math - mag moments now along crystal axes
1323        Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
1324        Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
1325        Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
1326        if SGData['SGInv'] and not SGData['SGFixed']:
1327            Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1328        Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
1329        Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
1330        Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
1331        if SGData['SGInv'] and not SGData['SGFixed']:
1332            Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
1333
1334
1335    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1336    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1337    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1338    FF = np.zeros(len(Tdata))
1339    if 'NC' in calcControls[hfx+'histType']:
1340        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1341    elif 'X' in calcControls[hfx+'histType']:
1342        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1343        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1344    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1345    bij = Mast*Uij
1346    blkSize = 32       #no. of reflections in a block
1347    nRef = refDict['RefList'].shape[0]
1348    SQ = 1./(2.*refDict['RefList'].T[5])**2
1349    if 'N' in calcControls[hfx+'histType']:
1350        dat = G2el.getBLvalues(BLtables)
1351        refDict['FF']['El'] = list(dat.keys())
1352        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1353        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1354        for iel,El in enumerate(refDict['FF']['El']):
1355            if El in MFtables:
1356                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1357    else:
1358        dat = G2el.getFFvalues(FFtables,0.)
1359        refDict['FF']['El'] = list(dat.keys())
1360        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1361        for iel,El in enumerate(refDict['FF']['El']):
1362            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1363    time0 = time.time()
1364#reflection processing begins here - big arrays!
1365    iBeg = 0
1366    while iBeg < nRef:
1367        iFin = min(iBeg+blkSize,nRef)
1368        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1369        H = refl.T[:4]                          #array(blkSize,4)
1370        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1371        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1372        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1373        Uniq = np.inner(H.T,SSGMT)
1374        UniqP = np.inner(HP.T,SGMT)
1375        Phi = np.inner(H.T,SSGT)
1376        if SGInv:   #if centro - expand HKL sets
1377            Uniq = np.hstack((Uniq,-Uniq))
1378            Phi = np.hstack((Phi,-Phi))
1379            UniqP = np.hstack((UniqP,-UniqP))
1380        if 'T' in calcControls[hfx+'histType']:
1381            if 'P' in calcControls[hfx+'histType']:
1382                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1383            else:
1384                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1385            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1386            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1387        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1388        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1389        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1390        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1391        sinp = np.sin(phase)
1392        cosp = np.cos(phase)
1393        biso = -SQfactor*Uisodata[:,nxs]
1394        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1395        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1396        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1397        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1398
1399        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: math here??
1400            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1401            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1402            if SGData['SGInv'] and not SGData['SGFixed']:
1403                mphase = np.hstack((phase,-phase))
1404            else:
1405                mphase = phase
1406            mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1407            mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
1408            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1409            cosm = np.cos(mphase)                               #ditto
1410            HM = np.inner(Bmat.T,H)                             #put into cartesian space
1411            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
1412            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
1413            Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1414            fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1415            fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1416            fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
1417            fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
1418            refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
1419            refl.T[7] = np.copy(refl.T[9])               
1420            refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
1421
1422
1423        else:
1424            if 'T' in calcControls[hfx+'histType']:
1425                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1426                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1427            else:
1428                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1429                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1430            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1431            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1432            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1433            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1434            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1435            if 'P' in calcControls[hfx+'histType']:
1436                refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1437                refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1438            else:
1439                refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1440                refl.T[8] = np.copy(refl.T[10])               
1441                refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1442        iBeg += blkSize
1443    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1444
1445def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1446    '''
1447    Compute super structure factors for all h,k,l,m for phase - twins only
1448    puts the result, F^2, in each ref[8+im] in refList
1449    works on blocks of 32 reflections for speed
1450    input:
1451   
1452    :param dict refDict: where
1453        'RefList' list where each ref = h,k,l,m,it,d,...
1454        'FF' dict of form factors - filed in below
1455    :param np.array G:      reciprocal metric tensor
1456    :param str pfx:    phase id string
1457    :param dict SGData: space group info. dictionary output from SpcGroup
1458    :param dict calcControls:
1459    :param dict ParmDict:
1460
1461    '''
1462    phfx = pfx.split(':')[0]+hfx
1463    ast = np.sqrt(np.diag(G))
1464    Mast = twopisq*np.multiply.outer(ast,ast)   
1465    SGInv = SGData['SGInv']
1466    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1467    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1468    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1469    FFtables = calcControls['FFtables']
1470    BLtables = calcControls['BLtables']
1471    MFtables = calcControls['MFtables']
1472    Flack = 1.0
1473    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1474        Flack = 1.-2.*parmDict[phfx+'Flack']
1475    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1476    TwDict = refDict.get('TwDict',{})           
1477    if 'S' in calcControls[hfx+'histType']:
1478        NTL = calcControls[phfx+'NTL']
1479        NM = calcControls[phfx+'TwinNMN']+1
1480        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1481        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1482        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1483    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1484        GetAtomFXU(pfx,calcControls,parmDict)
1485    if not Xdata.size:          #no atoms in phase!
1486        return
1487    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1488    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1489    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1490    FF = np.zeros(len(Tdata))
1491    if 'NC' in calcControls[hfx+'histType']:
1492        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1493    elif 'X' in calcControls[hfx+'histType']:
1494        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1495        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1496    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1497    bij = Mast*Uij
1498    blkSize = 32       #no. of reflections in a block
1499    nRef = refDict['RefList'].shape[0]
1500    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1501        SQ = 1./(2.*refDict['RefList'].T[5])**2
1502        if 'N' in calcControls[hfx+'histType']:
1503            dat = G2el.getBLvalues(BLtables)
1504            refDict['FF']['El'] = list(dat.keys())
1505            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1506            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1507            for iel,El in enumerate(refDict['FF']['El']):
1508                if El in MFtables:
1509                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1510        else:
1511            dat = G2el.getFFvalues(FFtables,0.)
1512            refDict['FF']['El'] = list(dat.keys())
1513            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1514            for iel,El in enumerate(refDict['FF']['El']):
1515                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1516    time0 = time.time()
1517#reflection processing begins here - big arrays!
1518    iBeg = 0
1519    while iBeg < nRef:
1520        iFin = min(iBeg+blkSize,nRef)
1521        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1522        H = refl[:,:4]                          #array(blkSize,4)
1523        H3 = refl[:,:3]
1524        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1525        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1526        H3 = np.inner(H3,TwinLaw)       
1527        TwMask = np.any(HP,axis=-1)
1528        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1529            for ir in range(blkSize):
1530                iref = ir+iBeg
1531                if iref in TwDict:
1532                    for i in TwDict[iref]:
1533                        for n in range(NTL):
1534                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1535                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1536            TwMask = np.any(HP,axis=-1)
1537        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1538        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1539        Uniq = np.inner(H,SSGMT)
1540        Uniq3 = np.inner(H3,SGMT)
1541        UniqP = np.inner(HP,SGMT)
1542        Phi = np.inner(H,SSGT)
1543        if SGInv:   #if centro - expand HKL sets
1544            Uniq = np.hstack((Uniq,-Uniq))
1545            Uniq3 = np.hstack((Uniq3,-Uniq3))
1546            Phi = np.hstack((Phi,-Phi))
1547            UniqP = np.hstack((UniqP,-UniqP))
1548        if 'T' in calcControls[hfx+'histType']:
1549            if 'P' in calcControls[hfx+'histType']:
1550                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1551            else:
1552                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1553            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1554            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1555        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1556        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1557        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1558        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1559        sinp = np.sin(phase)
1560        cosp = np.cos(phase)
1561        biso = -SQfactor*Uisodata[:,nxs]
1562        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1563        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1564        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1565        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1566#        GSASIIpath.IPyBreak()
1567        if 'T' in calcControls[hfx+'histType']:
1568            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1569            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1570        else:
1571            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1572            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1573        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1574        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1575        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1576        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1577        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1578        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1579        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1580            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1581        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1582        iBeg += blkSize
1583    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1584
1585def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1586    '''
1587    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1588    input:
1589   
1590    :param dict refDict: where
1591        'RefList' list where each ref = h,k,l,m,it,d,...
1592        'FF' dict of form factors - filled in below
1593    :param int im: = 1 (could be eliminated)
1594    :param np.array G:      reciprocal metric tensor
1595    :param str hfx:    histogram id string
1596    :param str pfx:    phase id string
1597    :param dict SGData: space group info. dictionary output from SpcGroup
1598    :param dict SSGData: super space group info.
1599    :param dict calcControls:
1600    :param dict ParmDict:
1601   
1602    :returns: dict dFdvDict: dictionary of derivatives
1603    '''
1604    phfx = pfx.split(':')[0]+hfx
1605    ast = np.sqrt(np.diag(G))
1606    Mast = twopisq*np.multiply.outer(ast,ast)
1607    SGInv = SGData['SGInv']
1608    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1609    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1610    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1611    FFtables = calcControls['FFtables']
1612    BLtables = calcControls['BLtables']
1613    nRef = len(refDict['RefList'])
1614    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1615        GetAtomFXU(pfx,calcControls,parmDict)
1616    if not Xdata.size:          #no atoms in phase!
1617        return {}
1618    mSize = len(Mdata)  #no. atoms
1619    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1620    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1621    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1622    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1623    FF = np.zeros(len(Tdata))
1624    if 'NC' in calcControls[hfx+'histType']:
1625        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1626    elif 'X' in calcControls[hfx+'histType']:
1627        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1628        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1629    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1630    bij = Mast*Uij
1631    if not len(refDict['FF']):
1632        if 'N' in calcControls[hfx+'histType']:
1633            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1634        else:
1635            dat = G2el.getFFvalues(FFtables,0.)       
1636        refDict['FF']['El'] = list(dat.keys())
1637        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1638    dFdvDict = {}
1639    dFdfr = np.zeros((nRef,mSize))
1640    dFdx = np.zeros((nRef,mSize,3))
1641    dFdui = np.zeros((nRef,mSize))
1642    dFdua = np.zeros((nRef,mSize,6))
1643    dFdbab = np.zeros((nRef,2))
1644    dFdfl = np.zeros((nRef))
1645    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1646    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1647    dFdGz = np.zeros((nRef,mSize,5))
1648    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1649    Flack = 1.0
1650    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1651        Flack = 1.-2.*parmDict[phfx+'Flack']
1652    time0 = time.time()
1653    nRef = len(refDict['RefList'])/100
1654    for iref,refl in enumerate(refDict['RefList']):
1655        if 'T' in calcControls[hfx+'histType']:
1656            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1657        H = np.array(refl[:4])
1658        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1659        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1660        SQfactor = 8.0*SQ*np.pi**2
1661        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1662        Bab = parmDict[phfx+'BabA']*dBabdA
1663        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1664        FF = refDict['FF']['FF'][iref].T[Tindx]
1665        Uniq = np.inner(H,SSGMT)
1666        Phi = np.inner(H,SSGT)
1667        UniqP = np.inner(HP,SGMT)
1668        if SGInv:   #if centro - expand HKL sets
1669            Uniq = np.vstack((Uniq,-Uniq))
1670            Phi = np.hstack((Phi,-Phi))
1671            UniqP = np.vstack((UniqP,-UniqP))
1672        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1673        sinp = np.sin(phase)
1674        cosp = np.cos(phase)
1675        occ = Mdata*Fdata/Uniq.shape[0]
1676        biso = -SQfactor*Uisodata[:,nxs]
1677        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1678        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1679        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1680        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1681        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1682        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1683        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1684        fotp = FPP*Tcorr            #ops x atoms
1685        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1686        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1687        # GfpuA is 2 x ops x atoms
1688        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1689        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1690        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1691        fag = fa*GfpuA[0]-fb*GfpuA[1]
1692        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1693       
1694        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1695        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1696        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1697        fbx = np.array([fot*cosp,-fotp*sinp])
1698        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1699        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1700        #sum below is over Uniq
1701        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1702        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1703        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1704        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1705        dfadui = np.sum(-SQfactor*fag,axis=1)
1706        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1707        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1708        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1709        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1710        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1711        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1712        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1713        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1714        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1715        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1716        dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1717        dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1)
1718        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1719        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1720        if not SGData['SGInv']:   #Flack derivative
1721            dfadfl = np.sum(-FPP*Tcorr*sinp)
1722            dfbdfl = np.sum(FPP*Tcorr*cosp)
1723        else:
1724            dfadfl = 1.0
1725            dfbdfl = 1.0
1726#        GSASIIpath.IPyBreak()
1727        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1728        SA = fas[0]+fas[1]      #float = A+A'
1729        SB = fbs[0]+fbs[1]      #float = B+B'
1730        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1731            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1732            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1733                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1734            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1735                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1736            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1737                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1738            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1739                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1740            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1741                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1742            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1743                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1744            dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
1745                2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
1746            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1747                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1748        else:                       #OK, I think
1749            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1750            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1751            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1752            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1753            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1754                           
1755            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1756            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1757            dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
1758            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1759#            GSASIIpath.IPyBreak()
1760        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1761            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1762        #loop over atoms - each dict entry is list of derivatives for all the reflections
1763        if not iref%100 :
1764            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1765    for i in range(len(Mdata)):     #loop over atoms
1766        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1767        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1768        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1769        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1770        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1771        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1772        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1773        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1774        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
1775        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
1776        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
1777        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1778            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1779            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1780        nx = 0
1781        if waveTypes[i] in ['Block','ZigZag']:
1782            nx = 1 
1783            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
1784            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
1785            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
1786            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
1787            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
1788        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1789            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1790            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1791            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1792            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1793            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1794            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1795        for j in range(USSdata.shape[1]):       #loop over waves
1796            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1797            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1798            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1799            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
1800            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
1801            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
1802            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1803            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1804            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1805            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
1806            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
1807            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
1808           
1809#        GSASIIpath.IPyBreak()
1810    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
1811    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1812    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1813    return dFdvDict
1814
1815def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1816    'Needs a doc string - no twins'
1817    phfx = pfx.split(':')[0]+hfx
1818    ast = np.sqrt(np.diag(G))
1819    Mast = twopisq*np.multiply.outer(ast,ast)
1820    SGInv = SGData['SGInv']
1821    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1822    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1823    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1824    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1825    FFtables = calcControls['FFtables']
1826    BLtables = calcControls['BLtables']
1827    nRef = len(refDict['RefList'])
1828    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1829        GetAtomFXU(pfx,calcControls,parmDict)
1830    if not Xdata.size:          #no atoms in phase!
1831        return {}
1832    mSize = len(Mdata)  #no. atoms
1833    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1834    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1835    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1836    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1837    FF = np.zeros(len(Tdata))
1838    if 'NC' in calcControls[hfx+'histType']:
1839        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1840    elif 'X' in calcControls[hfx+'histType']:
1841        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1842        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1843    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1844    bij = Mast*Uij
1845    if not len(refDict['FF']):
1846        if 'N' in calcControls[hfx+'histType']:
1847            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1848        else:
1849            dat = G2el.getFFvalues(FFtables,0.)       
1850        refDict['FF']['El'] = list(dat.keys())
1851        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1852    dFdvDict = {}
1853    dFdfr = np.zeros((nRef,mSize))
1854    dFdx = np.zeros((nRef,mSize,3))
1855    dFdui = np.zeros((nRef,mSize))
1856    dFdua = np.zeros((nRef,mSize,6))
1857    dFdbab = np.zeros((nRef,2))
1858    dFdfl = np.zeros((nRef))
1859    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1860    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1861    dFdGz = np.zeros((nRef,mSize,5))
1862    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1863    Flack = 1.0
1864    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1865        Flack = 1.-2.*parmDict[phfx+'Flack']
1866    time0 = time.time()
1867    iBeg = 0
1868    blkSize = 4       #no. of reflections in a block - optimized for speed
1869    while iBeg < nRef:
1870        iFin = min(iBeg+blkSize,nRef)
1871        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1872        H = refl.T[:4]
1873        HP = H[:3].T+modQ*H.T[:,3:]            #projected hklm to hkl
1874        SQ = 1./(2.*refl.T[4+im])**2             # or (sin(theta)/lambda)**2
1875        SQfactor = 8.0*SQ*np.pi**2
1876        if 'T' in calcControls[hfx+'histType']:
1877            if 'P' in calcControls[hfx+'histType']:
1878                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[15])
1879            else:
1880                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[13])
1881            FP = np.repeat(FP.T,len(SGT),axis=0)
1882            FPP = np.repeat(FPP.T,len(SGT),axis=0)
1883#        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1884        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
1885        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1886        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
1887        Uniq = np.inner(H.T,SSGMT)
1888        Phi = np.inner(H.T,SSGT)
1889        UniqP = np.inner(HP,SGMT)
1890        if SGInv:   #if centro - expand HKL sets
1891            Uniq = np.hstack((Uniq,-Uniq))
1892            Phi = np.hstack((Phi,-Phi))
1893            UniqP = np.hstack((UniqP,-UniqP))
1894            FF = np.vstack((FF,FF))
1895            Bab = np.concatenate((Bab,Bab))
1896        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata+Xdata).T)+Phi[:,:,nxs])
1897        sinp = np.sin(phase)
1898        cosp = np.cos(phase)
1899        occ = Mdata*Fdata/Uniq.shape[1]
1900        biso = -SQfactor*Uisodata[:,nxs]
1901        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T    #ops x atoms
1902        HbH = -np.sum(UniqP[:,:,nxs,:3]*np.inner(UniqP[:,:,:3],bij),axis=-1)  #ops x atoms
1903        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in np.reshape(UniqP,(-1,3))]) #atoms x 3x3
1904        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(iFin-iBeg,-1,6))                     #atoms x 6
1905        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1906#        GSASIIpath.IPyBreak()
1907        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1908        fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr     #ops x atoms
1909        fotp = FPP*Tcorr            #ops x atoms
1910        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1911        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1912        # GfpuA is 2 x ops x atoms
1913        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1914        fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1915        fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1916        fag = fa*GfpuA[0]-fb*GfpuA[1]
1917        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1918       
1919        fas = np.sum(np.sum(fag,axis=-1),axis=-1)     # 2 x refBlk
1920        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1921        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1922        fbx = np.array([fot*cosp,-fotp*sinp])
1923        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1924        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1925        #sum below is over Uniq
1926        dfadfr = np.sum(fag/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
1927        dfbdfr = np.sum(fbg/occ,axis=-2)        #Fdata != 0 avoids /0. problem
1928#        dfadba = np.sum(-cosp*Tcorr,axis=-2)
1929#        dfbdba = np.sum(-sinp*Tcorr,axis=-2)
1930        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fag,axis=-2)
1931        dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fbg,axis=-2)
1932        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fax[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
1933        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fbx[:,:,:,:,nxs],axis=-3)  #2 refBlk x x nAtom x 3xyz; sum nOps
1934        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fag[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
1935        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbg[:,:,:,:,nxs],axis=-3)             #2 x nAtom x 6Uij; sum nOps
1936        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1937        dfadGf = np.sum(fa[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
1938        dfbdGf = np.sum(fb[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2)
1939        dfadGx = np.sum(fa[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
1940        dfbdGx = np.sum(fb[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2)
1941        dfadGz = np.sum(fa[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
1942        dfbdGz = np.sum(fb[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2)
1943        dfadGu = np.sum(fa[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)
1944        dfbdGu = np.sum(fb[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2)   
1945        if not SGData['SGInv']:   #Flack derivative
1946            dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1)
1947            dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1)
1948        else:
1949            dfadfl = 1.0
1950            dfbdfl = 1.0
1951        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
1952        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
1953        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
1954        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1955            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1956            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(Uniq)+   \
1957                2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(Uniq)
1958            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
1959                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
1960            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
1961                2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
1962            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
1963                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
1964            dFdGf[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGf[0]+fas[1,:,nxs,nxs,nxs]*dfadGf[1])+  \
1965                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGf[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGf[1])
1966            dFdGx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGx[0]+fas[1,:,nxs,nxs,nxs]*dfadGx[1])+  \
1967                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGx[0]-fbs[1,:,nxs,nxs,nxs]*dfbdGx[1])
1968            dFdGz[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadGz[0]+fas[1,:,nxs,nxs]*dfadGz[1])+  \
1969                2.*(fbs[0,:,nxs,nxs]*dfbdGz[0]+fbs[1,:,nxs,nxs]*dfbdGz[1])
1970            dFdGu[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGu[0]+fas[1,:,nxs,nxs,nxs]*dfadGu[1])+  \
1971                2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGu[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGu[1])
1972        else:                       #OK, I think
1973            dFdfr[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadfr[0]+dfadfr[1])+SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) #array(nRef,nAtom)
1974            dFdx[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1]))    #array(nRef,nAtom,3)
1975            dFdui[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadui[0]+dfadui[1])+SB[:,nxs]*(dfbdui[0]+dfbdui[1]))   #array(nRef,nAtom)
1976            dFdua[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1]))    #array(nRef,nAtom,6)
1977            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1978                           
1979            dFdGf[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGf[0]+SB[:,nxs,nxs,nxs]*dfbdGf[1])      #array(nRef,natom,nwave,2)
1980            dFdGx[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGx[0]+SB[:,nxs,nxs,nxs]*dfbdGx[1])      #array(nRef,natom,nwave,6)
1981            dFdGz[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*dfadGz[0]+SB[:,nxs,nxs]*dfbdGz[1])      #array(nRef,natom,5)
1982            dFdGu[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGu[0]+SB[:,nxs,nxs,nxs]*dfbdGu[1])      #array(nRef,natom,nwave,12)
1983#            GSASIIpath.IPyBreak()
1984#        dFdbab[iBeg:iFin] = 2.*fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1985#            2.*fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1986        #loop over atoms - each dict entry is list of derivatives for all the reflections
1987        print (' %d derivative time %.4f\r'%(iBeg,time.time()-time0),end='')
1988        iBeg += blkSize
1989    for i in range(len(Mdata)):     #loop over atoms
1990        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1991        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1992        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1993        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1994        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1995        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1996        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1997        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1998        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
1999        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
2000        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
2001        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2002            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2003            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2004        nx = 0
2005        if waveTypes[i] in ['Block','ZigZag']:
2006            nx = 1 
2007            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2008            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2009            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2010            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2011            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2012        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2013            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2014            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2015            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2016            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2017            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2018            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2019        for j in range(USSdata.shape[1]):       #loop over waves
2020            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2021            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2022            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2023            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
2024            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
2025            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
2026            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2027            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2028            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2029            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
2030            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
2031            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
2032           
2033#        GSASIIpath.IPyBreak()
2034    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2035    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2036    return dFdvDict
2037   
2038def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2039    'Needs a doc string'
2040    phfx = pfx.split(':')[0]+hfx
2041    ast = np.sqrt(np.diag(G))
2042    Mast = twopisq*np.multiply.outer(ast,ast)
2043    SGInv = SGData['SGInv']
2044    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2045    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2046    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2047    FFtables = calcControls['FFtables']
2048    BLtables = calcControls['BLtables']
2049    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2050    TwDict = refDict.get('TwDict',{})           
2051    if 'S' in calcControls[hfx+'histType']:
2052        NTL = calcControls[phfx+'NTL']
2053        NM = calcControls[phfx+'TwinNMN']+1
2054        TwinLaw = calcControls[phfx+'TwinLaw']
2055        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2056    nTwin = len(TwinLaw)       
2057    nRef = len(refDict['RefList'])
2058    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2059        GetAtomFXU(pfx,calcControls,parmDict)
2060    if not Xdata.size:          #no atoms in phase!
2061        return {}
2062    mSize = len(Mdata)  #no. atoms
2063    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2064    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
2065    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
2066    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2067    FF = np.zeros(len(Tdata))
2068    if 'NC' in calcControls[hfx+'histType']:
2069        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2070    elif 'X' in calcControls[hfx+'histType']:
2071        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2072        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2073    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2074    bij = Mast*Uij
2075    if not len(refDict['FF']):
2076        if 'N' in calcControls[hfx+'histType']:
2077            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2078        else:
2079            dat = G2el.getFFvalues(FFtables,0.)       
2080        refDict['FF']['El'] = list(dat.keys())
2081        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2082    dFdvDict = {}
2083    dFdfr = np.zeros((nRef,nTwin,mSize))
2084    dFdx = np.zeros((nRef,nTwin,mSize,3))
2085    dFdui = np.zeros((nRef,nTwin,mSize))
2086    dFdua = np.zeros((nRef,nTwin,mSize,6))
2087    dFdbab = np.zeros((nRef,nTwin,2))
2088    dFdtw = np.zeros((nRef,nTwin))
2089    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2090    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2091    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2092    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2093    Flack = 1.0
2094    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2095        Flack = 1.-2.*parmDict[phfx+'Flack']
2096    time0 = time.time()
2097    nRef = len(refDict['RefList'])/100
2098    for iref,refl in enumerate(refDict['RefList']):
2099        if 'T' in calcControls[hfx+'histType']:
2100            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2101        H = np.array(refl[:4])
2102        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2103        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2104        TwMask = np.any(H,axis=-1)
2105        if TwinLaw.shape[0] > 1 and TwDict:
2106            if iref in TwDict:
2107                for i in TwDict[iref]:
2108                    for n in range(NTL):
2109                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2110            TwMask = np.any(H,axis=-1)
2111        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2112        SQfactor = 8.0*SQ*np.pi**2
2113        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2114        Bab = parmDict[phfx+'BabA']*dBabdA
2115        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2116        FF = refDict['FF']['FF'][iref].T[Tindx]
2117        Uniq = np.inner(H,SSGMT)
2118        Phi = np.inner(H,SSGT)
2119        UniqP = np.inner(HP,SGMT)
2120        if SGInv:   #if centro - expand HKL sets
2121            Uniq = np.vstack((Uniq,-Uniq))
2122            Phi = np.hstack((Phi,-Phi))
2123            UniqP = np.vstack((UniqP,-UniqP))
2124        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2125        sinp = np.sin(phase)
2126        cosp = np.cos(phase)
2127        occ = Mdata*Fdata/Uniq.shape[0]
2128        biso = -SQfactor*Uisodata[:,nxs]
2129        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2130        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2131        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2132        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2133        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2134        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2135        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2136        fotp = FPP*Tcorr            #ops x atoms
2137        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2138        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2139        # GfpuA is 2 x ops x atoms
2140        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2141        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2142        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2143        fag = fa*GfpuA[0]-fb*GfpuA[1]
2144        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2145       
2146        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2147        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2148        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2149        fbx = np.array([fot*cosp,-fotp*sinp])
2150        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2151        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2152        #sum below is over Uniq
2153        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2154        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2155        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2156        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2157        dfadui = np.sum(-SQfactor*fag,axis=1)
2158        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2159        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2160        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2161        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2162        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2163        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2164        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2165        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2166        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2167        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2168        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2169        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2170        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2171        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2172#        GSASIIpath.IPyBreak()
2173        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2174        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2175        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2176        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)]
2177        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)]
2178        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)]
2179        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)]
2180        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2181
2182        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2183        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2184        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2185        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2186#            GSASIIpath.IPyBreak()
2187        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2188            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2189        #loop over atoms - each dict entry is list of derivatives for all the reflections
2190        if not iref%100 :
2191            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2192    for i in range(len(Mdata)):     #loop over atoms
2193        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2194        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2195        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2196        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2197        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2198        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2199        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2200        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2201        dFdvDict[pfx+'AU12:'+str(i)] = dFdua.T[3][i]
2202        dFdvDict[pfx+'AU13:'+str(i)] = dFdua.T[4][i]
2203        dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
2204        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2205            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2206            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2207        nx = 0
2208        if waveTypes[i] in ['Block','ZigZag']:
2209            nx = 1 
2210            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2211            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2212            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2213            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2214            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2215        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2216            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2217            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2218            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2219            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2220            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2221            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2222        for j in range(USSdata.shape[1]):       #loop over waves
2223            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2224            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2225            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2226            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = dFdGu.T[3][j][i]
2227            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = dFdGu.T[4][j][i]
2228            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = dFdGu.T[5][j][i]
2229            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2230            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2231            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2232            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = dFdGu.T[9][j][i]
2233            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = dFdGu.T[10][j][i]
2234            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = dFdGu.T[11][j][i]
2235           
2236#        GSASIIpath.IPyBreak()
2237    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2238    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2239    return dFdvDict
2240   
2241def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2242    ''' Single crystal extinction function; returns extinction & derivative
2243    '''
2244    extCor = 1.0
2245    dervDict = {}
2246    dervCor = 1.0
2247    if calcControls[phfx+'EType'] != 'None':
2248        SQ = 1/(4.*ref[4+im]**2)
2249        if 'C' in parmDict[hfx+'Type']:           
2250            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2251        else:   #'T'
2252            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2253        if 'SXC' in parmDict[hfx+'Type']:
2254            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2255            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2256            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2257            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2258        elif 'SNT' in parmDict[hfx+'Type']:
2259            AV = 1.e7/parmDict[pfx+'Vol']**2
2260            PL = SQ
2261            PLZ = AV*ref[9+im]*ref[12+im]**2
2262        elif 'SNC' in parmDict[hfx+'Type']:
2263            AV = 1.e7/parmDict[pfx+'Vol']**2
2264            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2265            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2266           
2267        if 'Primary' in calcControls[phfx+'EType']:
2268            PLZ *= 1.5
2269        else:
2270            if 'C' in parmDict[hfx+'Type']:
2271                PLZ *= calcControls[phfx+'Tbar']
2272            else: #'T'
2273                PLZ *= ref[13+im]      #t-bar
2274        if 'Primary' in calcControls[phfx+'EType']:
2275            PLZ *= 1.5
2276            PSIG = parmDict[phfx+'Ep']
2277        elif 'I & II' in calcControls[phfx+'EType']:
2278            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2279        elif 'Type II' in calcControls[phfx+'EType']:
2280            PSIG = parmDict[phfx+'Es']
2281        else:       # 'Secondary Type I'
2282            PSIG = parmDict[phfx+'Eg']/PL
2283           
2284        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2285        AL = 0.025+0.285*cos2T
2286        BG = 0.02-0.025*cos2T
2287        BL = 0.15-0.2*(0.75-cos2T)**2
2288        if cos2T < 0.:
2289            BL = -0.45*cos2T
2290        CG = 2.
2291        CL = 2.
2292        PF = PLZ*PSIG
2293       
2294        if 'Gaussian' in calcControls[phfx+'EApprox']:
2295            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2296            extCor = np.sqrt(PF4)
2297            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2298        else:
2299            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2300            extCor = np.sqrt(PF4)
2301            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2302
2303        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2304        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2305            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2306        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2307            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2308        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2309            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2310               
2311    return 1./extCor,dervDict,dervCor
2312   
2313def Dict2Values(parmdict, varylist):
2314    '''Use before call to leastsq to setup list of values for the parameters
2315    in parmdict, as selected by key in varylist'''
2316    return [parmdict[key] for key in varylist] 
2317   
2318def Values2Dict(parmdict, varylist, values):
2319    ''' Use after call to leastsq to update the parameter dictionary with
2320    values corresponding to keys in varylist'''
2321    parmdict.update(zip(varylist,values))
2322   
2323def GetNewCellParms(parmDict,varyList):
2324    'Needs a doc string'
2325    newCellDict = {}
2326    Anames = ['A'+str(i) for i in range(6)]
2327    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2328    for item in varyList:
2329        keys = item.split(':')
2330        if keys[2] in Ddict:
2331            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2332            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2333            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2334    return newCellDict          # is e.g. {'0::D11':A0-D11}
2335   
2336def ApplyXYZshifts(parmDict,varyList):
2337    '''
2338    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2339   
2340    :param dict parmDict: parameter dictionary
2341    :param list varyList: list of variables (not used!)
2342    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2343
2344    '''
2345    newAtomDict = {}
2346    for item in parmDict:
2347        if 'dA' in item:
2348            parm = ''.join(item.split('d'))
2349            parmDict[parm] += parmDict[item]
2350            newAtomDict[item] = [parm,parmDict[parm]]
2351    return newAtomDict
2352   
2353def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2354    'Spherical harmonics texture'
2355    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2356    if 'T' in calcControls[hfx+'histType']:
2357        tth = parmDict[hfx+'2-theta']
2358    else:
2359        tth = refl[5+im]
2360    odfCor = 1.0
2361    H = refl[:3]
2362    cell = G2lat.Gmat2cell(g)
2363    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2364    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2365    phi,beta = G2lat.CrsAng(H,cell,SGData)
2366    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2367    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2368    for item in SHnames:
2369        L,M,N = eval(item.strip('C'))
2370        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2371        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2372        Lnorm = G2lat.Lnorm(L)
2373        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2374    return odfCor
2375   
2376def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2377    'Spherical harmonics texture derivatives'
2378    if 'T' in calcControls[hfx+'histType']:
2379        tth = parmDict[hfx+'2-theta']
2380    else:
2381        tth = refl[5+im]
2382    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2383    odfCor = 1.0
2384    dFdODF = {}
2385    dFdSA = [0,0,0]
2386    H = refl[:3]
2387    cell = G2lat.Gmat2cell(g)
2388    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2389    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2390    phi,beta = G2lat.CrsAng(H,cell,SGData)
2391    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2392    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2393    for item in SHnames:
2394        L,M,N = eval(item.strip('C'))
2395        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2396        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2397        Lnorm = G2lat.Lnorm(L)
2398        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2399        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2400        for i in range(3):
2401            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2402    return odfCor,dFdODF,dFdSA
2403   
2404def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2405    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2406    if 'T' in calcControls[hfx+'histType']:
2407        tth = parmDict[hfx+'2-theta']
2408    else:
2409        tth = refl[5+im]
2410    odfCor = 1.0
2411    H = refl[:3]
2412    cell = G2lat.Gmat2cell(g)
2413    Sangls = [0.,0.,0.]
2414    if 'Bragg' in calcControls[hfx+'instType']:
2415        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2416        IFCoup = True
2417    else:
2418        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2419        IFCoup = False
2420    phi,beta = G2lat.CrsAng(H,cell,SGData)
2421    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2422    SHnames = calcControls[phfx+'SHnames']
2423    for item in SHnames:
2424        L,N = eval(item.strip('C'))
2425        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2426        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2427        Lnorm = G2lat.Lnorm(L)
2428        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2429    return np.squeeze(odfCor)
2430   
2431def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2432    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2433    if 'T' in calcControls[hfx+'histType']:
2434        tth = parmDict[hfx+'2-theta']
2435    else:
2436        tth = refl[5+im]
2437    odfCor = 1.0
2438    dFdODF = {}
2439    H = refl[:3]
2440    cell = G2lat.Gmat2cell(g)
2441    Sangls = [0.,0.,0.]
2442    if 'Bragg' in calcControls[hfx+'instType']:
2443        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2444        IFCoup = True
2445    else:
2446        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2447        IFCoup = False
2448    phi,beta = G2lat.CrsAng(H,cell,SGData)
2449    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2450    SHnames = calcControls[phfx+'SHnames']
2451    for item in SHnames:
2452        L,N = eval(item.strip('C'))
2453        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2454        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2455        Lnorm = G2lat.Lnorm(L)
2456        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2457        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2458    return odfCor,dFdODF
2459   
2460def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2461    'March-Dollase preferred orientation correction'
2462    POcorr = 1.0
2463    MD = parmDict[phfx+'MD']
2464    if MD != 1.0:
2465        MDAxis = calcControls[phfx+'MDAxis']
2466        sumMD = 0
2467        for H in uniq:           
2468            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2469            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2470            sumMD += A**3
2471        POcorr = sumMD/len(uniq)
2472    return POcorr
2473   
2474def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2475    'Needs a doc string'
2476    POcorr = 1.0
2477    POderv = {}
2478    if calcControls[phfx+'poType'] == 'MD':
2479        MD = parmDict[phfx+'MD']
2480        MDAxis = calcControls[phfx+'MDAxis']
2481        sumMD = 0
2482        sumdMD = 0
2483        for H in uniq:           
2484            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2485            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2486            sumMD += A**3
2487            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2488        POcorr = sumMD/len(uniq)
2489        POderv[phfx+'MD'] = sumdMD/len(uniq)
2490    else:   #spherical harmonics
2491        if calcControls[phfx+'SHord']:
2492            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2493    return POcorr,POderv
2494   
2495def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2496    'Needs a doc string'
2497    if 'Debye' in calcControls[hfx+'instType']:
2498        if 'T' in calcControls[hfx+'histType']:
2499            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2500        else:
2501            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2502    else:
2503        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2504   
2505def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2506    'Needs a doc string'
2507    if 'Debye' in calcControls[hfx+'instType']:
2508        if 'T' in calcControls[hfx+'histType']:
2509            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2510        else:
2511            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2512    else:
2513        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2514       
2515def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2516    'Needs a doc string'
2517    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2518    pi2 = np.sqrt(2./np.pi)
2519    if 'T' in calcControls[hfx+'histType']:
2520        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2521        wave = refl[14+im]
2522    else:   #'C'W
2523        sth2 = sind(refl[5+im]/2.)**2
2524        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2525    c2th = 1.-2.0*sth2
2526    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2527    if 'X' in calcControls[hfx+'histType']:
2528        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2529    xfac = flv2*parmDict[phfx+'Extinction']
2530    exb = 1.0
2531    if xfac > -1.:
2532        exb = 1./np.sqrt(1.+xfac)
2533    exl = 1.0
2534    if 0 < xfac <= 1.:
2535        xn = np.array([xfac**(i+1) for i in range(6)])
2536        exl += np.sum(xn*coef)
2537    elif xfac > 1.:
2538        xfac2 = 1./np.sqrt(xfac)
2539        exl = pi2*(1.-0.125/xfac)*xfac2
2540    return exb*sth2+exl*(1.-sth2)
2541   
2542def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2543    'Needs a doc string'
2544    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2545    pi2 = np.sqrt(2./np.pi)
2546    if 'T' in calcControls[hfx+'histType']:
2547        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2548        wave = refl[14+im]
2549    else:   #'C'W
2550        sth2 = sind(refl[5+im]/2.)**2
2551        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2552    c2th = 1.-2.0*sth2
2553    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2554    if 'X' in calcControls[hfx+'histType']:
2555        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2556    xfac = flv2*parmDict[phfx+'Extinction']
2557    dbde = -500.*flv2
2558    if xfac > -1.:
2559        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2560    dlde = 0.
2561    if 0 < xfac <= 1.:
2562        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2563        dlde = np.sum(xn*coef)
2564    elif xfac > 1.:
2565        xfac2 = 1./np.sqrt(xfac)
2566        dlde = flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2567       
2568    return dbde*sth2+dlde*(1.-sth2)
2569   
2570def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2571    'Needs a doc string'    #need powder extinction!
2572    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2573    if 'X' in parmDict[hfx+'Type']:
2574        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2575    POcorr = 1.0
2576    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2577        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2578    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2579        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2580    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2581        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2582    Icorr *= POcorr
2583    AbsCorr = 1.0
2584    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2585    Icorr *= AbsCorr
2586    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2587    Icorr *= ExtCorr
2588    return Icorr,POcorr,AbsCorr,ExtCorr
2589   
2590def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2591    'Needs a doc string'    #need powder extinction derivs!
2592    dIdsh = 1./parmDict[hfx+'Scale']
2593    dIdsp = 1./parmDict[phfx+'Scale']
2594    if 'X' in parmDict[hfx+'Type']:
2595        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2596        dIdPola /= pola
2597    else:       #'N'
2598        dIdPola = 0.0
2599    dFdODF = {}
2600    dFdSA = [0,0,0]
2601    dIdPO = {}
2602    if pfx+'SHorder' in parmDict:
2603        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2604        for iSH in dFdODF:
2605            dFdODF[iSH] /= odfCor
2606        for i in range(3):
2607            dFdSA[i] /= odfCor
2608    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2609        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2610        for iPO in dIdPO:
2611            dIdPO[iPO] /= POcorr
2612    if 'T' in parmDict[hfx+'Type']:
2613        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2614        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2615    else:
2616        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2617        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2618    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2619       
2620def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2621    'Needs a doc string'
2622    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2623        costh = cosd(refl[5+im]/2.)
2624        #crystallite size
2625        if calcControls[phfx+'SizeType'] == 'isotropic':
2626            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2627        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2628            H = np.array(refl[:3])
2629            P = np.array(calcControls[phfx+'SizeAxis'])
2630            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2631            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2632            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2633        else:           #ellipsoidal crystallites
2634            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2635            H = np.array(refl[:3])
2636            lenR = G2pwd.ellipseSize(H,Sij,GB)
2637            Sgam = 1.8*wave/(np.pi*costh*lenR)
2638        #microstrain               
2639        if calcControls[phfx+'MustrainType'] == 'isotropic':
2640            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2641        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2642            H = np.array(refl[:3])
2643            P = np.array(calcControls[phfx+'MustrainAxis'])
2644            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2645            Si = parmDict[phfx+'Mustrain;i']
2646            Sa = parmDict[phfx+'Mustrain;a']
2647            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2648        else:       #generalized - P.W. Stephens model
2649            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2650            Sum = 0
2651            for i,strm in enumerate(Strms):
2652                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2653            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2654    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2655        #crystallite size
2656        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2657            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2658        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2659            H = np.array(refl[:3])
2660            P = np.array(calcControls[phfx+'SizeAxis'])
2661            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2662            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2663            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2664        else:           #ellipsoidal crystallites   #OK
2665            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2666            H = np.array(refl[:3])
2667            lenR = G2pwd.ellipseSize(H,Sij,GB)
2668            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2669        #microstrain               
2670        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2671            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2672        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2673            H = np.array(refl[:3])
2674            P = np.array(calcControls[phfx+'MustrainAxis'])
2675            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2676            Si = parmDict[phfx+'Mustrain;i']
2677            Sa = parmDict[phfx+'Mustrain;a']
2678            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2679        else:       #generalized - P.W. Stephens model  OK
2680            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2681            Sum = 0
2682            for i,strm in enumerate(Strms):
2683                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2684            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2685           
2686    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2687    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2688    sig /= ateln2
2689    return sig,gam
2690       
2691def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2692    'Needs a doc string'
2693    gamDict = {}
2694    sigDict = {}
2695    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2696        costh = cosd(refl[5+im]/2.)
2697        tanth = tand(refl[5+im]/2.)
2698        #crystallite size derivatives
2699        if calcControls[phfx+'SizeType'] == 'isotropic':
2700            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2701            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
2702            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2703        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2704            H = np.array(refl[:3])
2705            P = np.array(calcControls[phfx+'SizeAxis'])
2706            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2707            Si = parmDict[phfx+'Size;i']
2708            Sa = parmDict[phfx+'Size;a']
2709            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2710            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2711            Sgam = gami*sqtrm
2712            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2713            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2714            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2715            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2716            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2717            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2718        else:           #ellipsoidal crystallites
2719            const = 1.8*wave/(np.pi*costh)
2720            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2721            H = np.array(refl[:3])
2722            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2723            Sgam = const/lenR
2724            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2725                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2726                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2727        gamDict[phfx+'Size;mx'] = Sgam
2728        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2729               
2730        #microstrain derivatives               
2731        if calcControls[phfx+'MustrainType'] == 'isotropic':
2732            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2733            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2734            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2735        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2736            H = np.array(refl[:3])
2737            P = np.array(calcControls[phfx+'MustrainAxis'])
2738            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2739            Si = parmDict[phfx+'Mustrain;i']
2740            Sa = parmDict[phfx+'Mustrain;a']
2741            gami = 0.018*Si*Sa*tanth/np.pi
2742            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2743            Mgam = gami/sqtrm
2744            dsi = -gami*Si*cosP**2/sqtrm**3
2745            dsa = -gami*Sa*sinP**2/sqtrm**3
2746            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2747            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2748            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2749            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2750        else:       #generalized - P.W. Stephens model
2751            const = 0.018*refl[4+im]**2*tanth/np.pi
2752            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2753            Sum = 0
2754            for i,strm in enumerate(Strms):
2755                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2756                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2757                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2758            Mgam = const*np.sqrt(Sum)
2759            for i in range(len(Strms)):
2760                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2761                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2762        gamDict[phfx+'Mustrain;mx'] = Mgam
2763        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2764    else:   #'T'OF - All checked & OK
2765        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2766            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2767            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2768            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2769        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2770            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2771            H = np.array(refl[:3])
2772            P = np.array(calcControls[phfx+'SizeAxis'])
2773            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2774            Si = parmDict[phfx+'Size;i']
2775            Sa = parmDict[phfx+'Size;a']
2776            gami = const/(Si*Sa)
2777            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2778            Sgam = gami*sqtrm
2779            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2780            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2781            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2782            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2783            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2784            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2785        else:           #OK  ellipsoidal crystallites
2786            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2787            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2788            H = np.array(refl[:3])
2789            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2790            Sgam = const/lenR
2791            for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2792                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2793                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2794        gamDict[phfx+'Size;mx'] = Sgam  #OK
2795        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2796               
2797        #microstrain derivatives               
2798        if calcControls[phfx+'MustrainType'] == 'isotropic':
2799            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2800            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2801            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2802        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2803            H = np.array(refl[:3])
2804            P = np.array(calcControls[phfx+'MustrainAxis'])
2805            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2806            Si = parmDict[phfx+'Mustrain;i']
2807            Sa = parmDict[phfx+'Mustrain;a']
2808            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2809            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2810            Mgam = gami/sqtrm
2811            dsi = -gami*Si*cosP**2/sqtrm**3
2812            dsa = -gami*Sa*sinP**2/sqtrm**3
2813            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2814            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2815            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2816            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2817        else:       #generalized - P.W. Stephens model OK
2818            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2819            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2820            Sum = 0
2821            for i,strm in enumerate(Strms):
2822                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2823                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2824                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2825            Mgam = const*np.sqrt(Sum)
2826            for i in range(len(Strms)):
2827                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2828                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2829        gamDict[phfx+'Mustrain;mx'] = Mgam
2830        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2831       
2832    return sigDict,gamDict
2833       
2834def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2835    'Needs a doc string'
2836    if im:
2837        h,k,l,m = refl[:4]
2838        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2839        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2840    else:
2841        h,k,l = refl[:3]
2842        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2843    refl[4+im] = d
2844    if 'C' in calcControls[hfx+'histType']:
2845        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2846        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2847        if 'Bragg' in calcControls[hfx+'instType']:
2848            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2849                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2850        else:               #Debye-Scherrer - simple but maybe not right
2851            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2852    elif 'T' in calcControls[hfx+'histType']:
2853        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2854        #do I need sample position effects - maybe?
2855    return pos
2856
2857def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2858    'Needs a doc string'
2859    dpr = 180./np.pi
2860    if im:
2861        h,k,l,m = refl[:4]
2862        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2863        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2864        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
2865    else:
2866        m = 0
2867        h,k,l = refl[:3]       
2868        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2869    dst = np.sqrt(dstsq)
2870    dsp = 1./dst
2871    if 'C' in calcControls[hfx+'histType']:
2872        pos = refl[5+im]-parmDict[hfx+'Zero']
2873        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2874        dpdw = const*dst
2875        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2876        dpdZ = 1.0
2877        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2878            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2879        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2880        if 'Bragg' in calcControls[hfx+'instType']:
2881            dpdSh = -4.*shft*cosd(pos/2.0)
2882            dpdTr = -shft*sind(pos)*100.0
2883            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2884        else:               #Debye-Scherrer - simple but maybe not right
2885            dpdXd = -shft*cosd(pos)
2886            dpdYd = -shft*sind(pos)
2887            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2888    elif 'T' in calcControls[hfx+'histType']:
2889        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2890        dpdZ = 1.0
2891        dpdDC = dsp
2892        dpdDA = dsp**2
2893        dpdDB = 1./dsp
2894        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2895            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2896        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2897           
2898def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2899    'Needs a doc string'
2900    laue = SGData['SGLaue']
2901    uniq = SGData['SGUniq']
2902    h,k,l = refl[:3]
2903    if laue in ['m3','m3m']:
2904        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2905            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2906    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2907        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2908    elif laue in ['3R','3mR']:
2909        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2910    elif laue in ['4/m','4/mmm']:
2911        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2912    elif laue in ['mmm']:
2913        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2914    elif laue in ['2/m']:
2915        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2916        if uniq == 'a':
2917            Dij += parmDict[phfx+'D23']*k*l
2918        elif uniq == 'b':
2919            Dij += parmDict[phfx+'D13']*h*l
2920        elif uniq == 'c':
2921            Dij += parmDict[phfx+'D12']*h*k
2922    else:
2923        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2924            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2925    if 'C' in calcControls[hfx+'histType']:
2926        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2927    else:
2928        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2929           
2930def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2931    'Needs a doc string'
2932    laue = SGData['SGLaue']
2933    uniq = SGData['SGUniq']
2934    h,k,l = refl[:3]
2935    if laue in ['m3','m3m']:
2936        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2937            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2938    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2939        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2940    elif laue in ['3R','3mR']:
2941        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2942    elif laue in ['4/m','4/mmm']:
2943        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2944    elif laue in ['mmm']:
2945        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2946    elif laue in ['2/m']:
2947        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2948        if uniq == 'a':
2949            dDijDict[phfx+'D23'] = k*l
2950        elif uniq == 'b':
2951            dDijDict[phfx+'D13'] = h*l
2952        elif uniq == 'c':
2953            dDijDict[phfx+'D12'] = h*k
2954    else:
2955        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2956            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2957    if 'C' in calcControls[hfx+'histType']:
2958        for item in dDijDict:
2959            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2960    else:
2961        for item in dDijDict:
2962            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2963    return dDijDict
2964   
2965def GetDij(phfx,SGData,parmDict):
2966    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2967    return G2spc.HStrainVals(HSvals,SGData)
2968               
2969def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2970    '''Compute the observed structure factors for Powder histograms and store in reflection array
2971    Multiprocessing support added
2972    '''
2973    if GSASIIpath.GetConfigValue('Show_timing',False):
2974        starttime = time.time() #; print 'start GetFobsSq'
2975    histoList = list(Histograms.keys())
2976    histoList.sort()
2977    Ka2 = shl = lamRatio = kRatio = None
2978    for histogram in histoList:
2979        if 'PWDR' in histogram[:4]:
2980            Histogram = Histograms[histogram]
2981            hId = Histogram['hId']
2982            hfx = ':%d:'%(hId)
2983            Limits = calcControls[hfx+'Limits']
2984            if 'C' in calcControls[hfx+'histType']:
2985                shl = max(parmDict[hfx+'SH/L'],0.0005)
2986                Ka2 = False
2987                kRatio = 0.0
2988                if hfx+'Lam1' in list(parmDict.keys()):
2989                    Ka2 = True
2990                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2991                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
2992            x,y,w,yc,yb,yd = Histogram['Data']
2993            xMask = ma.getmaskarray(x)
2994            xB = np.searchsorted(x,Limits[0])
2995            xF = np.searchsorted(x,Limits[1])
2996            ymb = np.array(y-yb)
2997            ymb = np.where(ymb,ymb,1.0)
2998            ycmb = np.array(yc-yb)
2999            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3000            refLists = Histogram['Reflection Lists']
3001            for phase in refLists:
3002                if phase not in Phases:     #skips deleted or renamed phases silently!
3003                    continue
3004                Phase = Phases[phase]
3005                im = 0
3006                if Phase['General'].get('Modulated',False):
3007                    im = 1
3008                pId = Phase['pId']
3009                phfx = '%d:%d:'%(pId,hId)
3010                refDict = refLists[phase]
3011                sumFo = 0.0
3012                sumdF = 0.0
3013                sumFosq = 0.0
3014                sumdFsq = 0.0
3015                sumInt = 0.0
3016                nExcl = 0
3017                # test to see if we are using multiprocessing below
3018                useMP,ncores = G2mp.InitMP()
3019                if len(refDict['RefList']) < 100: useMP = False       
3020                if useMP: # multiprocessing: create a set of initialized Python processes
3021                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3022                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3023                    profArgs = [[] for i in range(G2mp.ncores)]
3024                else:
3025                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3026                if 'C' in calcControls[hfx+'histType']:
3027                    # are we multiprocessing?
3028                    for iref,refl in enumerate(refDict['RefList']):
3029                        if useMP: 
3030                            profArgs[iref%G2mp.ncores].append((refl,iref))
3031                        else:
3032                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3033                            if type(icod) is tuple:
3034                                refl[8+im] = icod[0]
3035                                sumInt += icod[1]
3036                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3037                            elif icod == -1:
3038                                refl[3+im] *= -1
3039                                nExcl += 1
3040                            elif icod == -2:
3041                                break
3042                    if useMP:
3043                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3044                            sumInt += sInt
3045                            for refl8im,irefl in resList:
3046                                if refl8im is None:
3047                                    refDict['RefList'][irefl][3+im] *= -1
3048                                    nExcl += 1
3049                                else:
3050                                    refDict['RefList'][irefl][8+im] = refl8im
3051                                    if parmDict[phfx+'LeBail']:
3052                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3053                elif 'T' in calcControls[hfx+'histType']:
3054                    for iref,refl in enumerate(refDict['RefList']):
3055                        if useMP: 
3056                            profArgs[iref%G2mp.ncores].append((refl,iref))
3057                        else:
3058                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3059                            if type(icod) is tuple:
3060                                refl[8+im] = icod[0]
3061                                sumInt += icod[1]
3062                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3063                            elif icod == -1:
3064                                refl[3+im] *= -1
3065                                nExcl += 1
3066                            elif icod == -2:
3067                                break
3068                    if useMP:
3069                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3070                            sumInt += sInt
3071                            for refl8im,irefl in resList:
3072                                if refl8im is None:
3073                                    refDict['RefList'][irefl][3+im] *= -1
3074                                    nExcl += 1
3075                                else:
3076                                    refDict['RefList'][irefl][8+im] = refl8im
3077                                    if parmDict[phfx+'LeBail']:
3078                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3079                if useMP: MPpool.terminate()
3080                sumFo = 0.0
3081                sumdF = 0.0
3082                sumFosq = 0.0
3083                sumdFsq = 0.0
3084                for iref,refl in enumerate(refDict['RefList']):
3085                    Fo = np.sqrt(np.abs(refl[8+im]))
3086                    Fc = np.sqrt(np.abs(refl[9]+im))
3087                    sumFo += Fo
3088                    sumFosq += refl[8+im]**2
3089                    sumdF += np.abs(Fo-Fc)
3090                    sumdFsq += (refl[8+im]-refl[9+im])**2
3091                if sumFo:
3092                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3093                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3094                else:
3095                    Histogram['Residuals'][phfx+'Rf'] = 100.
3096                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3097                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3098                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3099                Histogram['Residuals']['hId'] = hId
3100        elif 'HKLF' in histogram[:4]:
3101            Histogram = Histograms[histogram]
3102            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3103    if GSASIIpath.GetConfigValue('Show_timing',False):
3104        print ('GetFobsSq t=',time.time()-starttime)
3105               
3106def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3107    'Computes the powder pattern for a histogram based on contributions from all used phases'
3108    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3109   
3110    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3111        U = parmDict[hfx+'U']
3112        V = parmDict[hfx+'V']
3113        W = parmDict[hfx+'W']
3114        X = parmDict[hfx+'X']
3115        Y = parmDict[hfx+'Y']
3116        Z = parmDict[hfx+'Z']
3117        tanPos = tand(refl[5+im]/2.0)
3118        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3119        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3120        sig = max(0.001,sig)
3121        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3122        gam = max(0.001,gam)
3123        return sig,gam
3124               
3125    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3126        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3127            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3128        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3129        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3130        sig += Ssig
3131        gam += Sgam
3132        return sig,gam
3133       
3134    def GetReflAlpBet(refl,im,hfx,parmDict):
3135        alp = parmDict[hfx+'alpha']/refl[4+im]
3136        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3137        return alp,bet
3138       
3139    hId = Histogram['hId']
3140    hfx = ':%d:'%(hId)
3141    bakType = calcControls[hfx+'bakType']
3142    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
3143    yc = np.zeros_like(yb)
3144    cw = np.diff(ma.getdata(x))
3145    cw = np.append(cw,cw[-1])
3146       
3147    if 'C' in calcControls[hfx+'histType']:   
3148        shl = max(parmDict[hfx+'SH/L'],0.002)
3149        Ka2 = False
3150        if hfx+'Lam1' in (parmDict.keys()):
3151            wave = parmDict[hfx+'Lam1']
3152            Ka2 = True
3153            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3154            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3155        else:
3156            wave = parmDict[hfx+'Lam']
3157    else:
3158        shl = 0.
3159    for phase in Histogram['Reflection Lists']:
3160        refDict = Histogram['Reflection Lists'][phase]
3161        if phase not in Phases:     #skips deleted or renamed phases silently!