source: trunk/GSASIIstrMath.py @ 3808

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

latest attempt at incommensurate mag str fctr & some cleanup of dead ends

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