source: trunk/GSASIIstrMath.py @ 4824

Last change on this file since 4824 was 4824, checked in by vondreele, 9 months ago

add new refinable RB parameter for atom site fraction - always there; can be used in constraints.. Available for residue & vector style RBs
fixes to delete RB routines
fixes to RB naming when adding RB to structure
Implement use of view position for new RB origin

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