source: trunk/GSASIIstrMath.py @ 4051

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

Correct error in PWDR extinction derivatives

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