source: trunk/GSASIIstrMath.py @ 4129

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

minor IMSF change

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