source: trunk/GSASIIstrMath.py @ 4377

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

Add end-of-life warning for python 2.7
small fix for mag str. fctr calc

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