source: trunk/GSASIIstrMath.py @ 3805

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

fix normalization of mag moments --> Kvector
add save of map peaks as csv file

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