source: trunk/GSASIIstrMath.py @ 3802

Last change on this file since 3802 was 3802, checked in by toby, 4 years ago

generalized restraints completed

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