source: trunk/GSASIIstrMath.py @ 4145

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

mag. stuff again

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 218.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2019-09-13 18:57:37 +0000 (Fri, 13 Sep 2019) $
8# $Author: vondreele $
9# $Revision: 4145 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4145 2019-09-13 18:57:37Z 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: 4145 $")
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
1054    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1055    Kdata /= Kmean     #Cartesian unit vectors
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,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        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1435        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1436        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1437        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1438        if SGData['SGInv']: #centrosymmetric; B=0
1439            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1440        else:               
1441            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1442#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1443#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1444        iBeg += blkSize
1445    print (' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0))
1446    #loop over atoms - each dict entry is list of derivatives for all the reflections
1447    for i in range(len(Mdata)):     #these all OK
1448        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1449        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1450        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1451        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1452        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1453        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1454        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1455        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1456        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1457        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1458        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1459    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1460    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1461    for i in range(nTwin):
1462        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1463    return dFdvDict
1464   
1465def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1466    '''
1467    Compute super structure factors for all h,k,l,m for phase - no twins
1468    puts the result, F^2, in each ref[9] in refList
1469    works on blocks of 32 reflections for speed
1470    input:
1471   
1472    :param dict refDict: where
1473        'RefList' list where each ref = h,k,l,m,it,d,...
1474        'FF' dict of form factors - filed in below
1475    :param np.array G:      reciprocal metric tensor
1476    :param str pfx:    phase id string
1477    :param dict SGData: space group info. dictionary output from SpcGroup
1478    :param dict calcControls:
1479    :param dict ParmDict:
1480
1481    '''
1482    phfx = pfx.split(':')[0]+hfx
1483    ast = np.sqrt(np.diag(G))
1484    GS = G/np.outer(ast,ast)
1485    uAmat,uBmat = G2lat.Gmat2AB(GS)
1486    Mast = twopisq*np.multiply.outer(ast,ast)   
1487    SGInv = SGData['SGInv']
1488    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1489    Ncen = len(SGData['SGCen'])
1490    Nops = len(SGMT)*(1+SGData['SGInv'])
1491    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1492    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1493    SSCen = SSGData['SSGCen']
1494    FFtables = calcControls['FFtables']
1495    BLtables = calcControls['BLtables']
1496    MFtables = calcControls['MFtables']
1497    Amat,Bmat = G2lat.Gmat2AB(G)
1498    Flack = 1.0
1499    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1500        Flack = 1.-2.*parmDict[phfx+'Flack']
1501    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1502        GetAtomFXU(pfx,calcControls,parmDict)
1503    if not Xdata.size:          #no atoms in phase!
1504        return
1505    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1506    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1507    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1508
1509    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms
1510       
1511        mXYZ = np.array([[xyz[0] for xyz in list(G2spc.GenAtom(xyz,SGData,All=True,Move=True))] for xyz in (Xdata+dXdata).T])%1. #Natn,Nop,xyz
1512        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz sum matches drawing
1513        Mmod = MmodA+MmodB
1514       
1515        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1516            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1517            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1518                GSdata = np.hstack((GSdata,-GSdata))     
1519            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1520            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1521            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1522            Mmod += GSdata[nxs,:,:,:]
1523           
1524        Kdata = np.inner(Mmod,uAmat)    #Ntau,Nop,Natm,Mxyz
1525        Mag = np.sqrt(np.sum(Kdata**2,axis=-1))
1526        Kdata /= Mag[:,:,:,nxs]     #Cartesian unit vectors
1527        Kdata = np.nan_to_num(Kdata)    #Ntau,Nops,Natm,Mxyz       
1528
1529    FF = np.zeros(len(Tdata))
1530    if 'NC' in calcControls[hfx+'histType']:
1531        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1532    elif 'X' in calcControls[hfx+'histType']:
1533        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1534        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1535    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1536    bij = Mast*Uij
1537    blkSize = 48       #no. of reflections in a block
1538    nRef = refDict['RefList'].shape[0]
1539    SQ = 1./(2.*refDict['RefList'].T[5])**2
1540    if 'N' in calcControls[hfx+'histType']:
1541        dat = G2el.getBLvalues(BLtables)
1542        refDict['FF']['El'] = list(dat.keys())
1543        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1544        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1545        for iel,El in enumerate(refDict['FF']['El']):
1546            if El in MFtables:
1547                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1548    else:
1549        dat = G2el.getFFvalues(FFtables,0.)
1550        refDict['FF']['El'] = list(dat.keys())
1551        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1552        for iel,El in enumerate(refDict['FF']['El']):
1553            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1554#    time0 = time.time()
1555#reflection processing begins here - big arrays!
1556    iBeg = 0
1557    while iBeg < nRef:
1558        iFin = min(iBeg+blkSize,nRef)
1559        mRef= iFin-iBeg
1560        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1561        H = refl.T[:4]                          #array(blkSize,4)
1562        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1563        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1564        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1565        Uniq = np.inner(H.T,SSGMT)
1566        UniqP = np.inner(HP.T,SGMT)
1567        Phi = np.inner(H.T,SSGT)
1568        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1569            Uniq = np.hstack((Uniq,-Uniq))
1570            Phi = np.hstack((Phi,-Phi))
1571            UniqP = np.hstack((UniqP,-UniqP))
1572        if 'T' in calcControls[hfx+'histType']:
1573            if 'P' in calcControls[hfx+'histType']:
1574                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1575            else:
1576                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1577            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1578            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1579        Bab = 0.
1580        if phfx+'BabA' in parmDict:
1581            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1582        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1583        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1584        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1585        phase = np.hstack([phase for cen in SSCen])
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           
1596            phasem = twopi*np.inner(H.T[:,:3],mXYZ)
1597            phasem = np.swapaxes(phasem,1,2)            #Nref,Nops,Natm
1598            cosm = np.cos(phasem)
1599            sinm = np.sin(phasem)
1600            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1601            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Mdata*Fdata*MF/(2.*Nops)     #Nref,Natm
1602                     
1603            HM = np.inner(Bmat,HP.T)                            #put into cartesian space X||H,Z||H*L
1604            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T               #& normalize    Nref,hkl
1605            eDotK = np.sum(eM[:,nxs,nxs,nxs,:]*Kdata[nxs,:,:,:,:],axis=-1)    #Nref,Ntau,Nops,Natm
1606#for fixed moments --> m=0 reflections
1607            fam0 = 0.
1608            fbm0 = 0.
1609            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1610                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1611                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]   
1612#for modulated moments --> m != 0 reflections
1613       
1614            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]-    \
1615                np.sign(H[3,i])*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1616                       
1617            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+    \
1618                np.sign(H[3,i])*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
1619           
1620            if not SGData['SGGray']:
1621                fams *= 0.5
1622                fbms *= 0.5
1623                fams += fam0[:,nxs,:,:,:]
1624                fbms += fbm0[:,nxs,:,:,:]
1625# do sum on ops, atms 1st                       
1626            cos2 = np.sqrt(1.-eDotK**2)
1627           
1628            famqs = np.sum(np.sum(fams*cos2[:,:,:,:,nxs],axis=-2),axis=-2)      #Nref,Ntau,Mxyz; sum ops & atoms
1629            fbmqs = np.sum(np.sum(fbms*cos2[:,:,:,:,nxs],axis=-2),axis=-2)
1630           
1631            fass = np.sum(famqs**2,axis=-1)      #mag intensity calc F^2-(e.F)^2
1632            fbss = np.sum(fbmqs**2,axis=-1)
1633           
1634#do integration           
1635            fas = np.sum(glWt*fass,axis=1)
1636            fbs = np.sum(glWt*fbss,axis=1)
1637           
1638            refl.T[10] = fas+fbs
1639            refl.T[11] = atan2d(fbs,fas)
1640
1641        else:
1642            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1643            if 'T' in calcControls[hfx+'histType']:
1644                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1645                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1646            else:
1647                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1648                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1649            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1650            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1651            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1652            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1653           
1654            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1655            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1656        if 'P' not in calcControls[hfx+'histType']:
1657            refl.T[8] = np.copy(refl.T[10])               
1658        iBeg += blkSize
1659#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1660    return copy.deepcopy(refDict['RefList'])
1661
1662def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1663    '''
1664    Compute super structure factors for all h,k,l,m for phase - twins only
1665    puts the result, F^2, in each ref[8+im] in refList
1666    works on blocks of 32 reflections for speed
1667    input:
1668   
1669    :param dict refDict: where
1670        'RefList' list where each ref = h,k,l,m,it,d,...
1671        'FF' dict of form factors - filed in below
1672    :param np.array G:      reciprocal metric tensor
1673    :param str pfx:    phase id string
1674    :param dict SGData: space group info. dictionary output from SpcGroup
1675    :param dict calcControls:
1676    :param dict ParmDict:
1677
1678    '''
1679    phfx = pfx.split(':')[0]+hfx
1680    ast = np.sqrt(np.diag(G))
1681    Mast = twopisq*np.multiply.outer(ast,ast)   
1682    SGInv = SGData['SGInv']
1683    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1684    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1685    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1686    FFtables = calcControls['FFtables']
1687    BLtables = calcControls['BLtables']
1688    MFtables = calcControls['MFtables']
1689    Flack = 1.0
1690    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1691        Flack = 1.-2.*parmDict[phfx+'Flack']
1692    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1693    TwDict = refDict.get('TwDict',{})           
1694    if 'S' in calcControls[hfx+'histType']:
1695        NTL = calcControls[phfx+'NTL']
1696        NM = calcControls[phfx+'TwinNMN']+1
1697        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1698        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1699        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1700    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1701        GetAtomFXU(pfx,calcControls,parmDict)
1702    if not Xdata.size:          #no atoms in phase!
1703        return
1704    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1705    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1706    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1707    FF = np.zeros(len(Tdata))
1708    if 'NC' in calcControls[hfx+'histType']:
1709        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1710    elif 'X' in calcControls[hfx+'histType']:
1711        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1712        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1713    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1714    bij = Mast*Uij
1715    blkSize = 32       #no. of reflections in a block
1716    nRef = refDict['RefList'].shape[0]
1717    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1718        SQ = 1./(2.*refDict['RefList'].T[5])**2
1719        if 'N' in calcControls[hfx+'histType']:
1720            dat = G2el.getBLvalues(BLtables)
1721            refDict['FF']['El'] = list(dat.keys())
1722            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1723            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1724            for iel,El in enumerate(refDict['FF']['El']):
1725                if El in MFtables:
1726                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1727        else:
1728            dat = G2el.getFFvalues(FFtables,0.)
1729            refDict['FF']['El'] = list(dat.keys())
1730            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1731            for iel,El in enumerate(refDict['FF']['El']):
1732                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1733#    time0 = time.time()
1734#reflection processing begins here - big arrays!
1735    iBeg = 0
1736    while iBeg < nRef:
1737        iFin = min(iBeg+blkSize,nRef)
1738        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1739        H = refl[:,:4]                          #array(blkSize,4)
1740        H3 = refl[:,:3]
1741        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1742        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1743        H3 = np.inner(H3,TwinLaw)       
1744        TwMask = np.any(HP,axis=-1)
1745        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1746            for ir in range(blkSize):
1747                iref = ir+iBeg
1748                if iref in TwDict:
1749                    for i in TwDict[iref]:
1750                        for n in range(NTL):
1751                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1752                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1753            TwMask = np.any(HP,axis=-1)
1754        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1755        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1756        Uniq = np.inner(H,SSGMT)
1757        Uniq3 = np.inner(H3,SGMT)
1758        UniqP = np.inner(HP,SGMT)
1759        Phi = np.inner(H,SSGT)
1760        if SGInv:   #if centro - expand HKL sets
1761            Uniq = np.hstack((Uniq,-Uniq))
1762            Uniq3 = np.hstack((Uniq3,-Uniq3))
1763            Phi = np.hstack((Phi,-Phi))
1764            UniqP = np.hstack((UniqP,-UniqP))
1765        if 'T' in calcControls[hfx+'histType']:
1766            if 'P' in calcControls[hfx+'histType']:
1767                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1768            else:
1769                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1770            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1771            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1772        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1773        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1774        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1775        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1776        sinp = np.sin(phase)
1777        cosp = np.cos(phase)
1778        biso = -SQfactor*Uisodata[:,nxs]
1779        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1780        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1781        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1782        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1783        if 'T' in calcControls[hfx+'histType']:
1784            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1785            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1786        else:
1787            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1788            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1789        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1790        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1791        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1792        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1793        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1794        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1795        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1796            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1797        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1798        iBeg += blkSize
1799#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1800
1801def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1802    '''
1803    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1804    Only Fourier component are done analytically here
1805    input:
1806   
1807    :param dict refDict: where
1808        'RefList' list where each ref = h,k,l,m,it,d,...
1809        'FF' dict of form factors - filled in below
1810    :param int im: = 1 (could be eliminated)
1811    :param np.array G:      reciprocal metric tensor
1812    :param str hfx:    histogram id string
1813    :param str pfx:    phase id string
1814    :param dict SGData: space group info. dictionary output from SpcGroup
1815    :param dict SSGData: super space group info.
1816    :param dict calcControls:
1817    :param dict ParmDict:
1818   
1819    :returns: dict dFdvDict: dictionary of derivatives
1820    '''
1821    phfx = pfx.split(':')[0]+hfx
1822    ast = np.sqrt(np.diag(G))
1823    Mast = twopisq*np.multiply.outer(ast,ast)
1824    SGInv = SGData['SGInv']
1825    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1826    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1827    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1828    FFtables = calcControls['FFtables']
1829    BLtables = calcControls['BLtables']
1830    nRef = len(refDict['RefList'])
1831    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1832        GetAtomFXU(pfx,calcControls,parmDict)
1833    if not Xdata.size:          #no atoms in phase!
1834        return {}
1835    mSize = len(Mdata)  #no. atoms
1836    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1837    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1838    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1839    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1840    FF = np.zeros(len(Tdata))
1841    if 'NC' in calcControls[hfx+'histType']:
1842        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1843    elif 'X' in calcControls[hfx+'histType']:
1844        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1845        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1846    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1847    bij = Mast*Uij
1848    if not len(refDict['FF']):
1849        if 'N' in calcControls[hfx+'histType']:
1850            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1851        else:
1852            dat = G2el.getFFvalues(FFtables,0.)       
1853        refDict['FF']['El'] = list(dat.keys())
1854        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1855    dFdvDict = {}
1856    dFdfr = np.zeros((nRef,mSize))
1857    dFdx = np.zeros((nRef,mSize,3))
1858    dFdui = np.zeros((nRef,mSize))
1859    dFdua = np.zeros((nRef,mSize,6))
1860    dFdbab = np.zeros((nRef,2))
1861    dFdfl = np.zeros((nRef))
1862    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1863    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1864    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1865    Flack = 1.0
1866    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1867        Flack = 1.-2.*parmDict[phfx+'Flack']
1868    time0 = time.time()
1869    nRef = len(refDict['RefList'])/100
1870    for iref,refl in enumerate(refDict['RefList']):
1871        if 'T' in calcControls[hfx+'histType']:
1872            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1873        H = np.array(refl[:4])
1874        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1875        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1876        SQfactor = 8.0*SQ*np.pi**2
1877        Bab = 0.0
1878        if phfx+'BabA' in parmDict:
1879            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1880            Bab = parmDict[phfx+'BabA']*dBabdA
1881        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1882        FF = refDict['FF']['FF'][iref].T[Tindx]
1883        Uniq = np.inner(H,SSGMT)
1884        Phi = np.inner(H,SSGT)
1885        UniqP = np.inner(HP,SGMT)
1886        if SGInv:   #if centro - expand HKL sets
1887            Uniq = np.vstack((Uniq,-Uniq))
1888            Phi = np.hstack((Phi,-Phi))
1889            UniqP = np.vstack((UniqP,-UniqP))
1890        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1891        sinp = np.sin(phase)
1892        cosp = np.cos(phase)
1893        occ = Mdata*Fdata/Uniq.shape[0]
1894        biso = -SQfactor*Uisodata[:,nxs]
1895        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1896        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1897        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1898        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1899        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1900        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1901        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1902        fotp = FPP*Tcorr            #ops x atoms
1903        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1904        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1905        # GfpuA is 2 x ops x atoms
1906        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1907        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1908        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1909        fag = fa*GfpuA[0]-fb*GfpuA[1]
1910        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1911       
1912        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1913        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1914        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1915        fbx = np.array([fot*cosp,-fotp*sinp])
1916        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1917        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1918        #sum below is over Uniq
1919        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1920        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1921        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1922        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1923        dfadui = np.sum(-SQfactor*fag,axis=1)
1924        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1925        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1926        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1927        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1928        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1929        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1930        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1931        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1932        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1933        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1934        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1935        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1936        if not SGData['SGInv']:   #Flack derivative
1937            dfadfl = np.sum(-FPP*Tcorr*sinp)
1938            dfbdfl = np.sum(FPP*Tcorr*cosp)
1939        else:
1940            dfadfl = 1.0
1941            dfbdfl = 1.0
1942        SA = fas[0]+fas[1]      #float = A+A'
1943        SB = fbs[0]+fbs[1]      #float = B+B'
1944        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1945            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1946            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1947                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1948            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1949                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1950            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1951                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1952            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1953                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1954            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1955                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1956            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1957                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1958            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1959                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1960        else:                       #OK, I think
1961            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1962            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1963            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1964            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1965            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1966                           
1967            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1968            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1969            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1970        if phfx+'BabA' in parmDict:
1971            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1972                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1973        #loop over atoms - each dict entry is list of derivatives for all the reflections
1974        if not iref%100 :
1975            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1976    for i in range(len(Mdata)):     #loop over atoms
1977        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1978        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1979        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1980        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1981        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1982        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1983        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1984        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1985        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1986        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1987        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1988        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1989            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1990            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1991        nx = 0
1992        if waveTypes[i] in ['Block','ZigZag']:
1993            nx = 1 
1994        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1995            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1996            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1997            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1998            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1999            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2000            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2001        for j in range(USSdata.shape[1]):       #loop over waves
2002            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2003            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2004            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2005            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2006            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2007            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2008            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2009            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2010            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2011            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2012            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2013            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2014           
2015    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
2016    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2017    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2018    return dFdvDict
2019
2020def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2021    '''
2022    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
2023    input:
2024   
2025    :param dict refDict: where
2026        'RefList' list where each ref = h,k,l,m,it,d,...
2027        'FF' dict of form factors - filled in below
2028    :param int im: = 1 (could be eliminated)
2029    :param np.array G:      reciprocal metric tensor
2030    :param str hfx:    histogram id string
2031    :param str pfx:    phase id string
2032    :param dict SGData: space group info. dictionary output from SpcGroup
2033    :param dict SSGData: super space group info.
2034    :param dict calcControls:
2035    :param dict ParmDict:
2036   
2037    :returns: dict dFdvDict: dictionary of derivatives
2038    '''
2039
2040    trefDict = copy.deepcopy(refDict)
2041    dM = 1.e-4
2042    dFdvDict = {}
2043    for parm in parmDict:
2044        if ':' not in parm:
2045            continue
2046        if parm.split(':')[2] in ['Tmin','Tmax','Xmax','Ymax','Zmax','Fzero','Fwid',
2047            'MXsin','MXcos','MYsin','MYcos','MZsin','MZcos','AMx','AMy','AMz',]:
2048            parmDict[parm] += dM
2049            prefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2050            parmDict[parm] -= 2*dM
2051            mrefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2052            parmDict[parm] += dM
2053            dFdvDict[parm] = (prefList[:,9+im]-mrefList[:,9+im])/(2.*dM)
2054    return dFdvDict
2055   
2056def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2057    'Needs a doc string'
2058    phfx = pfx.split(':')[0]+hfx
2059    ast = np.sqrt(np.diag(G))
2060    Mast = twopisq*np.multiply.outer(ast,ast)
2061    SGInv = SGData['SGInv']
2062    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2063    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2064    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2065    FFtables = calcControls['FFtables']
2066    BLtables = calcControls['BLtables']
2067    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2068    TwDict = refDict.get('TwDict',{})           
2069    if 'S' in calcControls[hfx+'histType']:
2070        NTL = calcControls[phfx+'NTL']
2071        NM = calcControls[phfx+'TwinNMN']+1
2072        TwinLaw = calcControls[phfx+'TwinLaw']
2073        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2074    nTwin = len(TwinLaw)       
2075    nRef = len(refDict['RefList'])
2076    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2077        GetAtomFXU(pfx,calcControls,parmDict)
2078    if not Xdata.size:          #no atoms in phase!
2079        return {}
2080    mSize = len(Mdata)  #no. atoms
2081    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2082    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
2083    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2084    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2085    FF = np.zeros(len(Tdata))
2086    if 'NC' in calcControls[hfx+'histType']:
2087        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2088    elif 'X' in calcControls[hfx+'histType']:
2089        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2090        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2091    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2092    bij = Mast*Uij
2093    if not len(refDict['FF']):
2094        if 'N' in calcControls[hfx+'histType']:
2095            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2096        else:
2097            dat = G2el.getFFvalues(FFtables,0.)       
2098        refDict['FF']['El'] = list(dat.keys())
2099        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2100    dFdvDict = {}
2101    dFdfr = np.zeros((nRef,nTwin,mSize))
2102    dFdx = np.zeros((nRef,nTwin,mSize,3))
2103    dFdui = np.zeros((nRef,nTwin,mSize))
2104    dFdua = np.zeros((nRef,nTwin,mSize,6))
2105    dFdbab = np.zeros((nRef,nTwin,2))
2106    dFdtw = np.zeros((nRef,nTwin))
2107    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2108    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2109    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2110    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2111    Flack = 1.0
2112    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2113        Flack = 1.-2.*parmDict[phfx+'Flack']
2114    time0 = time.time()
2115    nRef = len(refDict['RefList'])/100
2116    for iref,refl in enumerate(refDict['RefList']):
2117        if 'T' in calcControls[hfx+'histType']:
2118            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2119        H = np.array(refl[:4])
2120        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2121        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2122        TwMask = np.any(H,axis=-1)
2123        if TwinLaw.shape[0] > 1 and TwDict:
2124            if iref in TwDict:
2125                for i in TwDict[iref]:
2126                    for n in range(NTL):
2127                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2128            TwMask = np.any(H,axis=-1)
2129        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2130        SQfactor = 8.0*SQ*np.pi**2
2131        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2132        Bab = parmDict[phfx+'BabA']*dBabdA
2133        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2134        FF = refDict['FF']['FF'][iref].T[Tindx]
2135        Uniq = np.inner(H,SSGMT)
2136        Phi = np.inner(H,SSGT)
2137        UniqP = np.inner(HP,SGMT)
2138        if SGInv:   #if centro - expand HKL sets
2139            Uniq = np.vstack((Uniq,-Uniq))
2140            Phi = np.hstack((Phi,-Phi))
2141            UniqP = np.vstack((UniqP,-UniqP))
2142        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2143        sinp = np.sin(phase)
2144        cosp = np.cos(phase)
2145        occ = Mdata*Fdata/Uniq.shape[0]
2146        biso = -SQfactor*Uisodata[:,nxs]
2147        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2148        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2149        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2150        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2151        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2152        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2153        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2154        fotp = FPP*Tcorr            #ops x atoms
2155        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2156        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2157        # GfpuA is 2 x ops x atoms
2158        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2159        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2160        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2161        fag = fa*GfpuA[0]-fb*GfpuA[1]
2162        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2163       
2164        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2165        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2166        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2167        fbx = np.array([fot*cosp,-fotp*sinp])
2168        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2169        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2170        #sum below is over Uniq
2171        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2172        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2173        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2174        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2175        dfadui = np.sum(-SQfactor*fag,axis=1)
2176        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2177        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2178        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2179        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2180        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2181        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2182        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2183        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2184        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2185        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2186        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2187        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2188        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2189        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2190#        GSASIIpath.IPyBreak()
2191        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2192        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2193        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2194        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)]
2195        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)]
2196        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)]
2197        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)]
2198        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2199
2200        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2201        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2202        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2203        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2204#            GSASIIpath.IPyBreak()
2205        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2206            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2207        #loop over atoms - each dict entry is list of derivatives for all the reflections
2208        if not iref%100 :
2209            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2210    for i in range(len(Mdata)):     #loop over atoms
2211        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2212        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2213        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2214        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2215        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2216        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2217        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2218        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2219        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2220        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2221        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2222        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2223            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2224            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2225        nx = 0
2226        if waveTypes[i] in ['Block','ZigZag']:
2227            nx = 1 
2228            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2229            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2230            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2231            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2232            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2233        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2234            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2235            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2236            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2237            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2238            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2239            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2240        for j in range(USSdata.shape[1]):       #loop over waves
2241            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2242            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2243            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2244            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2245            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2246            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2247            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2248            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2249            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2250            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2251            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2252            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2253           
2254#        GSASIIpath.IPyBreak()
2255    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2256    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2257    return dFdvDict
2258   
2259def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2260    ''' Single crystal extinction function; returns extinction & derivative
2261    '''
2262    extCor = 1.0
2263    dervDict = {}
2264    dervCor = 1.0
2265    if calcControls[phfx+'EType'] != 'None':
2266        SQ = 1/(4.*ref[4+im]**2)
2267        if 'C' in parmDict[hfx+'Type']:           
2268            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2269        else:   #'T'
2270            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2271        if 'SXC' in parmDict[hfx+'Type']:
2272            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2273            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2274            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2275            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2276        elif 'SNT' in parmDict[hfx+'Type']:
2277            AV = 1.e7/parmDict[pfx+'Vol']**2
2278            PL = SQ
2279            PLZ = AV*ref[9+im]*ref[12+im]**2
2280        elif 'SNC' in parmDict[hfx+'Type']:
2281            AV = 1.e7/parmDict[pfx+'Vol']**2
2282            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2283            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2284           
2285        if 'Primary' in calcControls[phfx+'EType']:
2286            PLZ *= 1.5
2287        else:
2288            if 'C' in parmDict[hfx+'Type']:
2289                PLZ *= calcControls[phfx+'Tbar']
2290            else: #'T'
2291                PLZ *= ref[13+im]      #t-bar
2292        if 'Primary' in calcControls[phfx+'EType']:
2293            PLZ *= 1.5
2294            PSIG = parmDict[phfx+'Ep']
2295        elif 'I & II' in calcControls[phfx+'EType']:
2296            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2297        elif 'Type II' in calcControls[phfx+'EType']:
2298            PSIG = parmDict[phfx+'Es']
2299        else:       # 'Secondary Type I'
2300            PSIG = parmDict[phfx+'Eg']/PL
2301           
2302        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2303        AL = 0.025+0.285*cos2T
2304        BG = 0.02-0.025*cos2T
2305        BL = 0.15-0.2*(0.75-cos2T)**2
2306        if cos2T < 0.:
2307            BL = -0.45*cos2T
2308        CG = 2.
2309        CL = 2.
2310        PF = PLZ*PSIG
2311       
2312        if 'Gaussian' in calcControls[phfx+'EApprox']:
2313            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2314            extCor = np.sqrt(PF4)
2315            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2316        else:
2317            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2318            extCor = np.sqrt(PF4)
2319            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2320
2321        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2322        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2323            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2324        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2325            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2326        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2327            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2328               
2329    return 1./extCor,dervDict,dervCor
2330   
2331def Dict2Values(parmdict, varylist):
2332    '''Use before call to leastsq to setup list of values for the parameters
2333    in parmdict, as selected by key in varylist'''
2334    return [parmdict[key] for key in varylist] 
2335   
2336def Values2Dict(parmdict, varylist, values):
2337    ''' Use after call to leastsq to update the parameter dictionary with
2338    values corresponding to keys in varylist'''
2339    parmdict.update(zip(varylist,values))
2340   
2341def GetNewCellParms(parmDict,varyList):
2342    '''Compute unit cell tensor terms from varied Aij and Dij values.
2343    Terms are included in the dict only if Aij or Dij is varied.
2344    '''
2345    newCellDict = {}
2346    Anames = ['A'+str(i) for i in range(6)]
2347    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2348    for item in varyList:
2349        keys = item.split(':')
2350        if keys[2] in Ddict:
2351            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2352            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2353            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2354    return newCellDict          # is e.g. {'0::D11':A0-D11}
2355   
2356def ApplyXYZshifts(parmDict,varyList):
2357    '''
2358    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2359   
2360    :param dict parmDict: parameter dictionary
2361    :param list varyList: list of variables (not used!)
2362    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2363
2364    '''
2365    newAtomDict = {}
2366    for item in parmDict:
2367        if 'dA' in item:
2368            parm = ''.join(item.split('d'))
2369            parmDict[parm] += parmDict[item]
2370            newAtomDict[item] = [parm,parmDict[parm]]
2371    return newAtomDict
2372   
2373def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2374    'Spherical harmonics texture'
2375    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2376    if 'T' in calcControls[hfx+'histType']:
2377        tth = parmDict[hfx+'2-theta']
2378    else:
2379        tth = refl[5+im]
2380    odfCor = 1.0
2381    H = refl[:3]
2382    cell = G2lat.Gmat2cell(g)
2383    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2384    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2385    phi,beta = G2lat.CrsAng(H,cell,SGData)
2386    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2387    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2388    for item in SHnames:
2389        L,M,N = eval(item.strip('C'))
2390        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2391        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2392        Lnorm = G2lat.Lnorm(L)
2393        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2394    return odfCor
2395   
2396def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2397    'Spherical harmonics texture derivatives'
2398    if 'T' in calcControls[hfx+'histType']:
2399        tth = parmDict[hfx+'2-theta']
2400    else:
2401        tth = refl[5+im]
2402    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2403    odfCor = 1.0
2404    dFdODF = {}
2405    dFdSA = [0,0,0]
2406    H = refl[:3]
2407    cell = G2lat.Gmat2cell(g)
2408    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2409    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2410    phi,beta = G2lat.CrsAng(H,cell,SGData)
2411    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2412    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2413    for item in SHnames:
2414        L,M,N = eval(item.strip('C'))
2415        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2416        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2417        Lnorm = G2lat.Lnorm(L)
2418        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2419        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2420        for i in range(3):
2421            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2422    return odfCor,dFdODF,dFdSA
2423   
2424def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2425    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2426    if 'T' in calcControls[hfx+'histType']:
2427        tth = parmDict[hfx+'2-theta']
2428    else:
2429        tth = refl[5+im]
2430    odfCor = 1.0
2431    H = refl[:3]
2432    cell = G2lat.Gmat2cell(g)
2433    Sangls = [0.,0.,0.]
2434    if 'Bragg' in calcControls[hfx+'instType']:
2435        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2436        IFCoup = True
2437    else:
2438        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2439        IFCoup = False
2440    phi,beta = G2lat.CrsAng(H,cell,SGData)
2441    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2442    SHnames = calcControls[phfx+'SHnames']
2443    for item in SHnames:
2444        L,N = eval(item.strip('C'))
2445        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2446        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2447        Lnorm = G2lat.Lnorm(L)
2448        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2449    return np.squeeze(odfCor)
2450   
2451def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2452    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2453    if 'T' in calcControls[hfx+'histType']:
2454        tth = parmDict[hfx+'2-theta']
2455    else:
2456        tth = refl[5+im]
2457    odfCor = 1.0
2458    dFdODF = {}
2459    H = refl[:3]
2460    cell = G2lat.Gmat2cell(g)
2461    Sangls = [0.,0.,0.]
2462    if 'Bragg' in calcControls[hfx+'instType']:
2463        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2464        IFCoup = True
2465    else:
2466        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2467        IFCoup = False
2468    phi,beta = G2lat.CrsAng(H,cell,SGData)
2469    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2470    SHnames = calcControls[phfx+'SHnames']
2471    for item in SHnames:
2472        L,N = eval(item.strip('C'))
2473        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2474        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2475        Lnorm = G2lat.Lnorm(L)
2476        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2477        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2478    return odfCor,dFdODF
2479   
2480def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2481    'March-Dollase preferred orientation correction'
2482    POcorr = 1.0
2483    MD = parmDict[phfx+'MD']
2484    if MD != 1.0:
2485        MDAxis = calcControls[phfx+'MDAxis']
2486        sumMD = 0
2487        for H in uniq:           
2488            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2489            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2490            sumMD += A**3
2491        POcorr = sumMD/len(uniq)
2492    return POcorr
2493   
2494def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2495    'Needs a doc string'
2496    POcorr = 1.0
2497    POderv = {}
2498    if calcControls[phfx+'poType'] == 'MD':
2499        MD = parmDict[phfx+'MD']
2500        MDAxis = calcControls[phfx+'MDAxis']
2501        sumMD = 0
2502        sumdMD = 0
2503        for H in uniq:           
2504            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2505            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2506            sumMD += A**3
2507            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2508        POcorr = sumMD/len(uniq)
2509        POderv[phfx+'MD'] = sumdMD/len(uniq)
2510    else:   #spherical harmonics
2511        if calcControls[phfx+'SHord']:
2512            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2513    return POcorr,POderv
2514   
2515def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2516    'Needs a doc string'
2517    if 'Debye' in calcControls[hfx+'instType']:
2518        if 'T' in calcControls[hfx+'histType']:
2519            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2520        else:
2521            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2522    else:
2523        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2524   
2525def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2526    'Needs a doc string'
2527    if 'Debye' in calcControls[hfx+'instType']:
2528        if 'T' in calcControls[hfx+'histType']:
2529            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2530        else:
2531            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2532    else:
2533        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2534       
2535def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2536    'Needs a doc string'
2537    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2538    pi2 = np.sqrt(2./np.pi)
2539    if 'T' in calcControls[hfx+'histType']:
2540        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2541        wave = refl[14+im]
2542    else:   #'C'W
2543        sth2 = sind(refl[5+im]/2.)**2
2544        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2545    c2th = 1.-2.0*sth2
2546    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2547    if 'X' in calcControls[hfx+'histType']:
2548        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2549    xfac = flv2*parmDict[phfx+'Extinction']
2550    exb = 1.0
2551    if xfac > -1.:
2552        exb = 1./np.sqrt(1.+xfac)
2553    exl = 1.0
2554    if 0 < xfac <= 1.:
2555        xn = np.array([xfac**(i+1) for i in range(6)])
2556        exl += np.sum(xn*coef)
2557    elif xfac > 1.:
2558        xfac2 = 1./np.sqrt(xfac)
2559        exl = pi2*(1.-0.125/xfac)*xfac2
2560    return exb*sth2+exl*(1.-sth2)
2561   
2562def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2563    'Needs a doc string'
2564    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2565    pi2 = np.sqrt(2./np.pi)
2566    if 'T' in calcControls[hfx+'histType']:
2567        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2568        wave = refl[14+im]
2569    else:   #'C'W
2570        sth2 = sind(refl[5+im]/2.)**2
2571        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2572    c2th = 1.-2.0*sth2
2573    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2574    if 'X' in calcControls[hfx+'histType']:
2575        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2576    xfac = flv2*parmDict[phfx+'Extinction']
2577    dbde = -500.*flv2
2578    if xfac > -1.:
2579        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2580    dlde = 0.
2581    if 0 < xfac <= 1.:
2582        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2583        dlde = np.sum(xn*coef)/xfac
2584    elif xfac > 1.:
2585        xfac2 = 1./np.sqrt(xfac)
2586        dlde = 0.5*flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2587       
2588    return dbde*sth2+dlde*(1.-sth2)
2589   
2590def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2591    'Needs a doc string'    #need powder extinction!
2592    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2593    if 'X' in parmDict[hfx+'Type']:
2594        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2595    POcorr = 1.0
2596    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2597        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2598    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2599        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2600    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2601        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2602    Icorr *= POcorr
2603    AbsCorr = 1.0
2604    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2605    Icorr *= AbsCorr
2606    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2607    Icorr *= ExtCorr
2608    return Icorr,POcorr,AbsCorr,ExtCorr
2609   
2610def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2611    'Needs a doc string'    #need powder extinction derivs!
2612    dIdsh = 1./parmDict[hfx+'Scale']
2613    dIdsp = 1./parmDict[phfx+'Scale']
2614    if 'X' in parmDict[hfx+'Type']:
2615        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2616        dIdPola /= pola
2617    else:       #'N'
2618        dIdPola = 0.0
2619    dFdODF = {}
2620    dFdSA = [0,0,0]
2621    dIdPO = {}
2622    if pfx+'SHorder' in parmDict:
2623        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2624        for iSH in dFdODF:
2625            dFdODF[iSH] /= odfCor
2626        for i in range(3):
2627            dFdSA[i] /= odfCor
2628    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2629        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2630        for iPO in dIdPO:
2631            dIdPO[iPO] /= POcorr
2632    if 'T' in parmDict[hfx+'Type']:
2633        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2634        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2635    else:
2636        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2637        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2638    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2639       
2640def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2641    'Needs a doc string'
2642    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2643        costh = cosd(refl[5+im]/2.)
2644        #crystallite size
2645        if calcControls[phfx+'SizeType'] == 'isotropic':
2646            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2647        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2648            H = np.array(refl[:3])
2649            P = np.array(calcControls[phfx+'SizeAxis'])
2650            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2651            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2652            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2653        else:           #ellipsoidal crystallites
2654            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2655            H = np.array(refl[:3])
2656            lenR = G2pwd.ellipseSize(H,Sij,GB)
2657            Sgam = 1.8*wave/(np.pi*costh*lenR)
2658        #microstrain               
2659        if calcControls[phfx+'MustrainType'] == 'isotropic':
2660            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2661        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2662            H = np.array(refl[:3])
2663            P = np.array(calcControls[phfx+'MustrainAxis'])
2664            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2665            Si = parmDict[phfx+'Mustrain;i']
2666            Sa = parmDict[phfx+'Mustrain;a']
2667            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2668        else:       #generalized - P.W. Stephens model
2669            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2670            Sum = 0
2671            for i,strm in enumerate(Strms):
2672                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2673            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2674    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2675        #crystallite size
2676        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2677            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2678        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2679            H = np.array(refl[:3])
2680            P = np.array(calcControls[phfx+'SizeAxis'])
2681            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2682            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2683            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2684        else:           #ellipsoidal crystallites   #OK
2685            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2686            H = np.array(refl[:3])
2687            lenR = G2pwd.ellipseSize(H,Sij,GB)
2688            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2689        #microstrain               
2690        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2691            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2692        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2693            H = np.array(refl[:3])
2694            P = np.array(calcControls[phfx+'MustrainAxis'])
2695            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2696            Si = parmDict[phfx+'Mustrain;i']
2697            Sa = parmDict[phfx+'Mustrain;a']
2698            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2699        else:       #generalized - P.W. Stephens model  OK
2700            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2701            Sum = 0
2702            for i,strm in enumerate(Strms):
2703                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2704            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2705           
2706    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2707    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2708    sig /= ateln2
2709    return sig,gam
2710       
2711def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2712    'Needs a doc string'
2713    gamDict = {}
2714    sigDict = {}
2715    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2716        costh = cosd(refl[5+im]/2.)
2717        tanth = tand(refl[5+im]/2.)
2718        #crystallite size derivatives
2719        if calcControls[phfx+'SizeType'] == 'isotropic':
2720            Sgam = 1.8*wave/(np.pi*costh*parmDict[phfx+'Size;i'])
2721            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh*parmDict[phfx+'Size;i']**2)
2722            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2723        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2724            H = np.array(refl[:3])
2725            P = np.array(calcControls[phfx+'SizeAxis'])
2726            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2727            Si = parmDict[phfx+'Size;i']
2728            Sa = parmDict[phfx+'Size;a']
2729            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2730            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2731            Sgam = gami*sqtrm
2732            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2733            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2734            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2735            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2736            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2737            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2738        else:           #ellipsoidal crystallites
2739            const = 1.8*wave/(np.pi*costh)
2740            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2741            H = np.array(refl[:3])
2742            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2743            Sgam = const/lenR
2744            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2745                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2746                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2747        gamDict[phfx+'Size;mx'] = Sgam
2748        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2749               
2750        #microstrain derivatives               
2751        if calcControls[phfx+'MustrainType'] == 'isotropic':
2752            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2753            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2754            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2755        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2756            H = np.array(refl[:3])
2757            P = np.array(calcControls[phfx+'MustrainAxis'])
2758            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2759            Si = parmDict[phfx+'Mustrain;i']
2760            Sa = parmDict[phfx+'Mustrain;a']
2761            gami = 0.018*Si*Sa*tanth/np.pi
2762            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2763            Mgam = gami/sqtrm
2764            dsi = -gami*Si*cosP**2/sqtrm**3
2765            dsa = -gami*Sa*sinP**2/sqtrm**3
2766            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2767            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2768            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2769            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2770        else:       #generalized - P.W. Stephens model
2771            const = 0.018*refl[4+im]**2*tanth/np.pi
2772            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2773            Sum = 0
2774            for i,strm in enumerate(Strms):
2775                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2776                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2777                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2778            Mgam = const*np.sqrt(Sum)
2779            for i in range(len(Strms)):
2780                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2781                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2782        gamDict[phfx+'Mustrain;mx'] = Mgam
2783        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2784    else:   #'T'OF - All checked & OK
2785        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2786            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2787            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2788            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2789        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2790            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2791            H = np.array(refl[:3])
2792            P = np.array(calcControls[phfx+'SizeAxis'])
2793            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2794            Si = parmDict[phfx+'Size;i']
2795            Sa = parmDict[phfx+'Size;a']
2796            gami = const/(Si*Sa)
2797            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2798            Sgam = gami*sqtrm
2799            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2800            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2801            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2802            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2803            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2804            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2805        else:           #OK  ellipsoidal crystallites
2806            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2807            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2808            H = np.array(refl[:3])
2809            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2810            Sgam = const/lenR
2811            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2812                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2813                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2814        gamDict[phfx+'Size;mx'] = Sgam  #OK
2815        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2816               
2817        #microstrain derivatives               
2818        if calcControls[phfx+'MustrainType'] == 'isotropic':
2819            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2820            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2821            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2822        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2823            H = np.array(refl[:3])
2824            P = np.array(calcControls[phfx+'MustrainAxis'])
2825            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2826            Si = parmDict[phfx+'Mustrain;i']
2827            Sa = parmDict[phfx+'Mustrain;a']
2828            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2829            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2830            Mgam = gami/sqtrm
2831            dsi = -gami*Si*cosP**2/sqtrm**3
2832            dsa = -gami*Sa*sinP**2/sqtrm**3
2833            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2834            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2835            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2836            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2837        else:       #generalized - P.W. Stephens model OK
2838            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2839            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2840            Sum = 0
2841            for i,strm in enumerate(Strms):
2842                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2843                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2844                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2845            Mgam = const*np.sqrt(Sum)
2846            for i in range(len(Strms)):
2847                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2848                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2849        gamDict[phfx+'Mustrain;mx'] = Mgam
2850        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2851       
2852    return sigDict,gamDict
2853       
2854def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2855    'Needs a doc string'
2856    if im:
2857        h,k,l,m = refl[:4]
2858        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2859        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2860    else:
2861        h,k,l = refl[:3]
2862        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2863    refl[4+im] = d
2864    if 'C' in calcControls[hfx+'histType']:
2865        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2866        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2867        if 'Bragg' in calcControls[hfx+'instType']:
2868            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2869                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2870        else:               #Debye-Scherrer - simple but maybe not right
2871            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2872    elif 'T' in calcControls[hfx+'histType']:
2873        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2874        #do I need sample position effects - maybe?
2875    return pos
2876
2877def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2878    'Needs a doc string'
2879    dpr = 180./np.pi
2880    if im:
2881        h,k,l,m = refl[:4]
2882        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2883        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2884        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
2885    else:
2886        m = 0
2887        h,k,l = refl[:3]       
2888        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2889    dst = np.sqrt(dstsq)
2890    dsp = 1./dst
2891    if 'C' in calcControls[hfx+'histType']:
2892        pos = refl[5+im]-parmDict[hfx+'Zero']
2893        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2894        dpdw = const*dst
2895        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2896        dpdZ = 1.0
2897        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2898            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2899        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2900        if 'Bragg' in calcControls[hfx+'instType']:
2901            dpdSh = -4.*shft*cosd(pos/2.0)
2902            dpdTr = -shft*sind(pos)*100.0
2903            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2904        else:               #Debye-Scherrer - simple but maybe not right
2905            dpdXd = -shft*cosd(pos)
2906            dpdYd = -shft*sind(pos)
2907            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2908    elif 'T' in calcControls[hfx+'histType']:
2909        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2910        dpdZ = 1.0
2911        dpdDC = dsp
2912        dpdDA = dsp**2
2913        dpdDB = 1./dsp
2914        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2915            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2916        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2917           
2918def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2919    'Needs a doc string'
2920    laue = SGData['SGLaue']
2921    uniq = SGData['SGUniq']
2922    h,k,l = refl[:3]
2923    if laue in ['m3','m3m']:
2924        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2925            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2926    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2927        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2928    elif laue in ['3R','3mR']:
2929        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2930    elif laue in ['4/m','4/mmm']:
2931        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2932    elif laue in ['mmm']:
2933        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2934    elif laue in ['2/m']:
2935        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2936        if uniq == 'a':
2937            Dij += parmDict[phfx+'D23']*k*l
2938        elif uniq == 'b':
2939            Dij += parmDict[phfx+'D13']*h*l
2940        elif uniq == 'c':
2941            Dij += parmDict[phfx+'D12']*h*k
2942    else:
2943        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2944            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2945    if 'C' in calcControls[hfx+'histType']:
2946        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2947    else:
2948        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2949           
2950def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2951    'Needs a doc string'
2952    laue = SGData['SGLaue']
2953    uniq = SGData['SGUniq']
2954    h,k,l = refl[:3]
2955    if laue in ['m3','m3m']:
2956        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2957            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2958    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2959        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2960    elif laue in ['3R','3mR']:
2961        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2962    elif laue in ['4/m','4/mmm']:
2963        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2964    elif laue in ['mmm']:
2965        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2966    elif laue in ['2/m']:
2967        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2968        if uniq == 'a':
2969            dDijDict[phfx+'D23'] = k*l
2970        elif uniq == 'b':
2971            dDijDict[phfx+'D13'] = h*l
2972        elif uniq == 'c':
2973            dDijDict[phfx+'D12'] = h*k
2974    else:
2975        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2976            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2977    if 'C' in calcControls[hfx+'histType']:
2978        for item in dDijDict:
2979            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2980    else:
2981        for item in dDijDict:
2982            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2983    return dDijDict
2984   
2985def GetDij(phfx,SGData,parmDict):
2986    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2987    return G2spc.HStrainVals(HSvals,SGData)
2988               
2989def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2990    '''Compute the observed structure factors for Powder histograms and store in reflection array
2991    Multiprocessing support added
2992    '''
2993    if GSASIIpath.GetConfigValue('Show_timing',False):
2994        starttime = time.time() #; print 'start GetFobsSq'
2995    histoList = list(Histograms.keys())
2996    histoList.sort()
2997    Ka2 = shl = lamRatio = kRatio = None
2998    for histogram in histoList:
2999        if 'PWDR' in histogram[:4]:
3000            Histogram = Histograms[histogram]
3001            hId = Histogram['hId']
3002            hfx = ':%d:'%(hId)
3003            Limits = calcControls[hfx+'Limits']
3004            if 'C' in calcControls[hfx+'histType']:
3005                shl = max(parmDict[hfx+'SH/L'],0.0005)
3006                Ka2 = False
3007                kRatio = 0.0
3008                if hfx+'Lam1' in list(parmDict.keys()):
3009                    Ka2 = True
3010                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3011                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
3012            x,y,w,yc,yb,yd = Histogram['Data']
3013            xMask = ma.getmaskarray(x)
3014            xB = np.searchsorted(x,Limits[0])
3015            xF = np.searchsorted(x,Limits[1])
3016            ymb = np.array(y-yb)
3017            ymb = np.where(ymb,ymb,1.0)
3018            ycmb = np.array(yc-yb)
3019            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3020            refLists = Histogram['Reflection Lists']
3021            for phase in refLists:
3022                if phase not in Phases:     #skips deleted or renamed phases silently!
3023                    continue
3024                Phase = Phases[phase]
3025                if histogram not in Phase['Histograms']:
3026                    continue
3027                im = 0
3028                if Phase['General'].get('Modulated',False):
3029                    im = 1
3030                pId = Phase['pId']
3031                phfx = '%d:%d:'%(pId,hId)
3032                refDict = refLists[phase]
3033                sumFo = 0.0
3034                sumdF = 0.0
3035                sumFosq = 0.0
3036                sumdFsq = 0.0
3037                sumInt = 0.0
3038                nExcl = 0
3039                # test to see if we are using multiprocessing below
3040                useMP,ncores = G2mp.InitMP()
3041                if len(refDict['RefList']) < 100: useMP = False       
3042                if useMP: # multiprocessing: create a set of initialized Python processes
3043                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3044                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3045                    profArgs = [[] for i in range(G2mp.ncores)]
3046                else:
3047                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3048                if 'C' in calcControls[hfx+'histType']:
3049                    # are we multiprocessing?
3050                    for iref,refl in enumerate(refDict['RefList']):
3051                        if useMP: 
3052                            profArgs[iref%G2mp.ncores].append((refl,iref))
3053                        else:
3054                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3055                            if type(icod) is tuple:
3056                                refl[8+im] = icod[0]
3057                                sumInt += icod[1]
3058                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3059                            elif icod == -1:
3060                                refl[3+im] *= -1
3061                                nExcl += 1
3062                            elif icod == -2:
3063                                break
3064                    if useMP:
3065                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3066                            sumInt += sInt
3067                            for refl8im,irefl in resList:
3068                                if refl8im is None:
3069                                    refDict['RefList'][irefl][3+im] *= -1
3070                                    nExcl += 1
3071                                else:
3072                                    refDict['RefList'][irefl][8+im] = refl8im
3073                                    if parmDict[phfx+'LeBail']:
3074                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3075                elif 'T' in calcControls[hfx+'histType']:
3076                    for iref,refl in enumerate(refDict['RefList']):
3077                        if useMP: 
3078                            profArgs[iref%G2mp.ncores].append((refl,iref))
3079                        else:
3080                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3081                            if type(icod) is tuple:
3082                                refl[8+im] = icod[0]
3083                                sumInt += icod[1]
3084                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3085                            elif icod == -1:
3086                                refl[3+im] *= -1
3087                                nExcl += 1
3088                            elif icod == -2:
3089                                break
3090                    if useMP:
3091                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3092                            sumInt += sInt
3093                            for refl8im,irefl in resList:
3094                                if refl8im is None:
3095                                    refDict['RefList'][irefl][3+im] *= -1
3096                                    nExcl += 1
3097                                else:
3098                                    refDict['RefList'][irefl][8+im] = refl8im
3099                                    if parmDict[phfx+'LeBail']:
3100                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3101                if useMP: MPpool.terminate()
3102                sumFo = 0.0
3103                sumdF = 0.0
3104                sumFosq = 0.0
3105                sumdFsq = 0.0
3106                for iref,refl in enumerate(refDict['RefList']):
3107                    Fo = np.sqrt(np.abs(refl[8+im]))
3108                    Fc = np.sqrt(np.abs(refl[9]+im))
3109                    sumFo += Fo
3110                    sumFosq += refl[8+im]**2
3111                    sumdF += np.abs(Fo-Fc)
3112                    sumdFsq += (refl[8+im]-refl[9+im])**2
3113                if sumFo:
3114                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3115                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3116                else:
3117                    Histogram['Residuals'][phfx+'Rf'] = 100.
3118                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3119                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3120                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3121                Histogram['Residuals']['hId'] = hId
3122        elif 'HKLF' in histogram[:4]:
3123            Histogram = Histograms[histogram]
3124            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3125    if GSASIIpath.GetConfigValue('Show_timing',False):
3126        print ('GetFobsSq t=',time.time()-starttime)
3127               
3128def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup,histogram=None):
3129    'Computes the powder pattern for a histogram based on contributions from all used phases'
3130    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3131   
3132    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3133        U = parmDict[hfx+'U']
3134        V = parmDict[hfx+'V']
3135        W = parmDict[hfx+'W']
3136        X = parmDict[hfx+'X']
3137        Y = parmDict[hfx+'Y']
3138        Z = parmDict[hfx+'Z']
3139        tanPos = tand(refl[5+im]/2.0)
3140        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3141        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3142        sig = max(0.001,sig)
3143        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3144        gam = max(0.001,gam)
3145        return sig,gam
3146               
3147    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3148        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3149            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3150        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3151        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3152        sig += Ssig
3153        gam += Sgam
3154        return sig,gam
3155       
3156    def GetReflAlpBet(refl,im,hfx,parmDict):
3157        alp = parmDict[hfx+'alpha']/refl[4+im]
3158        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3159        return alp,bet
3160       
3161    hId = Histogram['hId']
3162    hfx = ':%d:'%(hId)
3163    bakType = calcControls[hfx+'bakType']
3164    fixedBkg = {i:Histogram['Background'][1][i] for i in Histogram['Background'][1] if i.startswith("_")} 
3165    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixedBkg)
3166    yc = np.zeros_like(yb)
3167    cw = np.diff(ma.getdata(x))
3168    cw = np.append(cw,cw[-1])
3169       
3170    if 'C' in calcControls[hfx+'histType']:   
3171        shl = max(parmDict[hfx+'SH/L'],0.002)
3172        Ka2 = False
3173        if hfx+'Lam1' in (parmDict.keys()):
3174            wave = parmDict[hfx+'Lam1']
3175            Ka2 = True
3176            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3177            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3178        else:
3179            wave = parmDict[hfx+'Lam']
3180    else:
3181        shl = 0.
3182    for phase in Histogram['Reflection Lists']:
3183        refDict = Histogram['Reflection Lists'][phase]
3184        if phase not in Phases:     #skips deleted or renamed phases silently!
3185            continue
3186        Phase = Phases[phase]
3187        if histogram and not histogram in Phase['Histograms']:
3188            continue
3189        pId = Phase['pId']
3190        pfx = '%d::'%(pId)
3191        phfx = '%d:%d:'%(pId,hId)
3192        hfx = ':%d:'%(hId)
3193        SGData = Phase['General']['SGData']
3194        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3195        im = 0
3196        if Phase['General'].get('Modulated',False):
3197            SSGData = Phase['General']['SSGData']
3198            im = 1  #offset in SS reflection list
3199            #??
3200        Dij = GetDij(phfx,SGData,parmDict)
3201        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do someting if Dij << 0.
3202        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3203        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3204            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3205        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3206        Vst = np.sqrt(nl.det(G))    #V*
3207        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3208            if im:
3209                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3210            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3211                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3212            else:
3213                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3214        badPeak = False
3215        # test to see if we are using multiprocessing here
3216        useMP,ncores = G2mp.InitMP()
3217        if len(refDict['RefList']) < 100: useMP = False       
3218        if useMP: # multiprocessing: create a set of initialized Python processes
3219            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3220            profArgs = [[] for i in range(ncores)]
3221        if 'C' in calcControls[hfx+'histType']:
3222            for iref,refl in enumerate(refDict['RefList']):
3223                if im:
3224                    h,k,l,m = refl[:4]
3225                else:
3226                    h,k,l = refl[:3]
3227                Uniq = np.inner(refl[:3],SGMT)
3228                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
3229                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3230                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3231                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3232                refl[11+im] *= Vst*Lorenz
3233                 
3234                if Phase['General'].get('doPawley'):
3235                    try:
3236                        if im:
3237                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3238                        else:
3239                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3240                        refl[9+im] = parmDict[pInd]
3241                    except KeyError:
3242#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3243                        continue
3244                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3245                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3246                iFin = np.searchsorted(x,refl[5+im]+fmax)
3247                if not iBeg+iFin:       #peak below low limit - skip peak
3248                    continue
3249                elif not iBeg-iFin:     #peak above high limit - done
3250                    break
3251                elif iBeg > iFin:   #bad peak coeff - skip
3252                    badPeak = True
3253                    continue
3254                if useMP:
3255                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3256                else:
3257                    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
3258                if Ka2:
3259                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3260                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3261                    iBeg = np.searchsorted(x,pos2-fmin)
3262                    iFin = np.searchsorted(x,pos2+fmax)
3263                    if not iBeg+iFin:       #peak below low limit - skip peak
3264                        continue
3265                    elif not iBeg-iFin:     #peak above high limit - done
3266                        return yc,yb
3267                    elif iBeg > iFin:   #bad peak coeff - skip
3268                        continue
3269                    if useMP:
3270                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3271                    else:
3272                        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
3273        elif 'T' in calcControls[hfx+'histType']:
3274            for iref,refl in enumerate(refDict['RefList']):
3275                if im:
3276                    h,k,l,m = refl[:4]
3277                else:
3278                    h,k,l = refl[:3]
3279                Uniq = np.inner(refl[:3],SGMT)
3280                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
3281                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
3282#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
3283                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3284                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3285                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3286                refl[11+im] *= Vst*Lorenz
3287                if Phase['General'].get('doPawley'):
3288                    try:
3289                        if im:
3290                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3291                        else:
3292                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3293                        refl[9+im] = parmDict[pInd]
3294                    except KeyError:
3295#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3296                        continue
3297                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
3298                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3299                iFin = np.searchsorted(x,refl[5+im]+fmax)
3300                if not iBeg+iFin:       #peak below low limit - skip peak
3301                    continue
3302                elif not iBeg-iFin:     #peak above high limit - done
3303                    break
3304                elif iBeg > iFin:   #bad peak coeff - skip
3305                    badPeak = True
3306                    continue
3307                if useMP:
3308                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin))
3309                else:
3310                    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]
3311#        print 'profile calc time: %.3fs'%(time.time()-time0)
3312        if useMP and 'C' in calcControls[hfx+'histType']:
3313            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
3314                yc += y
3315            MPpool.terminate()
3316        elif useMP:
3317            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
3318                yc += y
3319            MPpool.terminate()
3320    if badPeak:
3321        print ('ouch #4 bad profile coefficients yield negative peak width; some reflections skipped')
3322    if GSASIIpath.GetConfigValue('Show_timing',False):
3323        print ('getPowderProfile t=%.3f'%(time.time()-starttime))
3324    return yc,yb
3325   
3326def getPowderProfileDervMP(args):
3327    '''Computes the derivatives of the computed powder pattern with respect to all
3328    refined parameters.
3329    Multiprocessing version.
3330    '''
3331    import pytexture as ptx
3332    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics for each processor
3333    parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars = args[:9]
3334    prc,tprc,histogram = 0,1,None
3335    if len(args) >= 10: prc=args[9]
3336    if len(args) >= 11: tprc=args[10]
3337    if len(args) >= 12: histogram=args[11]
3338    def cellVaryDerv(pfx,SGData,dpdA): 
3339        if SGData['SGLaue'] in ['-1',]:
3340            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
3341                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
3342        elif SGData['SGLaue']