source: trunk/GSASIIstrMath.py @ 3245

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

fix problem with incommensurate magnetic structures from CIFhklReader
make "Edit Body" as "Edit Residue Body" in RB menu
change protein validation bar plots to be vertical
change 'twin' to 'flag' in Reflection List table & add '-2' as 'free' (for proteins - not yet implemented)
fix PDB phase import to get atom type from 76:78 of ATOM record
change importer names for single crystal TOF data to be more explicit (SNS vs ISIS)
change cif structure factor importer - F2, sig(F2) & Fcalc now allowed

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