source: trunk/GSASIIstrMath.py @ 3865

Last change on this file since 3865 was 3865, checked in by vondreele, 4 years ago

a new attempt at incommensurate mag. str. fctrs.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 217.9 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMath - structure math routines*
4-----------------------------------------
5'''
6########### SVN repository information ###################
7# $Date: 2019-03-31 23:59:15 +0000 (Sun, 31 Mar 2019) $
8# $Author: vondreele $
9# $Revision: 3865 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 3865 2019-03-31 23:59:15Z 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: 3865 $")
23import GSASIIElem as G2el
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import GSASIIpwd as G2pwd
27import GSASIImapvars as G2mv
28import GSASIImath as G2mth
29import GSASIIobj as G2obj
30import GSASIImpsubs as G2mp
31#G2mp.InitMP(False)  # This disables multiprocessing
32
33sind = lambda x: np.sin(x*np.pi/180.)
34cosd = lambda x: np.cos(x*np.pi/180.)
35tand = lambda x: np.tan(x*np.pi/180.)
36asind = lambda x: 180.*np.arcsin(x)/np.pi
37acosd = lambda x: 180.*np.arccos(x)/np.pi
38atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
39   
40ateln2 = 8.0*np.log(2.0)
41twopi = 2.0*np.pi
42twopisq = 2.0*np.pi**2
43nxs = np.newaxis
44
45################################################################################
46##### Rigid Body Models
47################################################################################
48       
49def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
50    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
51    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
52    '''
53    atxIds = ['Ax:','Ay:','Az:']
54    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
55    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
56    if not RBIds['Vector'] and not RBIds['Residue']:
57        return
58    VRBIds = RBIds['Vector']
59    RRBIds = RBIds['Residue']
60    if Update:
61        RBData = rigidbodyDict
62    else:
63        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
64    if RBIds['Vector']:                       # first update the vector magnitudes
65        VRBData = RBData['Vector']
66        for i,rbId in enumerate(VRBIds):
67            if VRBData[rbId]['useCount']:
68                for j in range(len(VRBData[rbId]['VectMag'])):
69                    name = '::RBV;'+str(j)+':'+str(i)
70                    VRBData[rbId]['VectMag'][j] = parmDict[name]
71    for phase in Phases:
72        Phase = Phases[phase]
73        General = Phase['General']
74        cx,ct,cs,cia = General['AtomPtrs']
75        cell = General['Cell'][1:7]
76        Amat,Bmat = G2lat.cell2AB(cell)
77        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
78        pfx = str(Phase['pId'])+'::'
79        if Update:
80            RBModels = Phase['RBModels']
81        else:
82            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
83        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
84            jrb = VRBIds.index(RBObj['RBId'])
85            rbsx = str(irb)+':'+str(jrb)
86            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
87                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
88            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
89                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
90            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
91            TLS = RBObj['ThermalMotion']
92            if 'T' in TLS[0]:
93                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
94                    TLS[1][i] = parmDict[pfx+pt+rbsx]
95            if 'L' in TLS[0]:
96                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
97                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
98            if 'S' in TLS[0]:
99                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
100                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
101            if 'U' in TLS[0]:
102                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
103            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
104            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
105            for i,x in enumerate(XYZ):
106                atId = RBObj['Ids'][i]
107                for j in [0,1,2]:
108                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
109                if UIJ[i][0] == 'A':
110                    for j in range(6):
111                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
112                elif UIJ[i][0] == 'I':
113                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
114           
115        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
116            jrb = RRBIds.index(RBObj['RBId'])
117            rbsx = str(irb)+':'+str(jrb)
118            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
119                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
120            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
121                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
122            RBObj['Orient'][0] = G2mth.normQ(RBObj['Orient'][0])
123            TLS = RBObj['ThermalMotion']
124            if 'T' in TLS[0]:
125                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
126                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
127            if 'L' in TLS[0]:
128                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
129                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
130            if 'S' in TLS[0]:
131                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
132                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
133            if 'U' in TLS[0]:
134                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
135            for itors,tors in enumerate(RBObj['Torsions']):
136                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
137            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
138            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
139            for i,x in enumerate(XYZ):
140                atId = RBObj['Ids'][i]
141                for j in [0,1,2]:
142                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
143                if UIJ[i][0] == 'A':
144                    for j in range(6):
145                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
146                elif UIJ[i][0] == 'I':
147                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
148                   
149def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
150    'Needs a doc string'
151    atxIds = ['dAx:','dAy:','dAz:']
152    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
153    OIds = ['Oa:','Oi:','Oj:','Ok:']
154    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
155    if not RBIds['Vector'] and not RBIds['Residue']:
156        return
157    VRBIds = RBIds['Vector']
158    RRBIds = RBIds['Residue']
159    RBData = rigidbodyDict
160    for item in parmDict:
161        if 'RB' in item:
162            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
163    General = Phase['General']
164    cx,ct,cs,cia = General['AtomPtrs']
165    cell = General['Cell'][1:7]
166    Amat,Bmat = G2lat.cell2AB(cell)
167    rpd = np.pi/180.
168    rpd2 = rpd**2
169    g = nl.inv(np.inner(Bmat,Bmat))
170    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
171        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
172    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'],cia+8)
173    pfx = str(Phase['pId'])+'::'
174    RBModels =  Phase['RBModels']
175    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
176        VModel = RBData['Vector'][RBObj['RBId']]
177        Q = RBObj['Orient'][0]
178        jrb = VRBIds.index(RBObj['RBId'])
179        rbsx = str(irb)+':'+str(jrb)
180        dXdv = []
181        for iv in range(len(VModel['VectMag'])):
182            dCdv = []
183            for vec in VModel['rbVect'][iv]:
184                dCdv.append(G2mth.prodQVQ(Q,vec))
185            dXdv.append(np.inner(Bmat,np.array(dCdv)).T)
186        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
187        for ia,atId in enumerate(RBObj['Ids']):
188            atNum = AtLookup[atId]
189            dx = 0.00001
190            for iv in range(len(VModel['VectMag'])):
191                for ix in [0,1,2]:
192                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
193            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
194                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
195            for iv in range(4):
196                Q[iv] -= dx
197                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
198                Q[iv] += 2.*dx
199                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
200                Q[iv] -= dx
201                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
202                for ix in [0,1,2]:
203                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
204            X = G2mth.prodQVQ(Q,Cart[ia])
205            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
206            dFdu = G2lat.U6toUij(dFdu.T)
207            dFdu = np.tensordot(Amat,np.tensordot(Amat,dFdu,([1,0])),([0,1]))           
208            dFdu = G2lat.UijtoU6(dFdu)
209            atNum = AtLookup[atId]
210            if 'T' in RBObj['ThermalMotion'][0]:
211                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
212                    dFdvDict[pfx+name+rbsx] += dFdu[i]
213            if 'L' in RBObj['ThermalMotion'][0]:
214                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
215                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
216                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
217                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
218                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
219                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(-dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
220                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
221                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(-dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
222                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
223            if 'S' in RBObj['ThermalMotion'][0]:
224                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
225                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
226                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
227                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
228                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
229                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
230                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
231                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
232            if 'U' in RBObj['ThermalMotion'][0]:
233                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
234
235
236    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
237        Q = RBObj['Orient'][0]
238        jrb = RRBIds.index(RBObj['RBId'])
239        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
240        rbsx = str(irb)+':'+str(jrb)
241        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
242        for itors,tors in enumerate(RBObj['Torsions']):     #derivative error?
243            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
244            orId,pvId = torData[itors][:2]
245            pivotVec = Cart[orId]-Cart[pvId]
246            QA = G2mth.AVdeg2Q(-0.001,pivotVec)
247            QB = G2mth.AVdeg2Q(0.001,pivotVec)
248            for ir in torData[itors][3]:
249                atNum = AtLookup[RBObj['Ids'][ir]]
250                rVec = Cart[ir]-Cart[pvId]
251                dR = G2mth.prodQVQ(QB,rVec)-G2mth.prodQVQ(QA,rVec)
252                dRdT = np.inner(Bmat,G2mth.prodQVQ(Q,dR))/.002
253                for ix in [0,1,2]:
254                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
255        for ia,atId in enumerate(RBObj['Ids']):
256            atNum = AtLookup[atId]
257            dx = 0.00001
258            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
259                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
260            for iv in range(4):
261                Q[iv] -= dx
262                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
263                Q[iv] += 2.*dx
264                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
265                Q[iv] -= dx
266                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
267                for ix in [0,1,2]:
268                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
269            X = G2mth.prodQVQ(Q,Cart[ia])
270            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T/gvec
271            dFdu = G2lat.U6toUij(dFdu.T)
272            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
273            dFdu = G2lat.UijtoU6(dFdu)
274            atNum = AtLookup[atId]
275            if 'T' in RBObj['ThermalMotion'][0]:
276                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
277                    dFdvDict[pfx+name+rbsx] += dFdu[i]
278            if 'L' in RBObj['ThermalMotion'][0]:
279                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
280                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
281                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
282                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
283                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
284                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[4]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
285                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
286                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[5]*X[0]**2-2.*dFdu[0]*X[1]*X[2]+
287                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
288            if 'S' in RBObj['ThermalMotion'][0]:
289                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
290                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
291                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
292                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
293                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
294                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
295                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
296                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
297            if 'U' in RBObj['ThermalMotion'][0]:
298                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
299   
300################################################################################
301##### Penalty & restraint functions
302################################################################################
303
304def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
305    'Compute user-supplied and built-in restraint functions'
306    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
307    pNames = []
308    pVals = []
309    pWt = []
310    negWt = {}
311    pWsum = {}
312    pWnum = {}
313    for phase in Phases:
314        pId = Phases[phase]['pId']
315        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
316        General = Phases[phase]['General']
317        cx,ct,cs,cia = General['AtomPtrs']
318        textureData = General['SH Texture']
319        SGData = General['SGData']
320        Atoms = Phases[phase]['Atoms']
321        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
322        cell = General['Cell'][1:7]
323        Amat,Bmat = G2lat.cell2AB(cell)
324        if phase not in restraintDict:
325            continue
326        phaseRest = restraintDict[phase]
327        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
328            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
329            ['ChemComp','Sites'],['Texture','HKLs'],['General','General'],]
330        for name,rest in names:
331            pWsum[name] = 0.
332            pWnum[name] = 0
333            if name not in phaseRest:
334                continue
335            itemRest = phaseRest[name]
336            if itemRest[rest] and itemRest['Use']:
337                wt = itemRest.get('wtFactor',1.)
338                if name in ['Bond','Angle','Plane','Chiral']:
339                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
340                        pNames.append(str(pId)+':'+name+':'+str(i))
341                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
342                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
343                        if name == 'Bond':
344                            calc = G2mth.getRestDist(XYZ,Amat)
345                        elif name == 'Angle':
346                            calc = G2mth.getRestAngle(XYZ,Amat)
347                        elif name == 'Plane':
348                            calc = G2mth.getRestPlane(XYZ,Amat)
349                        elif name == 'Chiral':
350                            calc = G2mth.getRestChiral(XYZ,Amat)
351                        pVals.append(obs-calc)
352                        pWt.append(wt/esd**2)
353                        pWsum[name] += wt*((obs-calc)/esd)**2
354                        pWnum[name] += 1
355                elif name in ['Torsion','Rama']:
356                    coeffDict = itemRest['Coeff']
357                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
358                        pNames.append(str(pId)+':'+name+':'+str(i))
359                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
360                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
361                        if name == 'Torsion':
362                            tor = G2mth.getRestTorsion(XYZ,Amat)
363                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
364                        else:
365                            phi,psi = G2mth.getRestRama(XYZ,Amat)
366                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
367                        pVals.append(restr)
368                        pWt.append(wt/esd**2)
369                        pWsum[name] += wt*(restr/esd)**2
370                        pWnum[name] += 1
371                elif name == 'ChemComp':
372                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
373                        pNames.append(str(pId)+':'+name+':'+str(i))
374                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
375                        frac = np.array(G2mth.GetAtomFracByID(pId,parmDict,AtLookup,indx))
376                        calc = np.sum(mul*frac*factors)
377                        pVals.append(obs-calc)
378                        pWt.append(wt/esd**2)                   
379                        pWsum[name] += wt*((obs-calc)/esd)**2
380                        pWnum[name] += 1
381                elif name == 'Texture':
382                    SHkeys = list(textureData['SH Coeff'][1].keys())
383                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
384                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
385                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
386                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
387                        PH = np.array(hkl)
388                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
389                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
390                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
391                        Z1 = ma.masked_greater(Z,0.0)           #is this + or -?
392                        IndZ1 = np.array(ma.nonzero(Z1))
393                        for ind in IndZ1.T:
394                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
395                            pVals.append(Z1[ind[0]][ind[1]])
396                            pWt.append(wt/esd1**2)
397                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
398                            pWnum[name] += 1
399                        if ifesd2:
400                            Z2 = 1.-Z
401                            for ind in np.ndindex(grid,grid):
402                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
403                                pVals.append(Z2[ind[0]][ind[1]])
404                                pWt.append(wt/esd2**2)
405                                pWsum[name] += wt*(Z2/esd2)**2
406                                pWnum[name] += 1
407                elif name == 'General':
408                    for i,(eq,obs,esd) in enumerate(itemRest[rest]):
409                        pNames.append(str(pId)+':'+name+':'+str(i))
410                        calcobj = G2obj.ExpressionCalcObj(eq)
411                        calcobj.SetupCalc(parmDict)
412                        calc = calcobj.EvalExpression()
413                        pVals.append(obs-calc)
414                        pWt.append(wt/esd**2)                   
415                        pWsum[name] += wt*((obs-calc)/esd)**2
416                        pWnum[name] += 1
417       
418    for phase in Phases:
419        name = 'SH-Pref.Ori.'
420        pId = Phases[phase]['pId']
421        General = Phases[phase]['General']
422        SGData = General['SGData']
423        cell = General['Cell'][1:7]
424        pWsum[name] = 0.0
425        pWnum[name] = 0
426        for hist in Phases[phase]['Histograms']:
427            if not Phases[phase]['Histograms'][hist]['Use']:
428                continue
429            if hist in Histograms and 'PWDR' in hist:
430                hId = Histograms[hist]['hId']
431                phfx = '%d:%d:'%(pId,hId)
432                if calcControls[phfx+'poType'] == 'SH':
433                    toler = calcControls[phfx+'SHtoler']
434                    wt = 1./toler**2
435                    HKLs = np.array(calcControls[phfx+'SHhkl'])
436                    SHnames = calcControls[phfx+'SHnames']
437                    SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
438                    for i,PH in enumerate(HKLs):
439                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
440                        SH3Coef = {}
441                        for item in SHcof:
442                            L,N = eval(item.strip('C'))
443                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
444                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
445                        X = np.linspace(0,90.0,26)
446                        Y = ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)       #+ or -?
447                        IndY = ma.nonzero(Y)
448                        for ind in IndY[0]:
449                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
450                            pVals.append(Y[ind])
451                            pWt.append(wt)
452                            pWsum[name] += wt*(Y[ind])**2
453                            pWnum[name] += 1
454    pWsum['PWLref'] = 0.
455    pWnum['PWLref'] = 0
456    for item in varyList:
457        if 'PWLref' in item and parmDict[item] < 0.:
458            pId = int(item.split(':')[0])
459            if negWt[pId]:
460                pNames.append(item)
461                pVals.append(parmDict[item])
462                pWt.append(negWt[pId])
463                pWsum['PWLref'] += negWt[pId]*(parmDict[item])**2
464                pWnum['PWLref'] += 1
465    pVals = np.array(pVals)
466    pWt = np.array(pWt)         #should this be np.sqrt?
467    return pNames,pVals,pWt,pWsum,pWnum
468   
469def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
470    '''Compute derivatives on user-supplied and built-in restraint
471    (penalty) functions
472
473    where pNames is list of restraint labels
474
475    returns pDerv with partial derivatives by variable# in varList and
476       restraint# in pNames (pDerv[variable#][restraint#])
477    '''
478    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
479    pDerv = np.zeros((len(varyList),len(pVal)))
480    for pName in pNames: # loop over restraints
481        if 'General' == pName.split(':')[1]:
482            # initialize for General restraint(s) here
483            GeneralInit = True
484            parmDict0 = parmDict.copy()
485            # setup steps for each parameter
486            stepDict = {}
487            for parm in varyList:
488                stepDict[parm] = G2obj.getVarStep(parm,parmDict)
489            break
490    for phase in Phases:
491#        if phase not in restraintDict:
492#            continue
493        pId = Phases[phase]['pId']
494        General = Phases[phase]['General']
495        cx,ct,cs,cia = General['AtomPtrs']
496        SGData = General['SGData']
497        Atoms = Phases[phase]['Atoms']
498        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'],cia+8)
499        cell = General['Cell'][1:7]
500        Amat,Bmat = G2lat.cell2AB(cell)
501        textureData = General['SH Texture']
502
503        SHkeys = list(textureData['SH Coeff'][1].keys())
504        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
505        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
506        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
507        sam = SamSym[textureData['Model']]
508        phaseRest = restraintDict.get(phase,{})
509        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
510            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
511            'ChemComp':'Sites','Texture':'HKLs'}
512        lasthkl = np.array([0,0,0])
513        for ip,pName in enumerate(pNames): # loop over restraints
514            pnames = pName.split(':')
515            if pId == int(pnames[0]):
516                name = pnames[1]
517                if 'PWL' in pName:
518                    pDerv[varyList.index(pName)][ip] += 1.
519                    continue
520                elif 'SH-' in pName:
521                    continue
522                Id = int(pnames[2]) 
523                itemRest = phaseRest[name]
524                if name in ['Bond','Angle','Plane','Chiral']:
525                    indx,ops,obs,esd = itemRest[names[name]][Id]
526                    dNames = []
527                    for ind in indx:
528                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
529                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
530                    if name == 'Bond':
531                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
532                    elif name == 'Angle':
533                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
534                    elif name == 'Plane':
535                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
536                    elif name == 'Chiral':
537                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
538                elif name in ['Torsion','Rama']:
539                    coffDict = itemRest['Coeff']
540                    indx,ops,cofName,esd = itemRest[names[name]][Id]
541                    dNames = []
542                    for ind in indx:
543                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
544                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
545                    if name == 'Torsion':
546                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
547                    else:
548                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
549                elif name == 'ChemComp':
550                    indx,factors,obs,esd = itemRest[names[name]][Id]
551                    dNames = []
552                    for ind in indx:
553                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
554                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
555                        deriv = mul*factors
556                elif 'Texture' in name:
557                    deriv = []
558                    dNames = []
559                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][Id]
560                    hkl = np.array(hkl)
561                    if np.any(lasthkl-hkl):
562                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
563                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
564                        lasthkl = copy.copy(hkl)                       
565                    if 'unit' in name:
566                        pass
567                    else:
568                        gam = float(pnames[3])
569                        psi = float(pnames[4])
570                        for SHname in ODFln:
571                            l,m,n = eval(SHname[1:])
572                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
573                            dNames += [str(pId)+'::'+SHname]
574                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
575                elif name == 'General':
576                    deriv = []
577                    dNames = []
578                    eq,obs,esd = itemRest[name][Id]
579                    calcobj = G2obj.ExpressionCalcObj(eq)
580                    parmlist = list(eq.assgnVars.values()) # parameters used in this expression
581                    for parm in parmlist: # expand list if any parms are determined by constraints
582                        if parm in G2mv.dependentVars:
583                            parmlist += G2mv.independentVars
584                            break
585                    for ind,var in enumerate(varyList):
586                        drv = 0
587                        if var in parmlist:
588                            step = stepDict.get(var,1e-5)
589                            calc = []
590                            # apply step to parameter
591                            oneparm = True
592                            for s in -step,2*step:
593                                parmDict[var] += s
594                                # extend shift if needed to other parameters
595                                if var in G2mv.independentVars:
596                                    G2mv.Dict2Map(parmDict,[])
597                                    oneparm = False
598                                elif var in G2mv.dependentVars:
599                                    G2mv.Map2Dict(parmDict,[])
600                                    oneparm = False
601                                if 'RB' in var:
602                                    ApplyRBModels(parmDict,Phases,rigidbodyDict)
603# test
604                                    oneparm = False
605                                calcobj.SetupCalc(parmDict)
606                                calc.append(calcobj.EvalExpression())
607                            drv = (calc[1]-calc[0])*.5/step
608                            # restore the dict
609                            if oneparm:
610                                parmDict[var] = parmDict0[var]
611                            else:
612                                parmDict = parmDict0.copy()
613                        else:
614                            drv = 0
615                        pDerv[ind][ip] = drv
616                # Add derivatives into matrix, if needed
617                for dName,drv in zip(dNames,deriv):
618                    try:
619                        ind = varyList.index(dName)
620                        pDerv[ind][ip] += drv
621                    except ValueError:
622                        pass
623       
624        lasthkl = np.array([0,0,0])
625        for ip,pName in enumerate(pNames):
626            deriv = []
627            dNames = []
628            pnames = pName.split(':')
629            if 'SH-' in pName and pId == int(pnames[0]):
630                hId = int(pnames[1])
631                phfx = '%d:%d:'%(pId,hId)
632                psi = float(pnames[4])
633                HKLs = calcControls[phfx+'SHhkl']
634                SHnames = calcControls[phfx+'SHnames']
635                SHcof = dict(zip(SHnames,[parmDict[phfx+cof] for cof in SHnames]))
636                hkl = np.array(HKLs[int(pnames[3])])     
637                if np.any(lasthkl-hkl):
638                    phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
639                    SH3Coef = {}
640                    for item in SHcof:
641                        L,N = eval(item.strip('C'))
642                        SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
643                    ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
644                    lasthkl = copy.copy(hkl)                       
645                for SHname in SHnames:
646                    l,n = eval(SHname[1:])
647                    SH3name = 'C%d,0,%d'%(l,n)
648                    Ksl = G2lat.GetKsl(l,0,'0',psi,0.0)[0]
649                    dNames += [phfx+SHname]
650                    deriv.append(ODFln[SH3name][0]*Ksl/SHcof[SHname])
651            for dName,drv in zip(dNames,deriv):
652                try:
653                    ind = varyList.index(dName)
654                    pDerv[ind][ip] += drv
655                except ValueError:
656                    pass
657    return pDerv
658
659################################################################################
660##### Function & derivative calculations
661################################################################################       
662                   
663def GetAtomFXU(pfx,calcControls,parmDict):
664    'Needs a doc string'
665    Natoms = calcControls['Natoms'][pfx]
666    Tdata = Natoms*[' ',]
667    Mdata = np.zeros(Natoms)
668    IAdata = Natoms*[' ',]
669    Fdata = np.zeros(Natoms)
670    Xdata = np.zeros((3,Natoms))
671    dXdata = np.zeros((3,Natoms))
672    Uisodata = np.zeros(Natoms)
673    Uijdata = np.zeros((6,Natoms))
674    Gdata = np.zeros((3,Natoms))
675    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
676        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
677        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
678        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
679        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5],
680        'AMx:':Gdata[0],'AMy:':Gdata[1],'AMz:':Gdata[2],}
681    for iatm in range(Natoms):
682        for key in keys:
683            parm = pfx+key+str(iatm)
684            if parm in parmDict:
685                keys[key][iatm] = parmDict[parm]
686    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
687    Gdata = np.where(Gdata,Gdata,1.e-8)         #avoid divide by zero in derivative calc.
688   
689    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata
690   
691def GetAtomSSFXU(pfx,calcControls,parmDict):
692    'Needs a doc string'
693    Natoms = calcControls['Natoms'][pfx]
694    maxSSwave = calcControls['maxSSwave'][pfx]
695    Nwave = {'F':maxSSwave['Sfrac'],'X':maxSSwave['Spos'],'Y':maxSSwave['Spos'],'Z':maxSSwave['Spos'],
696        'U':maxSSwave['Sadp'],'M':maxSSwave['Smag'],'T':maxSSwave['Spos']}
697    XSSdata = np.zeros((6,maxSSwave['Spos'],Natoms))
698    FSSdata = np.zeros((2,maxSSwave['Sfrac'],Natoms))
699    USSdata = np.zeros((12,maxSSwave['Sadp'],Natoms))
700    MSSdata = np.zeros((6,maxSSwave['Smag'],Natoms))
701    waveTypes = []
702    keys = {'Fsin:':FSSdata[0],'Fcos:':FSSdata[1],'Fzero:':FSSdata[0],'Fwid:':FSSdata[1],
703        'Tmin:':XSSdata[0],'Tmax:':XSSdata[1],'Xmax:':XSSdata[2],'Ymax:':XSSdata[3],'Zmax:':XSSdata[4],
704        'Xsin:':XSSdata[0],'Ysin:':XSSdata[1],'Zsin:':XSSdata[2],'Xcos:':XSSdata[3],'Ycos:':XSSdata[4],'Zcos:':XSSdata[5],
705        'U11sin:':USSdata[0],'U22sin:':USSdata[1],'U33sin:':USSdata[2],'U12sin:':USSdata[3],'U13sin:':USSdata[4],'U23sin:':USSdata[5],
706        'U11cos:':USSdata[6],'U22cos:':USSdata[7],'U33cos:':USSdata[8],'U12cos:':USSdata[9],'U13cos:':USSdata[10],'U23cos:':USSdata[11],
707        'MXsin:':MSSdata[0],'MYsin:':MSSdata[1],'MZsin:':MSSdata[2],'MXcos:':MSSdata[3],'MYcos:':MSSdata[4],'MZcos:':MSSdata[5]}
708    for iatm in range(Natoms):
709        wavetype = [parmDict.get(pfx+kind+'waveType:'+str(iatm),'') for kind in ['F','P','A','M']]
710        waveTypes.append(wavetype)
711        for key in keys:
712            for m in range(Nwave[key[0]]):
713                parm = pfx+key+str(iatm)+':%d'%(m)
714                if parm in parmDict:
715                    keys[key][m][iatm] = parmDict[parm]
716    return waveTypes,FSSdata,XSSdata,USSdata,MSSdata
717   
718def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
719    ''' Compute structure factors for all h,k,l for phase
720    puts the result, F^2, in each ref[8] in refList
721    operates on blocks of 100 reflections for speed
722    input:
723   
724    :param dict refDict: where
725        'RefList' list where each ref = h,k,l,it,d,...
726        'FF' dict of form factors - filed in below
727    :param np.array G:      reciprocal metric tensor
728    :param str pfx:    phase id string
729    :param dict SGData: space group info. dictionary output from SpcGroup
730    :param dict calcControls:
731    :param dict ParmDict:
732
733    '''       
734    phfx = pfx.split(':')[0]+hfx
735    ast = np.sqrt(np.diag(G))
736    Mast = twopisq*np.multiply.outer(ast,ast)
737    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
738    SGT = np.array([ops[1] for ops in SGData['SGOps']])
739    FFtables = calcControls['FFtables']
740    BLtables = calcControls['BLtables']
741    Amat,Bmat = G2lat.Gmat2AB(G)
742    Flack = 1.0
743    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
744        Flack = 1.-2.*parmDict[phfx+'Flack']
745    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
746    TwDict = refDict.get('TwDict',{})           
747    if 'S' in calcControls[hfx+'histType']:
748        NTL = calcControls[phfx+'NTL']
749        NM = calcControls[phfx+'TwinNMN']+1
750        TwinLaw = calcControls[phfx+'TwinLaw']
751        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
752        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
753    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
754        GetAtomFXU(pfx,calcControls,parmDict)
755    if not Xdata.size:          #no atoms in phase!
756        return
757    if 'NC' in calcControls[hfx+'histType']:
758        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
759    elif 'X' in calcControls[hfx+'histType']:
760        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
761        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
762    Uij = np.array(G2lat.U6toUij(Uijdata))
763    bij = Mast*Uij.T
764    blkSize = 100       #no. of reflections in a block - size seems optimal
765    nRef = refDict['RefList'].shape[0]
766    SQ = 1./(2.*refDict['RefList'].T[4])**2
767    if 'N' in calcControls[hfx+'histType']:
768        dat = G2el.getBLvalues(BLtables)
769        refDict['FF']['El'] = list(dat.keys())
770        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
771    else:       #'X'
772        dat = G2el.getFFvalues(FFtables,0.)
773        refDict['FF']['El'] = list(dat.keys())
774        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
775        for iel,El in enumerate(refDict['FF']['El']):
776            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
777#reflection processing begins here - big arrays!
778    iBeg = 0
779    while iBeg < nRef:
780        iFin = min(iBeg+blkSize,nRef)
781        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
782        H = refl.T[:3]                          #array(blkSize,3)
783        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
784        TwMask = np.any(H,axis=-1)
785        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
786            for ir in range(blkSize):
787                iref = ir+iBeg
788                if iref in TwDict:
789                    for i in TwDict[iref]:
790                        for n in range(NTL):
791                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
792            TwMask = np.any(H,axis=-1)
793        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
794        SQfactor = 4.0*SQ*twopisq               #ditto prev.
795        if 'T' in calcControls[hfx+'histType']:
796            if 'P' in calcControls[hfx+'histType']:
797                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
798            else:
799                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
800            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
801            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
802        Uniq = np.inner(H,SGMT)
803        Phi = np.inner(H,SGT)
804        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
805        sinp = np.sin(phase)
806        cosp = np.cos(phase)
807        biso = -SQfactor*Uisodata[:,nxs]
808        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
809        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
810        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
811        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
812        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
813        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
814        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
815        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
816            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
817            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
818        else:
819            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
820            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
821        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
822        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
823        if SGData['SGInv']: #centrosymmetric; B=0
824            fbs[0] *= 0.
825            fas[1] *= 0.
826        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
827            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
828            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
829        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
830            if len(TwinLaw) > 1:
831                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
832                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
833                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
834                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
835            else:   # checked correct!!
836                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
837                refl.T[7] = np.copy(refl.T[9])               
838                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
839#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
840        iBeg += blkSize
841#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
842   
843def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
844    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
845    faster than StructureFactorDerv - correct for powders/nontwins!!
846    input:
847   
848    :param dict refDict: where
849        'RefList' list where each ref = h,k,l,it,d,...
850        'FF' dict of form factors - filled in below
851    :param np.array G:      reciprocal metric tensor
852    :param str hfx:    histogram id string
853    :param str pfx:    phase id string
854    :param dict SGData: space group info. dictionary output from SpcGroup
855    :param dict calcControls:
856    :param dict parmDict:
857   
858    :returns: dict dFdvDict: dictionary of derivatives
859    '''
860    phfx = pfx.split(':')[0]+hfx
861    ast = np.sqrt(np.diag(G))
862    Mast = twopisq*np.multiply.outer(ast,ast)
863    SGMT = np.array([ops[0] for ops in SGData['SGOps']])
864    SGT = np.array([ops[1] for ops in SGData['SGOps']])
865    FFtables = calcControls['FFtables']
866    BLtables = calcControls['BLtables']
867    Amat,Bmat = G2lat.Gmat2AB(G)
868    nRef = len(refDict['RefList'])
869    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
870        GetAtomFXU(pfx,calcControls,parmDict)
871    if not Xdata.size:          #no atoms in phase!
872        return {}
873    mSize = len(Mdata)
874    FF = np.zeros(len(Tdata))
875    if 'NC' in calcControls[hfx+'histType']:
876        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
877    elif 'X' in calcControls[hfx+'histType']:
878        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
879        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
880    Uij = np.array(G2lat.U6toUij(Uijdata))
881    bij = Mast*Uij.T
882    dFdvDict = {}
883    dFdfr = np.zeros((nRef,mSize))
884    dFdx = np.zeros((nRef,mSize,3))
885    dFdui = np.zeros((nRef,mSize))
886    dFdua = np.zeros((nRef,mSize,6))
887    dFdbab = np.zeros((nRef,2))
888    dFdfl = np.zeros((nRef))
889    Flack = 1.0
890    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
891        Flack = 1.-2.*parmDict[phfx+'Flack']
892#    time0 = time.time()
893#reflection processing begins here - big arrays!
894    iBeg = 0
895    blkSize = 32       #no. of reflections in a block - optimized for speed
896    while iBeg < nRef:
897        iFin = min(iBeg+blkSize,nRef)
898        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
899        H = refl.T[:3].T
900        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
901        SQfactor = 8.0*SQ*np.pi**2
902        if 'T' in calcControls[hfx+'histType']:
903            if 'P' in calcControls[hfx+'histType']:
904                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
905            else:
906                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
907            FP = np.repeat(FP.T,len(SGT),axis=0)
908            FPP = np.repeat(FPP.T,len(SGT),axis=0)
909        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
910        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
911        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
912        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
913        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
914        Phi = np.inner(H,SGT)
915        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
916        sinp = np.sin(phase)        #refBlk x nOps x nAtoms
917        cosp = np.cos(phase)
918        occ = Mdata*Fdata/len(SGT)
919        biso = -SQfactor*Uisodata[:,nxs]
920        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
921        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
922        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
923        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
924        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
925        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
926        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
927        if len(FPP.shape) > 1:
928            fotp = np.reshape(FPP,cosp.shape)*Tcorr
929        else:
930            fotp = FPP*Tcorr     
931        if 'T' in calcControls[hfx+'histType']:
932            fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
933            fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
934        else:
935            fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr])
936            fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr])
937        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,refBlk,nTwins)
938        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
939        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,refBlk,nEqv,nAtoms)
940        fbx = np.array([fot*cosp,-fotp*sinp])
941        #sum below is over Uniq
942        dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
943        dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
944        dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
945        dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
946        dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
947        # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6)
948        if not SGData['SGInv']:
949            dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
950            dfbdba = np.sum(-sinp*Tcorr,axis=-2)
951            dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
952            dfbdfl = np.sum(np.sum(fotp*cosp,axis=-1),axis=-1)
953            dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2)           
954            dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2)
955            dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2)
956        else:
957            dfbdfr = np.zeros_like(dfadfr)
958            dfbdx = np.zeros_like(dfadx)
959            dfbdui = np.zeros_like(dfadui)
960            dfbdua = np.zeros_like(dfadua)
961            dfbdba = np.zeros_like(dfadba)
962            dfadfl = 0.0
963            dfbdfl = 0.0
964        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
965        SA = fas[0]+fas[1]
966        SB = fbs[0]+fbs[1]
967        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
968            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
969            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
970            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
971            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
972        else:
973            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
974            dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
975            dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
976            dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
977            dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
978        dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \
979                            fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
980        iBeg += blkSize
981#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
982        #loop over atoms - each dict entry is list of derivatives for all the reflections
983    for i in range(len(Mdata)):
984        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
985        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
986        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
987        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
988        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
989        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
990        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
991        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
992        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
993        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
994        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
995    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
996    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
997    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
998    return dFdvDict
999   
1000def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1001    ''' Compute neutron magnetic structure factors for all h,k,l for phase
1002    puts the result, F^2, in each ref[8] in refList
1003    operates on blocks of 100 reflections for speed
1004    input:
1005   
1006    :param dict refDict: where
1007        'RefList' list where each ref = h,k,l,it,d,...
1008        'FF' dict of form factors - filed in below
1009    :param np.array G:      reciprocal metric tensor
1010    :param str pfx:    phase id string
1011    :param dict SGData: space group info. dictionary output from SpcGroup
1012    :param dict calcControls:
1013    :param dict ParmDict:
1014       
1015    :returns: copy of new refList - used in calculating numerical derivatives
1016
1017    '''       
1018    g = nl.inv(G)
1019    ast = np.sqrt(np.diag(G))
1020    ainv = np.sqrt(np.diag(g))
1021    GS = G/np.outer(ast,ast)
1022    Ginv = g/np.outer(ainv,ainv)
1023    uAmat = G2lat.Gmat2AB(GS)[0]
1024    Mast = twopisq*np.multiply.outer(ast,ast)
1025    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1026    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1027    Ncen = len(SGData['SGCen'])
1028    Nops = len(SGMT)*Ncen
1029    if not SGData['SGFixed']:
1030        Nops *= (1+SGData['SGInv'])
1031    MFtables = calcControls['MFtables']
1032    Bmat = G2lat.Gmat2AB(G)[1]
1033    TwinLaw = np.ones(1)
1034#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
1035#    TwDict = refDict.get('TwDict',{})           
1036#    if 'S' in calcControls[hfx+'histType']:
1037#        NTL = calcControls[phfx+'NTL']
1038#        NM = calcControls[phfx+'TwinNMN']+1
1039#        TwinLaw = calcControls[phfx+'TwinLaw']
1040#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1041#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1042    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1043        GetAtomFXU(pfx,calcControls,parmDict)
1044    if not Xdata.size:          #no atoms in phase!
1045        return
1046    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
1047    Gdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
1048    if SGData['SGInv'] and not SGData['SGFixed']:
1049        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1050    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering--> [Mxyz,nops,natms]
1051    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
1052    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1053    Kdata = np.inner(Gdata.T,uAmat).T     #Cartesian unit vectors
1054    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1055    Kdata /= Kmean
1056    Uij = np.array(G2lat.U6toUij(Uijdata))
1057    bij = Mast*Uij.T
1058    blkSize = 100       #no. of reflections in a block - size seems optimal
1059    nRef = refDict['RefList'].shape[0]
1060    SQ = 1./(2.*refDict['RefList'].T[4])**2
1061    refDict['FF']['El'] = list(MFtables.keys())
1062    refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
1063    for iel,El in enumerate(refDict['FF']['El']):
1064        refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1065#reflection processing begins here - big arrays!
1066    iBeg = 0
1067    while iBeg < nRef:
1068        iFin = min(iBeg+blkSize,nRef)
1069        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1070        H = refl.T[:3].T                          #array(blkSize,3)
1071#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
1072#        TwMask = np.any(H,axis=-1)
1073#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1074#            for ir in range(blkSize):
1075#                iref = ir+iBeg
1076#                if iref in TwDict:
1077#                    for i in TwDict[iref]:
1078#                        for n in range(NTL):
1079#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1080#            TwMask = np.any(H,axis=-1)
1081        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
1082        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1083        Uniq = np.inner(H,SGMT)
1084        Phi = np.inner(H,SGT)
1085        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1086        biso = -SQfactor*Uisodata[:,nxs]
1087        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
1088        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1089        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1090        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1091        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1092        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1093        if SGData['SGInv']:
1094            if not SGData['SGFixed']:
1095                mphase = np.hstack((phase,-phase))  #OK
1096            else:
1097                mphase = phase
1098        else:
1099            mphase = phase                    #
1100        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1101        mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
1102        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1103        cosm = np.cos(mphase)                               #ditto
1104        HM = np.inner(Bmat.T,H)                             #put into cartesian space
1105        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Kdata = MAGS & HM = UVEC in magstrfc.for both OK
1106        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
1107        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1108        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1109        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
1110        fams = np.sum(np.sum(fam,axis=-1),axis=-1)     #Mxyz,Nref  Sum(sum(fam,atoms),ops)
1111        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)     #ditto
1112        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)   #Sum(fams**2,Mxyz) Re + Im
1113        refl.T[7] = np.copy(refl.T[9])               
1114        refl.T[10] = atan2d(fbms[0],fams[0]) #- what is phase for mag refl?
1115#        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
1116#            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
1117#            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1118#        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
1119#            if len(TwinLaw) > 1:
1120#                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
1121#                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
1122#                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
1123#                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
1124#            else:   # checked correct!!
1125#                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
1126#                refl.T[7] = np.copy(refl.T[9])               
1127#                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1128##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
1129        iBeg += blkSize
1130#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
1131    return copy.deepcopy(refDict['RefList'])
1132
1133def MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1134    '''Compute magnetic structure factor derivatives numerically - for powders/nontwins only
1135    input:
1136   
1137    :param dict refDict: where
1138        'RefList' list where each ref = h,k,l,it,d,...
1139        'FF' dict of form factors - filled in below
1140    :param np.array G:      reciprocal metric tensor
1141    :param str hfx:    histogram id string
1142    :param str pfx:    phase id string
1143    :param dict SGData: space group info. dictionary output from SpcGroup
1144    :param dict calcControls:
1145    :param dict parmDict:
1146   
1147    :returns: dict dFdvDict: dictionary of magnetic derivatives
1148    '''
1149   
1150    trefDict = copy.deepcopy(refDict)
1151    dM = 1.e-6
1152    dFdvDict = {}
1153    for parm in parmDict:
1154        if 'AM' in parm:
1155            parmDict[parm] += dM
1156            prefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
1157            parmDict[parm] -= 2*dM
1158            mrefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
1159            parmDict[parm] += dM
1160            dFdvDict[parm] = (prefList[:,9]-mrefList[:,9])/(2.*dM)
1161    return dFdvDict
1162           
1163def MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1164    '''Compute nonmagnetic structure factor derivatives on blocks of reflections in magnetic structures - for powders/nontwins only
1165    input:
1166   
1167    :param dict refDict: where
1168        'RefList' list where each ref = h,k,l,it,d,...
1169        'FF' dict of form factors - filled in below
1170    :param np.array G:      reciprocal metric tensor
1171    :param str hfx:    histogram id string
1172    :param str pfx:    phase id string
1173    :param dict SGData: space group info. dictionary output from SpcGroup
1174    :param dict calcControls:
1175    :param dict parmDict:
1176   
1177    :returns: dict dFdvDict: dictionary of derivatives
1178    '''
1179   
1180    g = nl.inv(G)
1181    ast = np.sqrt(np.diag(G))
1182    ainv = np.sqrt(np.diag(g))
1183    GS = G/np.outer(ast,ast)
1184    Ginv = g/np.outer(ainv,ainv)
1185    uAmat = G2lat.Gmat2AB(GS)[0]
1186    Mast = twopisq*np.multiply.outer(ast,ast)
1187    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1188    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1189    Ncen = len(SGData['SGCen'])
1190    Nops = len(SGMT)*Ncen
1191    if not SGData['SGFixed']:
1192        Nops *= (1+SGData['SGInv'])
1193    Bmat = G2lat.Gmat2AB(G)[1]
1194    nRef = len(refDict['RefList'])
1195    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1196        GetAtomFXU(pfx,calcControls,parmDict)
1197    if not Xdata.size:          #no atoms in phase!
1198        return {}
1199    mSize = len(Mdata)
1200    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
1201    Gones = np.ones_like(Gdata)
1202    Gdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
1203    Gones = np.inner(Gones.T,SGMT).T
1204    if SGData['SGInv'] and not SGData['SGFixed']:
1205        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
1206        Gones = np.hstack((Gones,-Gones))       #inversion if any
1207    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
1208    Gones = np.hstack([Gones for icen in range(Ncen)])        #dup over cell centering
1209    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
1210    Gones = SGData['MagMom'][nxs,:,nxs]*Gones   #flip vectors according to spin flip
1211    Mag = np.tile(Mag[:,nxs],Nops).#make Mag same length as Gdata
1212    Kdata = np.inner(Gdata.T,uAmat).T     #Cartesian unit vectors
1213    Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
1214    Kdata /= Kmean
1215    Uij = np.array(G2lat.U6toUij(Uijdata))
1216    bij = Mast*Uij.T
1217    dFdvDict = {}
1218    dFdfr = np.zeros((nRef,mSize))
1219    dFdx = np.zeros((nRef,mSize,3))
1220    dFdui = np.zeros((nRef,mSize))
1221    dFdua = np.zeros((nRef,mSize,6))
1222    time0 = time.time()
1223#reflection processing begins here - big arrays!
1224    iBeg = 0
1225    blkSize = 5       #no. of reflections in a block - optimized for speed
1226    while iBeg < nRef:
1227        iFin = min(iBeg+blkSize,nRef)
1228        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1229        H = refl.T[:3].T
1230        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1231        SQfactor = 8.0*SQ*np.pi**2
1232        Uniq = np.inner(H,SGMT)             # array(nSGOp,3)
1233        Phi = np.inner(H,SGT)
1234        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1235        occ = Mdata*Fdata/Nops
1236        biso = -SQfactor*Uisodata[:,nxs]
1237        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
1238        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1239        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
1240        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1241        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))
1242        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1243        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1244        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1245        if SGData['SGInv']:
1246            if not SGData['SGFixed']:
1247                mphase = np.hstack((phase,-phase))  #OK
1248                Uniq = np.hstack((Uniq,-Uniq))      #Nref,Nops,hkl
1249                Hij = np.hstack((Hij,Hij))
1250            else:
1251                mphase = phase
1252        else:
1253            mphase = phase                    #
1254        Hij = np.concatenate(np.array([Hij for cen in SGData['SGCen']]),axis=1)
1255        Uniq = np.hstack([Uniq for cen in SGData['SGCen']])
1256        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
1257        mphase = np.concatenate(mphase,axis=1)              #Nref,Nop,Natm
1258        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
1259        cosm = np.cos(mphase)                               #ditto
1260        HM = np.inner(Bmat.T,H)                             #put into cartesian space
1261        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit cartesian vector for H
1262        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
1263        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
1264       
1265        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
1266        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]
1267        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                      #Mxyz,Nref
1268        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)
1269        famx = -Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*sinm[nxs,:,:,:]   #Mxyz,Nref,Nops,Natom
1270        fbmx = Q*TMcorr[nxs,:,nxs,:]*Mag[nxs,nxs,:,:]*cosm[nxs,:,:,:]
1271        #sums below are over Nops - real part
1272        dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK
1273        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK
1274        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
1275        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK
1276        # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
1277        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
1278        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2)
1279        dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
1280        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
1281        #accumulate derivatives   
1282        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
1283        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
1284        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
1285        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok
1286        iBeg += blkSize
1287    print (' %d derivative time %.4f\r'%(nRef,time.time()-time0))
1288        #loop over atoms - each dict entry is list of derivatives for all the reflections
1289    for i in range(len(Mdata)):
1290        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1291        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1292        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1293        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1294        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1295        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1296        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1297        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1298        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1299        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1300        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1301    return dFdvDict
1302       
1303def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
1304    '''Compute structure factor derivatives on blocks of reflections - for twins only
1305    faster than StructureFactorDervTw
1306    input:
1307   
1308    :param dict refDict: where
1309        'RefList' list where each ref = h,k,l,it,d,...
1310        'FF' dict of form factors - filled in below
1311    :param np.array G:      reciprocal metric tensor
1312    :param str hfx:    histogram id string
1313    :param str pfx:    phase id string
1314    :param dict SGData: space group info. dictionary output from SpcGroup
1315    :param dict calcControls:
1316    :param dict parmDict:
1317   
1318    :returns: dict dFdvDict: dictionary of derivatives
1319    '''
1320    phfx = pfx.split(':')[0]+hfx
1321    ast = np.sqrt(np.diag(G))
1322    Mast = twopisq*np.multiply.outer(ast,ast)
1323    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1324    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1325    FFtables = calcControls['FFtables']
1326    BLtables = calcControls['BLtables']
1327    TwDict = refDict.get('TwDict',{})           
1328    NTL = calcControls[phfx+'NTL']
1329    NM = calcControls[phfx+'TwinNMN']+1
1330    TwinLaw = calcControls[phfx+'TwinLaw']
1331    TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1332    TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1333    nTwin = len(TwinLaw)       
1334    nRef = len(refDict['RefList'])
1335    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1336        GetAtomFXU(pfx,calcControls,parmDict)
1337    if not Xdata.size:          #no atoms in phase!
1338        return {}
1339    mSize = len(Mdata)
1340    FF = np.zeros(len(Tdata))
1341    if 'NC' in calcControls[hfx+'histType']:
1342        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1343    elif 'X' in calcControls[hfx+'histType']:
1344        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1345        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1346    Uij = np.array(G2lat.U6toUij(Uijdata))
1347    bij = Mast*Uij.T
1348    dFdvDict = {}
1349    dFdfr = np.zeros((nRef,nTwin,mSize))
1350    dFdx = np.zeros((nRef,nTwin,mSize,3))
1351    dFdui = np.zeros((nRef,nTwin,mSize))
1352    dFdua = np.zeros((nRef,nTwin,mSize,6))
1353    dFdbab = np.zeros((nRef,nTwin,2))
1354    dFdtw = np.zeros((nRef,nTwin))
1355    time0 = time.time()
1356#reflection processing begins here - big arrays!
1357    iBeg = 0
1358    blkSize = 16       #no. of reflections in a block - optimized for speed
1359    while iBeg < nRef:
1360        iFin = min(iBeg+blkSize,nRef)
1361        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1362        H = refl.T[:3]
1363        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
1364        TwMask = np.any(H,axis=-1)
1365        for ir in range(blkSize):
1366            iref = ir+iBeg
1367            if iref in TwDict:
1368                for i in TwDict[iref]:
1369                    for n in range(NTL):
1370                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1371        TwMask = np.any(H,axis=-1)
1372        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
1373        SQfactor = 8.0*SQ*np.pi**2
1374        if 'T' in calcControls[hfx+'histType']:
1375            if 'P' in calcControls[hfx+'histType']:
1376                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1377            else:
1378                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1379            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
1380            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
1381        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1382        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
1383        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1384        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
1385        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
1386        Phi = np.inner(H,SGT)
1387        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
1388        sinp = np.sin(phase)
1389        cosp = np.cos(phase)
1390        occ = Mdata*Fdata/len(SGT)
1391        biso = -SQfactor*Uisodata[:,nxs]
1392        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)
1393        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
1394        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
1395        Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,nTwin,len(SGT),6))
1396        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1397        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
1398        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
1399        fotp = FPP*Tcorr       
1400        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1401            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
1402            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
1403        else:
1404            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
1405            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
1406        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
1407        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
1408        if SGData['SGInv']: #centrosymmetric; B=0
1409            fbs[0] *= 0.
1410            fas[1] *= 0.
1411        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
1412        fbx = np.array([fot*cosp,-fotp*sinp])
1413        #sum below is over Uniq
1414        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
1415        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1416        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
1417        dfadx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fax[:,:,:,:,:,nxs],axis=-3),axis=0) # nRef x nTwin x nAtoms x xyz; sum on ops & A,A'
1418        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 
1419        if not SGData['SGInv']:
1420            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
1421            dfadba /= 2.
1422#            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
1423            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
1424            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
1425            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
1426        else:
1427            dfbdfr = np.zeros_like(dfadfr)
1428            dfbdx = np.zeros_like(dfadx)
1429            dfbdui = np.zeros_like(dfadui)
1430            dfbdua = np.zeros_like(dfadua)
1431#            dfbdba = np.zeros_like(dfadba)
1432        SA = fas[0]+fas[1]
1433        SB = fbs[0]+fbs[1]
1434        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
1435        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
1436        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
1437        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
1438        if SGData['SGInv']: #centrosymmetric; B=0
1439            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
1440        else:               
1441            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
1442#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1443#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1444        iBeg += blkSize
1445    print (' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0))
1446    #loop over atoms - each dict entry is list of derivatives for all the reflections
1447    for i in range(len(Mdata)):     #these all OK
1448        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
1449        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
1450        dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)
1451        dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)
1452        dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)
1453        dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)
1454        dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)
1455        dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)
1456        dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)
1457        dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)
1458        dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)
1459    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
1460    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
1461    for i in range(nTwin):
1462        dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
1463    return dFdvDict
1464   
1465def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1466    '''
1467    Compute super structure factors for all h,k,l,m for phase - no twins
1468    puts the result, F^2, in each ref[9] in refList
1469    works on blocks of 32 reflections for speed
1470    input:
1471   
1472    :param dict refDict: where
1473        'RefList' list where each ref = h,k,l,m,it,d,...
1474        'FF' dict of form factors - filed in below
1475    :param np.array G:      reciprocal metric tensor
1476    :param str pfx:    phase id string
1477    :param dict SGData: space group info. dictionary output from SpcGroup
1478    :param dict calcControls:
1479    :param dict ParmDict:
1480
1481    '''
1482    phfx = pfx.split(':')[0]+hfx
1483    ast = np.sqrt(np.diag(G))
1484    GS = G/np.outer(ast,ast)
1485    uAmat = G2lat.Gmat2AB(GS)[0]
1486    Mast = twopisq*np.multiply.outer(ast,ast)   
1487    SGInv = SGData['SGInv']
1488    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1489    Ncen = len(SGData['SGCen'])
1490    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
1491    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1492    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1493    SSCen = SSGData['SSGCen']
1494    FFtables = calcControls['FFtables']
1495    BLtables = calcControls['BLtables']
1496    MFtables = calcControls['MFtables']
1497    Amat,Bmat = G2lat.Gmat2AB(G)
1498    Flack = 1.0
1499    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1500        Flack = 1.-2.*parmDict[phfx+'Flack']
1501    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1502        GetAtomFXU(pfx,calcControls,parmDict)
1503    if not Xdata.size:          #no atoms in phase!
1504        return
1505    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1506    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
1507    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1508
1509    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms
1510       
1511        GSdata0 = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1512        if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1513            GSdata0 = np.hstack((GSdata0,-GSdata0))     
1514        GSdata0 = np.hstack([GSdata0 for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1515        GSdata0 = SGData['MagMom'][nxs,:,nxs]*GSdata0   #flip vectors according to spin flip * det(opM)
1516        GSdata0 = np.swapaxes(GSdata0,0,1)    #Nop,Natm,Mxyz
1517        GSdata0 = np.inner(GSdata0,uAmat.T)   #--> cartesian
1518       
1519        mXYZ = np.array([[xyz[0] for xyz in list(G2spc.GenAtom(xyz,SGData,All=True,Move=True))] for xyz in (Xdata+dXdata).T])%1. #Natn,Nop,xyz
1520        Mmod,SMag = G2mth.MagMod(ngl,mXYZ,modQ,MSSdata,SGData,SSGData)   #Re cos/Im sin,Nops,Natm,Nwaves,Mxyz
1521#        MmodA = np.inner(MmodA,uAmat.T)   #make cartesian
1522#        MmodB = np.inner(MmodB,uAmat.T)
1523
1524#from #3812
1525#1st try at this
1526#        GSdata = Gdata.T[:,nxs,:]+Mmod    #Natm,ntau,Mxyz
1527#        GSdata = np.inner(GSdata,SGMT).T  #apply sym. ops.--> Mxyz,Nops,Ntau,Natm
1528#        if SGData['SGInv'] and not SGData['SGFixed']:
1529#            GSdata = np.hstack((GSdata,-GSdata))       #inversion if any
1530#        GSdata = np.hstack([GSdata for icen in range(Ncen)])        #dup over cell centering
1531#        GSdata = SGData['MagMom'][nxs,:,nxs,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1532#        Kdata = np.inner(GSdata.T,uAmat).T     #Cartesian unit vectors
1533#        SMag = np.sqrt(np.sum(Kdata**2,axis=0))         #Nops,Ntau,Natm
1534#        Kmean = np.mean(SMag,axis=0)    #normalization --> unit vectors
1535#        Kdata /= Kmean[nxs,nxs,:,:]      #mxyz,nops,ntau,natm
1536       
1537    FF = np.zeros(len(Tdata))
1538    if 'NC' in calcControls[hfx+'histType']:
1539        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1540    elif 'X' in calcControls[hfx+'histType']:
1541        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1542        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1543    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1544    bij = Mast*Uij
1545    blkSize = 48       #no. of reflections in a block
1546    nRef = refDict['RefList'].shape[0]
1547    SQ = 1./(2.*refDict['RefList'].T[5])**2
1548    if 'N' in calcControls[hfx+'histType']:
1549        dat = G2el.getBLvalues(BLtables)
1550        refDict['FF']['El'] = list(dat.keys())
1551        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1552        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1553        for iel,El in enumerate(refDict['FF']['El']):
1554            if El in MFtables:
1555                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1556    else:
1557        dat = G2el.getFFvalues(FFtables,0.)
1558        refDict['FF']['El'] = list(dat.keys())
1559        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1560        for iel,El in enumerate(refDict['FF']['El']):
1561            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1562#    time0 = time.time()
1563#reflection processing begins here - big arrays!
1564    iBeg = 0
1565    while iBeg < nRef:
1566        iFin = min(iBeg+blkSize,nRef)
1567        mRef = iFin-iBeg
1568        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1569        H = refl.T[:4]                          #array(blkSize,4)
1570        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1571        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1572        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1573        Uniq = np.inner(H.T,SSGMT)
1574        UniqP = np.inner(HP.T,SGMT)
1575        Phi = np.inner(H.T,SSGT)
1576        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1577            Uniq = np.hstack((Uniq,-Uniq))
1578            Phi = np.hstack((Phi,-Phi))
1579            UniqP = np.hstack((UniqP,-UniqP))
1580        if 'T' in calcControls[hfx+'histType']:
1581            if 'P' in calcControls[hfx+'histType']:
1582                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1583            else:
1584                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1585            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1586            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1587        Bab = 0.
1588        if phfx+'BabA' in parmDict:
1589            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1590        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1591        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1592        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1593        phase = np.hstack([phase for cen in SSCen])
1594        sinp = np.sin(phase)
1595        cosp = np.cos(phase)
1596        biso = -SQfactor*Uisodata[:,nxs]
1597        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1598        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1599        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1600        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1601
1602        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: mag math here??
1603            phasem = twopi*np.inner(H.T[:,:3],mXYZ)
1604            phasem = np.swapaxes(phasem,1,2)
1605            cosm = np.cos(phasem)
1606            sinm = np.sin(phasem)
1607            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1608            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
1609                     
1610            HM = np.inner(Bmat,HP.T)                            #put into cartesian space
1611            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #& normalize
1612#for fixed moments --> m=0 reflections                       
1613            fam0 = TMcorr[:,nxs,:,nxs]*GSdata0[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1614            fbm0 = TMcorr[:,nxs,:,nxs]*GSdata0[nxs,:,:,:]*sinm[:,:,:,nxs]   
1615                       
1616            famq0 = np.sum(np.sum(fam0,axis=-2),axis=-2)        #Nref,Mxyz; sum ops & atoms
1617            fbmq0 = np.sum(np.sum(fbm0,axis=-2),axis=-2)
1618           
1619            fas0 = np.sum(famq0,axis=-1)**2-np.sum(HM.T*famq0,axis=-1)**2   #mag intensity calc F^2-(e.F)^2
1620            fbs0 = np.sum(fbmq0,axis=-1)**2-np.sum(HM.T*fbmq0,axis=-1)**2
1621#for modulated moments --> m != 0 reflections
1622#            M = np.array(np.abs(H[3]),dtype=np.int)-1
1623#            fam = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(MmodA[:,:,M[i],:]*cosm[i,:,:,nxs]-np.sign(H[3])[i,nxs,nxs,nxs]*MmodB[:,:,M[i],:]*sinm[i,:,:,nxs]),0.) for i in range(mRef)])
1624#            fbm = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(MmodA[:,:,M[i],:]*sinm[i,:,:,nxs]+np.sign(H[3])[i,nxs,nxs,nxs]*MmodB[:,:,M[i],:]*cosm[i,:,:,nxs]),0.) for i in range(mRef)])
1625#                       
1626#            famq = np.sum(np.sum(fam/2.,axis=-2),axis=-2)      #Nref,Mxyz; sum ops & atoms
1627#            fbmq = np.sum(np.sum(fbm/2.,axis=-2),axis=-2)
1628#           
1629#            fas = np.sum(famq,axis=-1)**2-np.sum(HM.T*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
1630#            fbs = np.sum(fbmq,axis=-1)**2-np.sum(HM.T*fbmq,axis=-1)**2
1631#from #3812
1632            D = twopi*H.T[:,3:]*glTau[nxs,:]
1633            mphase = phase[:,:,nxs,:]+D[:,nxs,:,nxs]
1634            mphase = np.array([mphase+twopi*np.inner(cen,HP.T)[:,nxs,nxs,nxs] for cen in SGData['SGCen']])
1635            mphase = np.concatenate(mphase,axis=1)    #remove extra axis; Nref,Nop,Natm
1636            sinm = np.swapaxes(np.sin(mphase),2,3)    #--> Nref,Nop,Natm,Ntau
1637            cosm = np.swapaxes(np.cos(mphase),2,3)                               #ditto
1638           
1639            HM = np.inner(Bmat,HP.T)                             #put into cartesian space
1640            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
1641            eDotK = np.sum(HM[:,:,nxs,nxs,nxs]*Kdata[:,nxs,:,:,:],axis=0)
1642            Q = HM[:,:,nxs,nxs,nxs,nxs]*eDotK[nxs,:,:,:,:]-Kdata[:,nxs,:,:,:] #Mxyz,Nref,Nop,Ntau,Natm
1643
1644            fam = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*cosm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,nxs])   #Mxyz,Nref,Nop,Natm,Ntau,ReIm
1645            fbm = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*sinm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,nxs])
1646           
1647            fams = np.sum(np.sum(np.sum(fam,axis=-1),axis=2),axis=2)      #xyz,Nref,ntau; sum ops & atoms
1648            fbms = np.sum(np.sum(np.sum(fbm,axis=-1),axis=2),axis=2)      #ditto
1649           
1650            fas = np.sum(fams*glWt[nxs,nxs,:],axis=-1)
1651            fbs = np.sum(fbms*glWt[nxs,nxs,:],axis=-1)
1652                       
1653            refl.T[10] = np.where(H[3],fas+fbs,fas0+fbs0)
1654            refl.T[11] = np.where(H[3],atan2d(fas,fbs),atan2d(fas0,fbs0))
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])  #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 parm == '0':
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)
2599    elif xfac > 1.:
2600        xfac2 = 1./np.sqrt(xfac)
2601        dlde = 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,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']*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,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                im = 0
3041                if Phase['General'].get('Modulated',False):
3042                    im = 1
3043                pId = Phase['pId']
3044                phfx = '%d:%d:'%(pId,hId)
3045                refDict = refLists[phase]
3046                sumFo = 0.0
3047                sumdF = 0.0
3048                sumFosq = 0.0
3049                sumdFsq = 0.0
3050                sumInt = 0.0
3051                nExcl = 0
3052                # test to see if we are using multiprocessing below
3053                useMP,ncores = G2mp.InitMP()
3054                if len(refDict['RefList']) < 100: useMP = False       
3055                if useMP: # multiprocessing: create a set of initialized Python processes
3056                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3057                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3058                    profArgs = [[] for i in range(G2mp.ncores)]
3059                else:
3060                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3061                if 'C' in calcControls[hfx+'histType']:
3062                    # are we multiprocessing?
3063                    for iref,refl in enumerate(refDict['RefList']):
3064                        if useMP: 
3065                            profArgs[iref%G2mp.ncores].append((refl,iref))
3066                        else:
3067                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3068                            if type(icod) is tuple:
3069                                refl[8+im] = icod[0]
3070                                sumInt += icod[1]
3071                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3072                            elif icod == -1:
3073                                refl[3+im] *= -1
3074                                nExcl += 1
3075                            elif icod == -2:
3076                                break
3077                    if useMP:
3078                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3079                            sumInt += sInt
3080                            for refl8im,irefl in resList:
3081                                if refl8im is None:
3082                                    refDict['RefList'][irefl][3+im] *= -1
3083                                    nExcl += 1
3084                                else:
3085                                    refDict['RefList'][irefl][8+im] = refl8im
3086                                    if parmDict[phfx+'LeBail']:
3087                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3088                elif 'T' in calcControls[hfx+'histType']:
3089                    for iref,refl in enumerate(refDict['RefList']):
3090                        if useMP: 
3091                            profArgs[iref%G2mp.ncores].append((refl,iref))
3092                        else:
3093                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3094                            if type(icod) is tuple:
3095                                refl[8+im] = icod[0]
3096                                sumInt += icod[1]
3097                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3098                            elif icod == -1:
3099                                refl[3+im] *= -1
3100                                nExcl += 1
3101                            elif icod == -2:
3102                                break
3103                    if useMP:
3104                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3105                            sumInt += sInt
3106                            for refl8im,irefl in resList:
3107                                if refl8im is None:
3108                                    refDict['RefList'][irefl][3+im] *= -1
3109                                    nExcl += 1
3110                                else:
3111                                    refDict['RefList'][irefl][8+im] = refl8im
3112                                    if parmDict[phfx+'LeBail']:
3113                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3114                if useMP: MPpool.terminate()
3115                sumFo = 0.0
3116                sumdF = 0.0
3117                sumFosq = 0.0
3118                sumdFsq = 0.0
3119                for iref,refl in enumerate(refDict['RefList']):
3120                    Fo = np.sqrt(np.abs(refl[8+im]))
3121                    Fc = np.sqrt(np.abs(refl[9]+im))
3122                    sumFo += Fo
3123                    sumFosq += refl[8+im]**2
3124                    sumdF += np.abs(Fo-Fc)
3125                    sumdFsq += (refl[8+im]-refl[9+im])**2
3126                if sumFo:
3127                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3128                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3129                else:
3130                    Histogram['Residuals'][phfx+'Rf'] = 100.
3131                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3132                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3133                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3134                Histogram['Residuals']['hId'] = hId
3135        elif 'HKLF' in histogram[:4]:
3136            Histogram = Histograms[histogram]
3137            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3138    if GSASIIpath.GetConfigValue('Show_timing',False):
3139        print ('GetFobsSq t=',time.time()-starttime)
3140               
3141def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3142    'Computes the powder pattern for a histogram based on contributions from all used phases'
3143    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3144   
3145    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3146        U = parmDict[hfx+'U']
3147        V = parmDict[hfx+'V']
3148        W = parmDict[hfx+'W']
3149        X = parmDict[hfx+'X']
3150        Y = parmDict[hfx+'Y']
3151        Z = parmDict[hfx+'Z']
3152        tanPos = tand(refl[5+im]/2.0)
3153        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3154        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3155        sig = max(0.001,sig)
3156        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3157        gam = max(0.001,gam)
3158        return sig,gam
3159               
3160    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3161        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3162            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3163        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3164        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3165        sig += Ssig
3166        gam += Sgam
3167        return sig,gam
3168       
3169    def GetReflAlpBet(refl,im,hfx,parmDict):
3170        alp = parmDict[hfx+'alpha']/refl[4+im]
3171        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3172        return alp,bet
3173       
3174    hId = Histogram['hId']
3175    hfx = ':%d:'%(hId)
3176    bakType = calcControls[hfx+'bakType']
3177    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
3178    yc = np.zeros_like(yb)
3179    cw = np.diff(ma.getdata(x))
3180    cw = np.append(cw,cw[-1])
3181       
3182    if 'C' in calcControls[hfx+'histType']:   
3183        shl = max(parmDict[hfx+'SH/L'],0.002)
3184        Ka2 = False
3185        if hfx+'Lam1' in (parmDict.keys()):
3186            wave = parmDict[hfx+'Lam1']
3187            Ka2 = True
3188            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3189            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3190        else:
3191            wave = parmDict[hfx+'Lam']
3192    else:
3193        shl = 0.
3194    for phase in Histogram['Reflection Lists']:
3195        refDict = Histogram['Reflection Lists'][phase]
3196        if phase not in Phases:     #skips deleted or renamed phases silently!
3197            continue
3198        Phase = Phases[phase]
3199        pId = Phase['pId']
3200        pfx = '%d::'%(pId)
3201        phfx = '%d:%d:'%(pId,hId)
3202        hfx = ':%d:'%(hId)
3203        SGData = Phase['General']['SGData']
3204        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3205        im = 0
3206        if Phase['General'].get('Modulated',False):
3207            SSGData = Phase['General']['SSGData']
3208            im = 1  #offset in SS reflection list
3209            #??
3210        Dij = GetDij(phfx,SGData,parmDict)
3211        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do someting if Dij << 0.
3212        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3213        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3214            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3215        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3216        Vst = np.sqrt(nl.det(G))    #V*
3217        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3218            if im:
3219                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3220            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3221                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3222            else:
3223                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3224        badPeak = False
3225        # test to see if we are using multiprocessing here
3226        useMP,ncores = G2mp.InitMP()
3227        if len(refDict['RefList']) < 100: useMP = False       
3228        if useMP: # multiprocessing: create a set of initialized Python processes
3229            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3230            profArgs = [[] for i in range(ncores)]
3231        if 'C' in calcControls[hfx+'histType']:
3232            for iref,refl in enumerate(refDict['RefList']):
3233                if im:
3234                    h,k,l,m = refl[:4]
3235                else:
3236                    h,k,l = refl[:3]
3237                Uniq = np.inner(refl[:3],SGMT)
3238                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
3239                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3240                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3241                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3242                refl[11+im] *= Vst*Lorenz
3243                 
3244                if Phase['General'].get('doPawley'):
3245                    try:
3246                        if im:
3247                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3248                        else:
3249                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3250                        refl[9+im] = parmDict[pInd]
3251                    except KeyError:
3252#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3253                        continue
3254                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3255                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3256                iFin = np.searchsorted(x,refl[5+im]+fmax)
3257                if not iBeg+iFin:       #peak below low limit - skip peak
3258                    continue
3259                elif not iBeg-iFin:     #peak above high limit - done
3260                    break
3261                elif iBeg > iFin:   #bad peak coeff - skip
3262                    badPeak = True
3263                    continue
3264                if useMP:
3265                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3266                else:
3267                    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
3268                if Ka2:
3269                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3270                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3271                    iBeg = np.searchsorted(x,pos2-fmin)
3272                    iFin = np.searchsorted(x,pos2+fmax)
3273                    if not iBeg+iFin:       #peak below low limit - skip peak
3274                        continue
3275                    elif not iBeg-iFin:     #peak above high limit - done
3276                        return yc,yb
3277                    elif iBeg > iFin:   #bad peak coeff - skip
3278                        continue
3279                    if useMP:
3280                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3281                    else:
3282                        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
3283        elif 'T' in calcControls[hfx+'histType']:
3284            for iref,refl in enumerate(refDict['RefList']):
3285                if im:
3286                    h,k,l,m = refl[:4]
3287                else:
3288                    h,k,l = refl[:3]
3289                Uniq = np.inner(refl[:3],SGMT)
3290                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
3291                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
3292#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
3293                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3294                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3295                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3296                refl[11+im] *= Vst*Lorenz
3297                if Phase['General'].get('doPawley'):
3298                    try:
3299                        if im:
3300                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3301                        else:
3302                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3303                        refl[9+im] = parmDict[pInd]
3304                    except KeyError:
3305#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3306                        continue
3307                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
3308                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3309                iFin = np.searchsorted(x,refl[5+im]+fmax)
3310                if not iBeg+iFin:       #peak below low limit - skip peak
3311                    continue
3312                elif not iBeg-iFin:     #peak above high limit - done
3313                    break
3314                elif iBeg > iFin:   #bad peak coeff - skip
3315                    badPeak = True
3316                    continue
3317                if useMP:
3318                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin))
3319                else:
3320                    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]
3321#        print 'profile calc time: %.3fs'%(time.time()-time0)
3322        if useMP and 'C' in calcControls[hfx+'histType']:
3323            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
3324                yc += y
3325            MPpool.terminate()
3326        elif useMP:
3327            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
3328                yc += y
3329            MPpool.terminate()
3330    if badPeak:
3331        print ('ouch #4 bad profile coefficients yield negative peak width; some reflections skipped')
3332    if GSASIIpath.GetConfigValue('Show_timing',False):
3333        print ('getPowderProfile t=%.3f'%(time.time()-starttime))
3334    return yc,yb
3335   
3336def getPowderProfileDervMP(args):
3337    '''Computes the derivatives of the computed powder pattern with respect to all
3338    refined parameters.
3339    Multiprocessing version.
3340    '''
3341    import pytexture as ptx
3342    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics for each processor
3343    parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars = args[:9]
3344    prc=0
3345    tprc=1
3346    if len(args)