source: trunk/GSASIIstrMath.py @ 3859

Last change on this file since 3859 was 3859, checked in by toby, 3 years ago

Use Dij values for SeqRef? table even if not varied

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