source: trunk/GSASIIstrMath.py @ 4068

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

fix bug in CalcBack? when no background file specified

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