source: trunk/GSASIIstrMath.py @ 3795

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

fix space group & incommensurate space group displays for mag. incomensurate cases
disable cell refine in Unit Cell List if no peaks
another shot at incommensurate mag. structure factor calc.

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