source: trunk/GSASIIstrMath.py @ 3825

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

create csv exporter for small angle data - includes size distribution histogram if any.
fix formatting issue with size distribution axis label
a bit different incomm. mag str fctr. math - still wrong

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