source: trunk/GSASIIstrMath.py @ 3755

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

change delt-M for mag. derivative calcs to 1.e-6 - same results.
use this code for single crystal magnetism as well (untestes)

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