source: trunk/GSASIIstrMath.py @ 4514

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

remove non gpx entries from recent projects
continue work on magnetism

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