source: trunk/GSASIIstrMath.py @ 4054

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

changes to incommensurate magnetism math

  • 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-07-07 00:00:41 +0000 (Sun, 07 Jul 2019) $
8# $Author: vondreele $
9# $Revision: 4054 $
10# $URL: trunk/GSASIIstrMath.py $
11# $Id: GSASIIstrMath.py 4054 2019-07-07 00:00:41Z 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: 4054 $")
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        Tmag,TmagA,TmagB = G2mth.MagMod2(mXYZ,modQ,MSSdata,SGData,SSGData)   #Nops,Natm,Mxyz-Tmag matches drawing moments
1513       
1514        if not SGData['SGGray']:    #for fixed Mx,My,Mz
1515            Mmod += Gdata.T[:,nxs,:]
1516            GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
1517            if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
1518                GSdata = np.hstack((GSdata,-GSdata))     
1519            GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
1520            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
1521            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
1522
1523    FF = np.zeros(len(Tdata))
1524    if 'NC' in calcControls[hfx+'histType']:
1525        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1526    elif 'X' in calcControls[hfx+'histType']:
1527        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1528        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1529    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1530    bij = Mast*Uij
1531    blkSize = 48       #no. of reflections in a block
1532    nRef = refDict['RefList'].shape[0]
1533    SQ = 1./(2.*refDict['RefList'].T[5])**2
1534    if 'N' in calcControls[hfx+'histType']:
1535        dat = G2el.getBLvalues(BLtables)
1536        refDict['FF']['El'] = list(dat.keys())
1537        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1538        refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1539        for iel,El in enumerate(refDict['FF']['El']):
1540            if El in MFtables:
1541                refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1542    else:
1543        dat = G2el.getFFvalues(FFtables,0.)
1544        refDict['FF']['El'] = list(dat.keys())
1545        refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1546        for iel,El in enumerate(refDict['FF']['El']):
1547            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1548#    time0 = time.time()
1549#reflection processing begins here - big arrays!
1550    iBeg = 0
1551    while iBeg < nRef:
1552        iFin = min(iBeg+blkSize,nRef)
1553        mRef= iFin-iBeg
1554        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1555        H = refl.T[:4]                          #array(blkSize,4)
1556        HP = H[:3]+modQ[:,nxs]*H[3:]            #projected hklm to hkl
1557        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1558        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1559        Uniq = np.inner(H.T,SSGMT)
1560        UniqP = np.inner(HP.T,SGMT)
1561        Phi = np.inner(H.T,SSGT)
1562        if SGInv and not SGData['SGFixed']:   #if centro - expand HKL sets
1563            Uniq = np.hstack((Uniq,-Uniq))
1564            Phi = np.hstack((Phi,-Phi))
1565            UniqP = np.hstack((UniqP,-UniqP))
1566        if 'T' in calcControls[hfx+'histType']:
1567            if 'P' in calcControls[hfx+'histType']:
1568                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1569            else:
1570                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1571            FP = np.repeat(FP.T,Uniq.shape[1],axis=0)
1572            FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0)
1573        Bab = 0.
1574        if phfx+'BabA' in parmDict:
1575            Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1])
1576        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1577        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0)
1578        phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])
1579        phase = np.hstack([phase for cen in SSCen])
1580        sinp = np.sin(phase)
1581        cosp = np.cos(phase)
1582        biso = -SQfactor*Uisodata[:,nxs]
1583        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T
1584        HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1585        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1586        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1587
1588        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: mag math here??           
1589           
1590            phasem = twopi*np.inner(H.T[:,:3],mXYZ)
1591            phasem = np.swapaxes(phasem,1,2)
1592            cosm = np.cos(phasem)
1593            sinm = np.sin(phasem)
1594            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
1595            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2.*Nops)     #Nref,Natm
1596#                     
1597            HM = np.inner(Bmat,HP.T)                            #put into cartesian space
1598            eM = HM/np.sqrt(np.sum(HM**2,axis=0))               #& normalize
1599#for fixed moments --> m=0 reflections
1600            fas0 = 0.
1601            fbs0 = 0.
1602            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
1603                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
1604                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]   
1605                           
1606                famq0 = np.sum(np.sum(fam0,axis=-2),axis=-2)        #Nref,Mxyz; sum ops & atoms
1607                fbmq0 = np.sum(np.sum(fbm0,axis=-2),axis=-2)
1608               
1609                fas0 = np.sum(famq0,axis=-1)**2-np.sum(eM.T*famq0,axis=-1)**2   #mag intensity calc F^2-(e.F)^2
1610                fbs0 = np.sum(fbmq0,axis=-1)**2-np.sum(eM.T*fbmq0,axis=-1)**2
1611#for modulated moments --> m != 0 reflections
1612            M = np.array(np.abs(H[3]),dtype=np.int)-1
1613       
1614            fam = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*cosm[i,:,:,nxs]-    \
1615                np.sign(H[3,i])*TmagA*sinm[i,:,:,nxs]),0.) for i in range(mRef)])
1616            fbm = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*sinm[i,:,:,nxs]+    \
1617                np.sign(H[3,i])*TmagA*cosm[i,:,:,nxs]),0.) for i in range(mRef)])
1618                       
1619            famq = np.sum(np.sum(fam/2.,axis=-2),axis=-2)      #Nref,Mxyz; sum ops & atoms
1620            fbmq = np.sum(np.sum(fbm/2.,axis=-2),axis=-2)
1621           
1622            fas = np.sum(famq,axis=-1)**2-np.sum(eM.T*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
1623            fbs = np.sum(fbmq,axis=-1)**2-np.sum(eM.T*fbmq,axis=-1)**2
1624           
1625            refl.T[10] = np.where(H[3],fas+fbs,fas0+fbs0)
1626            refl.T[11] = np.where(H[3],atan2d(fbs,fas),atan2d(fbs0,fas0))
1627
1628        else:
1629            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1630            if 'T' in calcControls[hfx+'histType']:
1631                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1632                fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1633            else:
1634                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1635                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1636            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1637            fbg = fb*GfpuA[0]+fa*GfpuA[1]
1638            fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1639            fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1640           
1641            refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2    #square of sums
1642            refl.T[11] = atan2d(fbs[0],fas[0])  #ignore f' & f"
1643        if 'P' not in calcControls[hfx+'histType']:
1644            refl.T[8] = np.copy(refl.T[10])               
1645        iBeg += blkSize
1646#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1647    return copy.deepcopy(refDict['RefList'])
1648
1649def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1650    '''
1651    Compute super structure factors for all h,k,l,m for phase - twins only
1652    puts the result, F^2, in each ref[8+im] in refList
1653    works on blocks of 32 reflections for speed
1654    input:
1655   
1656    :param dict refDict: where
1657        'RefList' list where each ref = h,k,l,m,it,d,...
1658        'FF' dict of form factors - filed in below
1659    :param np.array G:      reciprocal metric tensor
1660    :param str pfx:    phase id string
1661    :param dict SGData: space group info. dictionary output from SpcGroup
1662    :param dict calcControls:
1663    :param dict ParmDict:
1664
1665    '''
1666    phfx = pfx.split(':')[0]+hfx
1667    ast = np.sqrt(np.diag(G))
1668    Mast = twopisq*np.multiply.outer(ast,ast)   
1669    SGInv = SGData['SGInv']
1670    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1671    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1672    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1673    FFtables = calcControls['FFtables']
1674    BLtables = calcControls['BLtables']
1675    MFtables = calcControls['MFtables']
1676    Flack = 1.0
1677    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1678        Flack = 1.-2.*parmDict[phfx+'Flack']
1679    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])    #4D?
1680    TwDict = refDict.get('TwDict',{})           
1681    if 'S' in calcControls[hfx+'histType']:
1682        NTL = calcControls[phfx+'NTL']
1683        NM = calcControls[phfx+'TwinNMN']+1
1684        TwinLaw = calcControls[phfx+'TwinLaw']  #this'll have to be 4D also...
1685        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
1686        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
1687    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1688        GetAtomFXU(pfx,calcControls,parmDict)
1689    if not Xdata.size:          #no atoms in phase!
1690        return
1691    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1692    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
1693    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1694    FF = np.zeros(len(Tdata))
1695    if 'NC' in calcControls[hfx+'histType']:
1696        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1697    elif 'X' in calcControls[hfx+'histType']:
1698        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1699        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1700    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1701    bij = Mast*Uij
1702    blkSize = 32       #no. of reflections in a block
1703    nRef = refDict['RefList'].shape[0]
1704    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
1705        SQ = 1./(2.*refDict['RefList'].T[5])**2
1706        if 'N' in calcControls[hfx+'histType']:
1707            dat = G2el.getBLvalues(BLtables)
1708            refDict['FF']['El'] = list(dat.keys())
1709            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
1710            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
1711            for iel,El in enumerate(refDict['FF']['El']):
1712                if El in MFtables:
1713                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
1714        else:
1715            dat = G2el.getFFvalues(FFtables,0.)
1716            refDict['FF']['El'] = list(dat.keys())
1717            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
1718            for iel,El in enumerate(refDict['FF']['El']):
1719                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
1720#    time0 = time.time()
1721#reflection processing begins here - big arrays!
1722    iBeg = 0
1723    while iBeg < nRef:
1724        iFin = min(iBeg+blkSize,nRef)
1725        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
1726        H = refl[:,:4]                          #array(blkSize,4)
1727        H3 = refl[:,:3]
1728        HP = H[:,:3]+modQ[nxs,:]*H[:,3:]        #projected hklm to hkl
1729        HP = np.inner(HP,TwinLaw)             #array(blkSize,nTwins,4)
1730        H3 = np.inner(H3,TwinLaw)       
1731        TwMask = np.any(HP,axis=-1)
1732        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
1733            for ir in range(blkSize):
1734                iref = ir+iBeg
1735                if iref in TwDict:
1736                    for i in TwDict[iref]:
1737                        for n in range(NTL):
1738                            HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1739                            H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
1740            TwMask = np.any(HP,axis=-1)
1741        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
1742        SQfactor = 4.0*SQ*twopisq               #ditto prev.
1743        Uniq = np.inner(H,SSGMT)
1744        Uniq3 = np.inner(H3,SGMT)
1745        UniqP = np.inner(HP,SGMT)
1746        Phi = np.inner(H,SSGT)
1747        if SGInv:   #if centro - expand HKL sets
1748            Uniq = np.hstack((Uniq,-Uniq))
1749            Uniq3 = np.hstack((Uniq3,-Uniq3))
1750            Phi = np.hstack((Phi,-Phi))
1751            UniqP = np.hstack((UniqP,-UniqP))
1752        if 'T' in calcControls[hfx+'histType']:
1753            if 'P' in calcControls[hfx+'histType']:
1754                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
1755            else:
1756                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
1757            FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1758            FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0)
1759        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw))
1760        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1761        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0)
1762        phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs])
1763        sinp = np.sin(phase)
1764        cosp = np.cos(phase)
1765        biso = -SQfactor*Uisodata[:,nxs]
1766        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T
1767        HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1)  #use hklt proj to hkl
1768        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1769        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
1770        if 'T' in calcControls[hfx+'histType']:
1771            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
1772            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1773        else:
1774            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
1775            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
1776        GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
1777        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
1778        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1779        fas = np.sum(np.sum(fag,axis=-1),axis=-1)   #2 x refBlk; sum sym & atoms
1780        fbs = np.sum(np.sum(fbg,axis=-1),axis=-1)
1781        refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2                  #FcT from primary twin element
1782        refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+   \
1783            np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1)                 #Fc sum over twins
1784        refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
1785        iBeg += blkSize
1786#    print ('nRef %d time %.4f\r'%(nRef,time.time()-time0))
1787
1788def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
1789    '''
1790    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
1791    Only Fourier component are done analytically here
1792    input:
1793   
1794    :param dict refDict: where
1795        'RefList' list where each ref = h,k,l,m,it,d,...
1796        'FF' dict of form factors - filled in below
1797    :param int im: = 1 (could be eliminated)
1798    :param np.array G:      reciprocal metric tensor
1799    :param str hfx:    histogram id string
1800    :param str pfx:    phase id string
1801    :param dict SGData: space group info. dictionary output from SpcGroup
1802    :param dict SSGData: super space group info.
1803    :param dict calcControls:
1804    :param dict ParmDict:
1805   
1806    :returns: dict dFdvDict: dictionary of derivatives
1807    '''
1808    phfx = pfx.split(':')[0]+hfx
1809    ast = np.sqrt(np.diag(G))
1810    Mast = twopisq*np.multiply.outer(ast,ast)
1811    SGInv = SGData['SGInv']
1812    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1813    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1814    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1815    FFtables = calcControls['FFtables']
1816    BLtables = calcControls['BLtables']
1817    nRef = len(refDict['RefList'])
1818    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
1819        GetAtomFXU(pfx,calcControls,parmDict)
1820    if not Xdata.size:          #no atoms in phase!
1821        return {}
1822    mSize = len(Mdata)  #no. atoms
1823    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
1824    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
1825    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
1826    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
1827    FF = np.zeros(len(Tdata))
1828    if 'NC' in calcControls[hfx+'histType']:
1829        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
1830    elif 'X' in calcControls[hfx+'histType']:
1831        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1832        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1833    Uij = np.array(G2lat.U6toUij(Uijdata)).T
1834    bij = Mast*Uij
1835    if not len(refDict['FF']):
1836        if 'N' in calcControls[hfx+'histType']:
1837            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
1838        else:
1839            dat = G2el.getFFvalues(FFtables,0.)       
1840        refDict['FF']['El'] = list(dat.keys())
1841        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
1842    dFdvDict = {}
1843    dFdfr = np.zeros((nRef,mSize))
1844    dFdx = np.zeros((nRef,mSize,3))
1845    dFdui = np.zeros((nRef,mSize))
1846    dFdua = np.zeros((nRef,mSize,6))
1847    dFdbab = np.zeros((nRef,2))
1848    dFdfl = np.zeros((nRef))
1849    dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2))
1850    dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6))
1851    dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12))
1852    Flack = 1.0
1853    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
1854        Flack = 1.-2.*parmDict[phfx+'Flack']
1855    time0 = time.time()
1856    nRef = len(refDict['RefList'])/100
1857    for iref,refl in enumerate(refDict['RefList']):
1858        if 'T' in calcControls[hfx+'histType']:
1859            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
1860        H = np.array(refl[:4])
1861        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
1862        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
1863        SQfactor = 8.0*SQ*np.pi**2
1864        Bab = 0.0
1865        if phfx+'BabA' in parmDict:
1866            dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
1867            Bab = parmDict[phfx+'BabA']*dBabdA
1868        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
1869        FF = refDict['FF']['FF'][iref].T[Tindx]
1870        Uniq = np.inner(H,SSGMT)
1871        Phi = np.inner(H,SSGT)
1872        UniqP = np.inner(HP,SGMT)
1873        if SGInv:   #if centro - expand HKL sets
1874            Uniq = np.vstack((Uniq,-Uniq))
1875            Phi = np.hstack((Phi,-Phi))
1876            UniqP = np.vstack((UniqP,-UniqP))
1877        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
1878        sinp = np.sin(phase)
1879        cosp = np.cos(phase)
1880        occ = Mdata*Fdata/Uniq.shape[0]
1881        biso = -SQfactor*Uisodata[:,nxs]
1882        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T    #ops x atoms
1883        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
1884        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
1885        Hij = np.array([G2lat.UijtoU6(uij) for uij in Hij])                     #atoms x 6
1886        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
1887        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
1888        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
1889        fotp = FPP*Tcorr            #ops x atoms
1890        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
1891        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
1892        # GfpuA is 2 x ops x atoms
1893        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
1894        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms)
1895        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
1896        fag = fa*GfpuA[0]-fb*GfpuA[1]
1897        fbg = fb*GfpuA[0]+fa*GfpuA[1]
1898       
1899        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
1900        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
1901        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x ops x atoms
1902        fbx = np.array([fot*cosp,-fotp*sinp])
1903        fax = fax*GfpuA[0]-fbx*GfpuA[1]
1904        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
1905        #sum below is over Uniq
1906        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
1907        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
1908        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
1909        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
1910        dfadui = np.sum(-SQfactor*fag,axis=1)
1911        dfbdui = np.sum(-SQfactor*fbg,axis=1)
1912        dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)  #2 x nAtom x 3xyz; sum nOps
1913        dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
1914        dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)             #2 x nAtom x 6Uij; sum nOps
1915        dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)         #these are correct also for twins above
1916        # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps
1917        dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1918        dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1)
1919        dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1920        dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
1921        dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
1922        dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
1923        if not SGData['SGInv']:   #Flack derivative
1924            dfadfl = np.sum(-FPP*Tcorr*sinp)
1925            dfbdfl = np.sum(FPP*Tcorr*cosp)
1926        else:
1927            dfadfl = 1.0
1928            dfbdfl = 1.0
1929        SA = fas[0]+fas[1]      #float = A+A'
1930        SB = fbs[0]+fbs[1]      #float = B+B'
1931        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro?
1932            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1933            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
1934                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1935            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
1936                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1937            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
1938                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1939            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
1940                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1941            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
1942                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
1943            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
1944                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
1945            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
1946                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
1947        else:                       #OK, I think
1948            dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
1949            dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
1950            dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
1951            dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
1952            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
1953                           
1954            dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
1955            dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
1956            dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
1957        if phfx+'BabA' in parmDict:
1958            dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
1959                2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
1960        #loop over atoms - each dict entry is list of derivatives for all the reflections
1961        if not iref%100 :
1962            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
1963    for i in range(len(Mdata)):     #loop over atoms
1964        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1965        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1966        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1967        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1968        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1969        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1970        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1971        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1972        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1973        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1974        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1975        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
1976            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
1977            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
1978        nx = 0
1979        if waveTypes[i] in ['Block','ZigZag']:
1980            nx = 1 
1981        for j in range(XSSdata.shape[1]-nx):       #loop over waves
1982            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
1983            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
1984            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
1985            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
1986            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
1987            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
1988        for j in range(USSdata.shape[1]):       #loop over waves
1989            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
1990            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
1991            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
1992            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
1993            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
1994            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
1995            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
1996            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
1997            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
1998            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
1999            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2000            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2001           
2002    dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
2003    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2004    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2005    return dFdvDict
2006
2007def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2008    '''
2009    Compute super structure factor derivatives for all h,k,l,m for phase - no twins
2010    input:
2011   
2012    :param dict refDict: where
2013        'RefList' list where each ref = h,k,l,m,it,d,...
2014        'FF' dict of form factors - filled in below
2015    :param int im: = 1 (could be eliminated)
2016    :param np.array G:      reciprocal metric tensor
2017    :param str hfx:    histogram id string
2018    :param str pfx:    phase id string
2019    :param dict SGData: space group info. dictionary output from SpcGroup
2020    :param dict SSGData: super space group info.
2021    :param dict calcControls:
2022    :param dict ParmDict:
2023   
2024    :returns: dict dFdvDict: dictionary of derivatives
2025    '''
2026
2027    trefDict = copy.deepcopy(refDict)
2028    dM = 1.e-4
2029    dFdvDict = {}
2030    for parm in parmDict:
2031        if parm == '0':
2032            continue
2033        if parm.split(':')[2] in ['Tmin','Tmax','Xmax','Ymax','Zmax','Fzero','Fwid',
2034            'MXsin','MXcos','MYsin','MYcos','MZsin','MZcos','AMx','AMy','AMz',]:
2035            parmDict[parm] += dM
2036            prefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2037            parmDict[parm] -= 2*dM
2038            mrefList = SStructureFactor(trefDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
2039            parmDict[parm] += dM
2040            dFdvDict[parm] = (prefList[:,9+im]-mrefList[:,9+im])/(2.*dM)
2041    return dFdvDict
2042   
2043def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
2044    'Needs a doc string'
2045    phfx = pfx.split(':')[0]+hfx
2046    ast = np.sqrt(np.diag(G))
2047    Mast = twopisq*np.multiply.outer(ast,ast)
2048    SGInv = SGData['SGInv']
2049    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2050    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2051    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2052    FFtables = calcControls['FFtables']
2053    BLtables = calcControls['BLtables']
2054    TwinLaw = np.array([[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]],])
2055    TwDict = refDict.get('TwDict',{})           
2056    if 'S' in calcControls[hfx+'histType']:
2057        NTL = calcControls[phfx+'NTL']
2058        NM = calcControls[phfx+'TwinNMN']+1
2059        TwinLaw = calcControls[phfx+'TwinLaw']
2060        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
2061    nTwin = len(TwinLaw)       
2062    nRef = len(refDict['RefList'])
2063    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
2064        GetAtomFXU(pfx,calcControls,parmDict)
2065    if not Xdata.size:          #no atoms in phase!
2066        return {}
2067    mSize = len(Mdata)  #no. atoms
2068    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
2069    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
2070    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
2071    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2072    FF = np.zeros(len(Tdata))
2073    if 'NC' in calcControls[hfx+'histType']:
2074        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
2075    elif 'X' in calcControls[hfx+'histType']:
2076        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2077        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2078    Uij = np.array(G2lat.U6toUij(Uijdata)).T
2079    bij = Mast*Uij
2080    if not len(refDict['FF']):
2081        if 'N' in calcControls[hfx+'histType']:
2082            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
2083        else:
2084            dat = G2el.getFFvalues(FFtables,0.)       
2085        refDict['FF']['El'] = list(dat.keys())
2086        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
2087    dFdvDict = {}
2088    dFdfr = np.zeros((nRef,nTwin,mSize))
2089    dFdx = np.zeros((nRef,nTwin,mSize,3))
2090    dFdui = np.zeros((nRef,nTwin,mSize))
2091    dFdua = np.zeros((nRef,nTwin,mSize,6))
2092    dFdbab = np.zeros((nRef,nTwin,2))
2093    dFdtw = np.zeros((nRef,nTwin))
2094    dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
2095    dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
2096    dFdGz = np.zeros((nRef,nTwin,mSize,5))
2097    dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
2098    Flack = 1.0
2099    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
2100        Flack = 1.-2.*parmDict[phfx+'Flack']
2101    time0 = time.time()
2102    nRef = len(refDict['RefList'])/100
2103    for iref,refl in enumerate(refDict['RefList']):
2104        if 'T' in calcControls[hfx+'histType']:
2105            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im])
2106        H = np.array(refl[:4])
2107        HP = H[:3]+modQ*H[3:]            #projected hklm to hkl
2108        H = np.inner(H.T,TwinLaw)   #maybe array(4,nTwins) or (4)
2109        TwMask = np.any(H,axis=-1)
2110        if TwinLaw.shape[0] > 1 and TwDict:
2111            if iref in TwDict:
2112                for i in TwDict[iref]:
2113                    for n in range(NTL):
2114                        H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
2115            TwMask = np.any(H,axis=-1)
2116        SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2
2117        SQfactor = 8.0*SQ*np.pi**2
2118        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2119        Bab = parmDict[phfx+'BabA']*dBabdA
2120        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
2121        FF = refDict['FF']['FF'][iref].T[Tindx]
2122        Uniq = np.inner(H,SSGMT)
2123        Phi = np.inner(H,SSGT)
2124        UniqP = np.inner(HP,SGMT)
2125        if SGInv:   #if centro - expand HKL sets
2126            Uniq = np.vstack((Uniq,-Uniq))
2127            Phi = np.hstack((Phi,-Phi))
2128            UniqP = np.vstack((UniqP,-UniqP))
2129        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
2130        sinp = np.sin(phase)
2131        cosp = np.cos(phase)
2132        occ = Mdata*Fdata/Uniq.shape[0]
2133        biso = -SQfactor*Uisodata[:,nxs]
2134        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
2135        HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1)  #ops x atoms
2136        Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3
2137        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(nTwin,-1,6)))
2138        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
2139        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
2140        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
2141        fotp = FPP*Tcorr            #ops x atoms
2142        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
2143        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
2144        # GfpuA is 2 x ops x atoms
2145        # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts
2146        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
2147        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])  #or array(2,nEqv,nAtoms)
2148        fag = fa*GfpuA[0]-fb*GfpuA[1]
2149        fbg = fb*GfpuA[0]+fa*GfpuA[1]
2150       
2151        fas = np.sum(np.sum(fag,axis=1),axis=1)     # 2 x twin
2152        fbs = np.sum(np.sum(fbg,axis=1),axis=1)
2153        fax = np.array([-fot*sinp,-fotp*cosp])   #positions; 2 x twin x ops x atoms
2154        fbx = np.array([fot*cosp,-fotp*sinp])
2155        fax = fax*GfpuA[0]-fbx*GfpuA[1]
2156        fbx = fbx*GfpuA[0]+fax*GfpuA[1]
2157        #sum below is over Uniq
2158        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
2159        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
2160        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
2161        dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)
2162        dfadui = np.sum(-SQfactor*fag,axis=1)
2163        dfbdui = np.sum(-SQfactor*fbg,axis=1)
2164        dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2165        dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
2166        dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2167        dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
2168        # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12)
2169        dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2170        dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1)
2171        dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2172        dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
2173        dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2174        dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1)
2175        dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2176        dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1)
2177#        GSASIIpath.IPyBreak()
2178        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
2179        SA = fas[0]+fas[1]      #float = A+A' (might be array[nTwin])
2180        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
2181        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)]
2182        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)]
2183        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)]
2184        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)]
2185        dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2
2186
2187        dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)]
2188        dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)]
2189        dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)]
2190        dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)]               
2191#            GSASIIpath.IPyBreak()
2192        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
2193            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
2194        #loop over atoms - each dict entry is list of derivatives for all the reflections
2195        if not iref%100 :
2196            print (' %d derivative time %.4f\r'%(iref,time.time()-time0),end='')
2197    for i in range(len(Mdata)):     #loop over atoms
2198        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2199        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2200        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2201        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2202        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2203        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2204        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2205        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2206        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2207        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2208        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2209        for j in range(FSSdata.shape[1]):        #loop over waves Fzero & Fwid?
2210            dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i]
2211            dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i]
2212        nx = 0
2213        if waveTypes[i] in ['Block','ZigZag']:
2214            nx = 1 
2215            dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i]   #ZigZag/Block waves (if any)
2216            dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i]
2217            dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i]
2218            dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i]
2219            dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i]
2220        for j in range(XSSdata.shape[1]-nx):       #loop over waves
2221            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i]
2222            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i]
2223            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i]
2224            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i]
2225            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i]
2226            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i]
2227        for j in range(USSdata.shape[1]):       #loop over waves
2228            dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
2229            dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i]
2230            dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i]
2231            dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i]
2232            dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i]
2233            dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i]
2234            dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i]
2235            dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i]
2236            dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i]
2237            dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i]
2238            dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i]
2239            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
2240           
2241#        GSASIIpath.IPyBreak()
2242    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
2243    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
2244    return dFdvDict
2245   
2246def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
2247    ''' Single crystal extinction function; returns extinction & derivative
2248    '''
2249    extCor = 1.0
2250    dervDict = {}
2251    dervCor = 1.0
2252    if calcControls[phfx+'EType'] != 'None':
2253        SQ = 1/(4.*ref[4+im]**2)
2254        if 'C' in parmDict[hfx+'Type']:           
2255            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
2256        else:   #'T'
2257            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
2258        if 'SXC' in parmDict[hfx+'Type']:
2259            AV = 7.9406e5/parmDict[pfx+'Vol']**2
2260            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2261            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
2262            PLZ = AV*P12*ref[9+im]*parmDict[hfx+'Lam']**2
2263        elif 'SNT' in parmDict[hfx+'Type']:
2264            AV = 1.e7/parmDict[pfx+'Vol']**2
2265            PL = SQ
2266            PLZ = AV*ref[9+im]*ref[12+im]**2
2267        elif 'SNC' in parmDict[hfx+'Type']:
2268            AV = 1.e7/parmDict[pfx+'Vol']**2
2269            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
2270            PLZ = AV*ref[9+im]*parmDict[hfx+'Lam']**2
2271           
2272        if 'Primary' in calcControls[phfx+'EType']:
2273            PLZ *= 1.5
2274        else:
2275            if 'C' in parmDict[hfx+'Type']:
2276                PLZ *= calcControls[phfx+'Tbar']
2277            else: #'T'
2278                PLZ *= ref[13+im]      #t-bar
2279        if 'Primary' in calcControls[phfx+'EType']:
2280            PLZ *= 1.5
2281            PSIG = parmDict[phfx+'Ep']
2282        elif 'I & II' in calcControls[phfx+'EType']:
2283            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
2284        elif 'Type II' in calcControls[phfx+'EType']:
2285            PSIG = parmDict[phfx+'Es']
2286        else:       # 'Secondary Type I'
2287            PSIG = parmDict[phfx+'Eg']/PL
2288           
2289        AG = 0.58+0.48*cos2T+0.24*cos2T**2
2290        AL = 0.025+0.285*cos2T
2291        BG = 0.02-0.025*cos2T
2292        BL = 0.15-0.2*(0.75-cos2T)**2
2293        if cos2T < 0.:
2294            BL = -0.45*cos2T
2295        CG = 2.
2296        CL = 2.
2297        PF = PLZ*PSIG
2298       
2299        if 'Gaussian' in calcControls[phfx+'EApprox']:
2300            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
2301            extCor = np.sqrt(PF4)
2302            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
2303        else:
2304            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
2305            extCor = np.sqrt(PF4)
2306            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
2307
2308        dervCor = (1.+PF)*PF3   #extinction corr for other derivatives
2309        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
2310            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
2311        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
2312            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
2313        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
2314            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
2315               
2316    return 1./extCor,dervDict,dervCor
2317   
2318def Dict2Values(parmdict, varylist):
2319    '''Use before call to leastsq to setup list of values for the parameters
2320    in parmdict, as selected by key in varylist'''
2321    return [parmdict[key] for key in varylist] 
2322   
2323def Values2Dict(parmdict, varylist, values):
2324    ''' Use after call to leastsq to update the parameter dictionary with
2325    values corresponding to keys in varylist'''
2326    parmdict.update(zip(varylist,values))
2327   
2328def GetNewCellParms(parmDict,varyList):
2329    '''Compute unit cell tensor terms from varied Aij and Dij values.
2330    Terms are included in the dict only if Aij or Dij is varied.
2331    '''
2332    newCellDict = {}
2333    Anames = ['A'+str(i) for i in range(6)]
2334    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
2335    for item in varyList:
2336        keys = item.split(':')
2337        if keys[2] in Ddict:
2338            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
2339            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
2340            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2341    return newCellDict          # is e.g. {'0::D11':A0-D11}
2342   
2343def ApplyXYZshifts(parmDict,varyList):
2344    '''
2345    takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2346   
2347    :param dict parmDict: parameter dictionary
2348    :param list varyList: list of variables (not used!)
2349    :returns: newAtomDict - dictionary of new atomic coordinate names & values; key is parameter shift name
2350
2351    '''
2352    newAtomDict = {}
2353    for item in parmDict:
2354        if 'dA' in item:
2355            parm = ''.join(item.split('d'))
2356            parmDict[parm] += parmDict[item]
2357            newAtomDict[item] = [parm,parmDict[parm]]
2358    return newAtomDict
2359   
2360def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2361    'Spherical harmonics texture'
2362    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2363    if 'T' in calcControls[hfx+'histType']:
2364        tth = parmDict[hfx+'2-theta']
2365    else:
2366        tth = refl[5+im]
2367    odfCor = 1.0
2368    H = refl[:3]
2369    cell = G2lat.Gmat2cell(g)
2370    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2371    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2372    phi,beta = G2lat.CrsAng(H,cell,SGData)
2373    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2374    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2375    for item in SHnames:
2376        L,M,N = eval(item.strip('C'))
2377        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2378        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2379        Lnorm = G2lat.Lnorm(L)
2380        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2381    return odfCor
2382   
2383def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
2384    'Spherical harmonics texture derivatives'
2385    if 'T' in calcControls[hfx+'histType']:
2386        tth = parmDict[hfx+'2-theta']
2387    else:
2388        tth = refl[5+im]
2389    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2390    odfCor = 1.0
2391    dFdODF = {}
2392    dFdSA = [0,0,0]
2393    H = refl[:3]
2394    cell = G2lat.Gmat2cell(g)
2395    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2396    Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2397    phi,beta = G2lat.CrsAng(H,cell,SGData)
2398    psi,gam,dPSdA,dGMdA = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup)
2399    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2400    for item in SHnames:
2401        L,M,N = eval(item.strip('C'))
2402        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2403        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2404        Lnorm = G2lat.Lnorm(L)
2405        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2406        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2407        for i in range(3):
2408            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2409    return odfCor,dFdODF,dFdSA
2410   
2411def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2412    'spherical harmonics preferred orientation (cylindrical symmetry only)'
2413    if 'T' in calcControls[hfx+'histType']:
2414        tth = parmDict[hfx+'2-theta']
2415    else:
2416        tth = refl[5+im]
2417    odfCor = 1.0
2418    H = refl[:3]
2419    cell = G2lat.Gmat2cell(g)
2420    Sangls = [0.,0.,0.]
2421    if 'Bragg' in calcControls[hfx+'instType']:
2422        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2423        IFCoup = True
2424    else:
2425        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2426        IFCoup = False
2427    phi,beta = G2lat.CrsAng(H,cell,SGData)
2428    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2429    SHnames = calcControls[phfx+'SHnames']
2430    for item in SHnames:
2431        L,N = eval(item.strip('C'))
2432        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2433        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2434        Lnorm = G2lat.Lnorm(L)
2435        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2436    return np.squeeze(odfCor)
2437   
2438def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
2439    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
2440    if 'T' in calcControls[hfx+'histType']:
2441        tth = parmDict[hfx+'2-theta']
2442    else:
2443        tth = refl[5+im]
2444    odfCor = 1.0
2445    dFdODF = {}
2446    H = refl[:3]
2447    cell = G2lat.Gmat2cell(g)
2448    Sangls = [0.,0.,0.]
2449    if 'Bragg' in calcControls[hfx+'instType']:
2450        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2451        IFCoup = True
2452    else:
2453        Gangls = [parmDict[hfx+'Phi'],parmDict[hfx+'Chi'],parmDict[hfx+'Omega'],parmDict[hfx+'Azimuth']]
2454        IFCoup = False
2455    phi,beta = G2lat.CrsAng(H,cell,SGData)
2456    psi,gam,x,x = G2lat.SamAng(tth/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2457    SHnames = calcControls[phfx+'SHnames']
2458    for item in SHnames:
2459        L,N = eval(item.strip('C'))
2460        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2461        Ksl,x,x = G2lat.GetKsl(L,0,'0',psi,gam)
2462        Lnorm = G2lat.Lnorm(L)
2463        odfCor += parmDict[phfx+item]*Lnorm*Kcl*Ksl
2464        dFdODF[phfx+item] = Kcl*Ksl*Lnorm
2465    return odfCor,dFdODF
2466   
2467def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2468    'March-Dollase preferred orientation correction'
2469    POcorr = 1.0
2470    MD = parmDict[phfx+'MD']
2471    if MD != 1.0:
2472        MDAxis = calcControls[phfx+'MDAxis']
2473        sumMD = 0
2474        for H in uniq:           
2475            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2476            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2477            sumMD += A**3
2478        POcorr = sumMD/len(uniq)
2479    return POcorr
2480   
2481def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
2482    'Needs a doc string'
2483    POcorr = 1.0
2484    POderv = {}
2485    if calcControls[phfx+'poType'] == 'MD':
2486        MD = parmDict[phfx+'MD']
2487        MDAxis = calcControls[phfx+'MDAxis']
2488        sumMD = 0
2489        sumdMD = 0
2490        for H in uniq:           
2491            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2492            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2493            sumMD += A**3
2494            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2495        POcorr = sumMD/len(uniq)
2496        POderv[phfx+'MD'] = sumdMD/len(uniq)
2497    else:   #spherical harmonics
2498        if calcControls[phfx+'SHord']:
2499            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2500    return POcorr,POderv
2501   
2502def GetAbsorb(refl,im,hfx,calcControls,parmDict):
2503    'Needs a doc string'
2504    if 'Debye' in calcControls[hfx+'instType']:
2505        if 'T' in calcControls[hfx+'histType']:
2506            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2507        else:
2508            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2509    else:
2510        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
2511   
2512def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
2513    'Needs a doc string'
2514    if 'Debye' in calcControls[hfx+'instType']:
2515        if 'T' in calcControls[hfx+'histType']:
2516            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],abs(parmDict[hfx+'2-theta']),0,0)
2517        else:
2518            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
2519    else:
2520        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
2521       
2522def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2523    'Needs a doc string'
2524    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2525    pi2 = np.sqrt(2./np.pi)
2526    if 'T' in calcControls[hfx+'histType']:
2527        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2528        wave = refl[14+im]
2529    else:   #'C'W
2530        sth2 = sind(refl[5+im]/2.)**2
2531        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2532    c2th = 1.-2.0*sth2
2533    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2534    if 'X' in calcControls[hfx+'histType']:
2535        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2536    xfac = flv2*parmDict[phfx+'Extinction']
2537    exb = 1.0
2538    if xfac > -1.:
2539        exb = 1./np.sqrt(1.+xfac)
2540    exl = 1.0
2541    if 0 < xfac <= 1.:
2542        xn = np.array([xfac**(i+1) for i in range(6)])
2543        exl += np.sum(xn*coef)
2544    elif xfac > 1.:
2545        xfac2 = 1./np.sqrt(xfac)
2546        exl = pi2*(1.-0.125/xfac)*xfac2
2547    return exb*sth2+exl*(1.-sth2)
2548   
2549def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
2550    'Needs a doc string'
2551    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
2552    pi2 = np.sqrt(2./np.pi)
2553    if 'T' in calcControls[hfx+'histType']:
2554        sth2 = sind(abs(parmDict[hfx+'2-theta'])/2.)**2
2555        wave = refl[14+im]
2556    else:   #'C'W
2557        sth2 = sind(refl[5+im]/2.)**2
2558        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
2559    c2th = 1.-2.0*sth2
2560    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
2561    if 'X' in calcControls[hfx+'histType']:
2562        flv2 *= 0.079411*(1.0+c2th**2)/2.0
2563    xfac = flv2*parmDict[phfx+'Extinction']
2564    dbde = -500.*flv2
2565    if xfac > -1.:
2566        dbde = -0.5*flv2/np.sqrt(1.+xfac)**3
2567    dlde = 0.
2568    if 0 < xfac <= 1.:
2569        xn = np.array([i*flv2*xfac**i for i in [1,2,3,4,5,6]])
2570        dlde = np.sum(xn*coef)/xfac
2571    elif xfac > 1.:
2572        xfac2 = 1./np.sqrt(xfac)
2573        dlde = 0.5*flv2*pi2*xfac2*(-1./xfac+0.375/xfac**2)
2574       
2575    return dbde*sth2+dlde*(1.-sth2)
2576   
2577def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2578    'Needs a doc string'    #need powder extinction!
2579    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
2580    if 'X' in parmDict[hfx+'Type']:
2581        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
2582    POcorr = 1.0
2583    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture - takes precidence
2584        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2585    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
2586        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
2587    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
2588        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
2589    Icorr *= POcorr
2590    AbsCorr = 1.0
2591    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
2592    Icorr *= AbsCorr
2593    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
2594    Icorr *= ExtCorr
2595    return Icorr,POcorr,AbsCorr,ExtCorr
2596   
2597def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2598    'Needs a doc string'    #need powder extinction derivs!
2599    dIdsh = 1./parmDict[hfx+'Scale']
2600    dIdsp = 1./parmDict[phfx+'Scale']
2601    if 'X' in parmDict[hfx+'Type']:
2602        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
2603        dIdPola /= pola
2604    else:       #'N'
2605        dIdPola = 0.0
2606    dFdODF = {}
2607    dFdSA = [0,0,0]
2608    dIdPO = {}
2609    if pfx+'SHorder' in parmDict:
2610        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
2611        for iSH in dFdODF:
2612            dFdODF[iSH] /= odfCor
2613        for i in range(3):
2614            dFdSA[i] /= odfCor
2615    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
2616        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
2617        for iPO in dIdPO:
2618            dIdPO[iPO] /= POcorr
2619    if 'T' in parmDict[hfx+'Type']:
2620        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
2621        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
2622    else:
2623        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
2624        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
2625    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
2626       
2627def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2628    'Needs a doc string'
2629    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
2630        costh = cosd(refl[5+im]/2.)
2631        #crystallite size
2632        if calcControls[phfx+'SizeType'] == 'isotropic':
2633            Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
2634        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2635            H = np.array(refl[:3])
2636            P = np.array(calcControls[phfx+'SizeAxis'])
2637            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2638            Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
2639            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2640        else:           #ellipsoidal crystallites
2641            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2642            H = np.array(refl[:3])
2643            lenR = G2pwd.ellipseSize(H,Sij,GB)
2644            Sgam = 1.8*wave/(np.pi*costh*lenR)
2645        #microstrain               
2646        if calcControls[phfx+'MustrainType'] == 'isotropic':
2647            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2648        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2649            H = np.array(refl[:3])
2650            P = np.array(calcControls[phfx+'MustrainAxis'])
2651            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2652            Si = parmDict[phfx+'Mustrain;i']
2653            Sa = parmDict[phfx+'Mustrain;a']
2654            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2655        else:       #generalized - P.W. Stephens model
2656            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2657            Sum = 0
2658            for i,strm in enumerate(Strms):
2659                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2660            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
2661    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
2662        #crystallite size
2663        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2664            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2665        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2666            H = np.array(refl[:3])
2667            P = np.array(calcControls[phfx+'SizeAxis'])
2668            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2669            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
2670            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
2671        else:           #ellipsoidal crystallites   #OK
2672            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2673            H = np.array(refl[:3])
2674            lenR = G2pwd.ellipseSize(H,Sij,GB)
2675            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
2676        #microstrain               
2677        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
2678            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2679        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
2680            H = np.array(refl[:3])
2681            P = np.array(calcControls[phfx+'MustrainAxis'])
2682            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2683            Si = parmDict[phfx+'Mustrain;i']
2684            Sa = parmDict[phfx+'Mustrain;a']
2685            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2686        else:       #generalized - P.W. Stephens model  OK
2687            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2688            Sum = 0
2689            for i,strm in enumerate(Strms):
2690                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2691            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
2692           
2693    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
2694    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
2695    sig /= ateln2
2696    return sig,gam
2697       
2698def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
2699    'Needs a doc string'
2700    gamDict = {}
2701    sigDict = {}
2702    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
2703        costh = cosd(refl[5+im]/2.)
2704        tanth = tand(refl[5+im]/2.)
2705        #crystallite size derivatives
2706        if calcControls[phfx+'SizeType'] == 'isotropic':
2707            Sgam = 1.8*wave/(np.pi*costh*parmDict[phfx+'Size;i'])
2708            gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh*parmDict[phfx+'Size;i']**2)
2709            sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
2710        elif calcControls[phfx+'SizeType'] == 'uniaxial':
2711            H = np.array(refl[:3])
2712            P = np.array(calcControls[phfx+'SizeAxis'])
2713            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2714            Si = parmDict[phfx+'Size;i']
2715            Sa = parmDict[phfx+'Size;a']
2716            gami = 1.8*wave/(costh*np.pi*Si*Sa)
2717            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2718            Sgam = gami*sqtrm
2719            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2720            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2721            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2722            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2723            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2724            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2725        else:           #ellipsoidal crystallites
2726            const = 1.8*wave/(np.pi*costh)
2727            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2728            H = np.array(refl[:3])
2729            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2730            Sgam = const/lenR
2731            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2732                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2733                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2734        gamDict[phfx+'Size;mx'] = Sgam
2735        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
2736               
2737        #microstrain derivatives               
2738        if calcControls[phfx+'MustrainType'] == 'isotropic':
2739            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
2740            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
2741            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
2742        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2743            H = np.array(refl[:3])
2744            P = np.array(calcControls[phfx+'MustrainAxis'])
2745            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2746            Si = parmDict[phfx+'Mustrain;i']
2747            Sa = parmDict[phfx+'Mustrain;a']
2748            gami = 0.018*Si*Sa*tanth/np.pi
2749            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2750            Mgam = gami/sqtrm
2751            dsi = -gami*Si*cosP**2/sqtrm**3
2752            dsa = -gami*Sa*sinP**2/sqtrm**3
2753            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2754            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2755            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2756            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2757        else:       #generalized - P.W. Stephens model
2758            const = 0.018*refl[4+im]**2*tanth/np.pi
2759            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2760            Sum = 0
2761            for i,strm in enumerate(Strms):
2762                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2763                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2764                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2765            Mgam = const*np.sqrt(Sum)
2766            for i in range(len(Strms)):
2767                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2768                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2
2769        gamDict[phfx+'Mustrain;mx'] = Mgam
2770        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2771    else:   #'T'OF - All checked & OK
2772        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
2773            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
2774            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
2775            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
2776        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
2777            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2778            H = np.array(refl[:3])
2779            P = np.array(calcControls[phfx+'SizeAxis'])
2780            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2781            Si = parmDict[phfx+'Size;i']
2782            Sa = parmDict[phfx+'Size;a']
2783            gami = const/(Si*Sa)
2784            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2785            Sgam = gami*sqtrm
2786            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
2787            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
2788            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
2789            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
2790            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2791            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2792        else:           #OK  ellipsoidal crystallites
2793            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
2794            Sij =[parmDict[phfx+'Size;%d'%(i)] for i in range(6)]
2795            H = np.array(refl[:3])
2796            lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2797            Sgam = const/lenR
2798            for i,item in enumerate([phfx+'Size;%d'%(j) for j in range(6)]):
2799                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
2800                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
2801        gamDict[phfx+'Size;mx'] = Sgam  #OK
2802        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
2803               
2804        #microstrain derivatives               
2805        if calcControls[phfx+'MustrainType'] == 'isotropic':
2806            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
2807            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
2808            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
2809        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2810            H = np.array(refl[:3])
2811            P = np.array(calcControls[phfx+'MustrainAxis'])
2812            cosP,sinP = G2lat.CosSinAngle(H,P,G)
2813            Si = parmDict[phfx+'Mustrain;i']
2814            Sa = parmDict[phfx+'Mustrain;a']
2815            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
2816            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2817            Mgam = gami/sqtrm
2818            dsi = -gami*Si*cosP**2/sqtrm**3
2819            dsa = -gami*Sa*sinP**2/sqtrm**3
2820            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
2821            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
2822            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
2823            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
2824        else:       #generalized - P.W. Stephens model OK
2825            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
2826            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
2827            Sum = 0
2828            for i,strm in enumerate(Strms):
2829                Sum += parmDict[phfx+'Mustrain;'+str(i)]*strm
2830                gamDict[phfx+'Mustrain;'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
2831                sigDict[phfx+'Mustrain;'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
2832            Mgam = const*np.sqrt(Sum)
2833            for i in range(len(Strms)):
2834                gamDict[phfx+'Mustrain;'+str(i)] *= Mgam/Sum
2835                sigDict[phfx+'Mustrain;'+str(i)] *= const**2/ateln2       
2836        gamDict[phfx+'Mustrain;mx'] = Mgam
2837        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
2838       
2839    return sigDict,gamDict
2840       
2841def GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2842    'Needs a doc string'
2843    if im:
2844        h,k,l,m = refl[:4]
2845        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2846        d = 1./np.sqrt(G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec))
2847    else:
2848        h,k,l = refl[:3]
2849        d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
2850    refl[4+im] = d
2851    if 'C' in calcControls[hfx+'histType']:
2852        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2853        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2854        if 'Bragg' in calcControls[hfx+'instType']:
2855            pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2856                parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2857        else:               #Debye-Scherrer - simple but maybe not right
2858            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2859    elif 'T' in calcControls[hfx+'histType']:
2860        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
2861        #do I need sample position effects - maybe?
2862    return pos
2863
2864def GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict):
2865    'Needs a doc string'
2866    dpr = 180./np.pi
2867    if im:
2868        h,k,l,m = refl[:4]
2869        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
2870        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
2871        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
2872    else:
2873        m = 0
2874        h,k,l = refl[:3]       
2875        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2876    dst = np.sqrt(dstsq)
2877    dsp = 1./dst
2878    if 'C' in calcControls[hfx+'histType']:
2879        pos = refl[5+im]-parmDict[hfx+'Zero']
2880        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2881        dpdw = const*dst
2882        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*const*wave/(2.0*dst)
2883        dpdZ = 1.0
2884        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2885            2*l*A[2]+h*A[4]+k*A[5]])*m*const*wave/(2.0*dst)
2886        shft = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2887        if 'Bragg' in calcControls[hfx+'instType']:
2888            dpdSh = -4.*shft*cosd(pos/2.0)
2889            dpdTr = -shft*sind(pos)*100.0
2890            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.,dpdV
2891        else:               #Debye-Scherrer - simple but maybe not right
2892            dpdXd = -shft*cosd(pos)
2893            dpdYd = -shft*sind(pos)
2894            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
2895    elif 'T' in calcControls[hfx+'histType']:
2896        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
2897        dpdZ = 1.0
2898        dpdDC = dsp
2899        dpdDA = dsp**2
2900        dpdDB = 1./dsp
2901        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
2902            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
2903        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
2904           
2905def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2906    'Needs a doc string'
2907    laue = SGData['SGLaue']
2908    uniq = SGData['SGUniq']
2909    h,k,l = refl[:3]
2910    if laue in ['m3','m3m']:
2911        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2912            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2913    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2914        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2915    elif laue in ['3R','3mR']:
2916        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2917    elif laue in ['4/m','4/mmm']:
2918        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2919    elif laue in ['mmm']:
2920        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2921    elif laue in ['2/m']:
2922        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2923        if uniq == 'a':
2924            Dij += parmDict[phfx+'D23']*k*l
2925        elif uniq == 'b':
2926            Dij += parmDict[phfx+'D13']*h*l
2927        elif uniq == 'c':
2928            Dij += parmDict[phfx+'D12']*h*k
2929    else:
2930        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2931            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2932    if 'C' in calcControls[hfx+'histType']:
2933        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2934    else:
2935        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
2936           
2937def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
2938    'Needs a doc string'
2939    laue = SGData['SGLaue']
2940    uniq = SGData['SGUniq']
2941    h,k,l = refl[:3]
2942    if laue in ['m3','m3m']:
2943        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2944            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2945    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2946        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2947    elif laue in ['3R','3mR']:
2948        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2949    elif laue in ['4/m','4/mmm']:
2950        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2951    elif laue in ['mmm']:
2952        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2953    elif laue in ['2/m']:
2954        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2955        if uniq == 'a':
2956            dDijDict[phfx+'D23'] = k*l
2957        elif uniq == 'b':
2958            dDijDict[phfx+'D13'] = h*l
2959        elif uniq == 'c':
2960            dDijDict[phfx+'D12'] = h*k
2961    else:
2962        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2963            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2964    if 'C' in calcControls[hfx+'histType']:
2965        for item in dDijDict:
2966            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
2967    else:
2968        for item in dDijDict:
2969            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
2970    return dDijDict
2971   
2972def GetDij(phfx,SGData,parmDict):
2973    HSvals = [parmDict[phfx+name] for name in G2spc.HStrainNames(SGData)]
2974    return G2spc.HStrainVals(HSvals,SGData)
2975               
2976def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2977    '''Compute the observed structure factors for Powder histograms and store in reflection array
2978    Multiprocessing support added
2979    '''
2980    if GSASIIpath.GetConfigValue('Show_timing',False):
2981        starttime = time.time() #; print 'start GetFobsSq'
2982    histoList = list(Histograms.keys())
2983    histoList.sort()
2984    Ka2 = shl = lamRatio = kRatio = None
2985    for histogram in histoList:
2986        if 'PWDR' in histogram[:4]:
2987            Histogram = Histograms[histogram]
2988            hId = Histogram['hId']
2989            hfx = ':%d:'%(hId)
2990            Limits = calcControls[hfx+'Limits']
2991            if 'C' in calcControls[hfx+'histType']:
2992                shl = max(parmDict[hfx+'SH/L'],0.0005)
2993                Ka2 = False
2994                kRatio = 0.0
2995                if hfx+'Lam1' in list(parmDict.keys()):
2996                    Ka2 = True
2997                    lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2998                    kRatio = parmDict[hfx+'I(L2)/I(L1)']
2999            x,y,w,yc,yb,yd = Histogram['Data']
3000            xMask = ma.getmaskarray(x)
3001            xB = np.searchsorted(x,Limits[0])
3002            xF = np.searchsorted(x,Limits[1])
3003            ymb = np.array(y-yb)
3004            ymb = np.where(ymb,ymb,1.0)
3005            ycmb = np.array(yc-yb)
3006            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3007            refLists = Histogram['Reflection Lists']
3008            for phase in refLists:
3009                if phase not in Phases:     #skips deleted or renamed phases silently!
3010                    continue
3011                Phase = Phases[phase]
3012                if histogram not in Phase['Histograms']:
3013                    continue
3014                im = 0
3015                if Phase['General'].get('Modulated',False):
3016                    im = 1
3017                pId = Phase['pId']
3018                phfx = '%d:%d:'%(pId,hId)
3019                refDict = refLists[phase]
3020                sumFo = 0.0
3021                sumdF = 0.0
3022                sumFosq = 0.0
3023                sumdFsq = 0.0
3024                sumInt = 0.0
3025                nExcl = 0
3026                # test to see if we are using multiprocessing below
3027                useMP,ncores = G2mp.InitMP()
3028                if len(refDict['RefList']) < 100: useMP = False       
3029                if useMP: # multiprocessing: create a set of initialized Python processes
3030                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
3031                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
3032                    profArgs = [[] for i in range(G2mp.ncores)]
3033                else:
3034                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
3035                if 'C' in calcControls[hfx+'histType']:
3036                    # are we multiprocessing?
3037                    for iref,refl in enumerate(refDict['RefList']):
3038                        if useMP: 
3039                            profArgs[iref%G2mp.ncores].append((refl,iref))
3040                        else:
3041                            icod= G2mp.ComputeFobsSqCW(refl,iref)
3042                            if type(icod) is tuple:
3043                                refl[8+im] = icod[0]
3044                                sumInt += icod[1]
3045                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3046                            elif icod == -1:
3047                                refl[3+im] *= -1
3048                                nExcl += 1
3049                            elif icod == -2:
3050                                break
3051                    if useMP:
3052                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
3053                            sumInt += sInt
3054                            for refl8im,irefl in resList:
3055                                if refl8im is None:
3056                                    refDict['RefList'][irefl][3+im] *= -1
3057                                    nExcl += 1
3058                                else:
3059                                    refDict['RefList'][irefl][8+im] = refl8im
3060                                    if parmDict[phfx+'LeBail']:
3061                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3062                elif 'T' in calcControls[hfx+'histType']:
3063                    for iref,refl in enumerate(refDict['RefList']):
3064                        if useMP: 
3065                            profArgs[iref%G2mp.ncores].append((refl,iref))
3066                        else:
3067                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
3068                            if type(icod) is tuple:
3069                                refl[8+im] = icod[0]
3070                                sumInt += icod[1]
3071                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
3072                            elif icod == -1:
3073                                refl[3+im] *= -1
3074                                nExcl += 1
3075                            elif icod == -2:
3076                                break
3077                    if useMP:
3078                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
3079                            sumInt += sInt
3080                            for refl8im,irefl in resList:
3081                                if refl8im is None:
3082                                    refDict['RefList'][irefl][3+im] *= -1
3083                                    nExcl += 1
3084                                else:
3085                                    refDict['RefList'][irefl][8+im] = refl8im
3086                                    if parmDict[phfx+'LeBail']:
3087                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
3088                if useMP: MPpool.terminate()
3089                sumFo = 0.0
3090                sumdF = 0.0
3091                sumFosq = 0.0
3092                sumdFsq = 0.0
3093                for iref,refl in enumerate(refDict['RefList']):
3094                    Fo = np.sqrt(np.abs(refl[8+im]))
3095                    Fc = np.sqrt(np.abs(refl[9]+im))
3096                    sumFo += Fo
3097                    sumFosq += refl[8+im]**2
3098                    sumdF += np.abs(Fo-Fc)
3099                    sumdFsq += (refl[8+im]-refl[9+im])**2
3100                if sumFo:
3101                    Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3102                    Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3103                else:
3104                    Histogram['Residuals'][phfx+'Rf'] = 100.
3105                    Histogram['Residuals'][phfx+'Rf^2'] = 100.
3106                Histogram['Residuals'][phfx+'sumInt'] = sumInt
3107                Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList'])-nExcl
3108                Histogram['Residuals']['hId'] = hId
3109        elif 'HKLF' in histogram[:4]:
3110            Histogram = Histograms[histogram]
3111            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
3112    if GSASIIpath.GetConfigValue('Show_timing',False):
3113        print ('GetFobsSq t=',time.time()-starttime)
3114               
3115def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup,histogram=None):
3116    'Computes the powder pattern for a histogram based on contributions from all used phases'
3117    if GSASIIpath.GetConfigValue('Show_timing',False): starttime = time.time()
3118   
3119    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
3120        U = parmDict[hfx+'U']
3121        V = parmDict[hfx+'V']
3122        W = parmDict[hfx+'W']
3123        X = parmDict[hfx+'X']
3124        Y = parmDict[hfx+'Y']
3125        Z = parmDict[hfx+'Z']
3126        tanPos = tand(refl[5+im]/2.0)
3127        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3128        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3129        sig = max(0.001,sig)
3130        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam+Z     #save peak gamma
3131        gam = max(0.001,gam)
3132        return sig,gam
3133               
3134    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
3135        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
3136            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']*refl[4+im]
3137        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2+parmDict[hfx+'Z']
3138        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
3139        sig += Ssig
3140        gam += Sgam
3141        return sig,gam
3142       
3143    def GetReflAlpBet(refl,im,hfx,parmDict):
3144        alp = parmDict[hfx+'alpha']/refl[4+im]
3145        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
3146        return alp,bet
3147       
3148    hId = Histogram['hId']
3149    hfx = ':%d:'%(hId)
3150    bakType = calcControls[hfx+'bakType']
3151    fixedBkg = {i:Histogram['Background'][1][i] for i in Histogram['Background'][1] if i.startswith("_")} 
3152    yb,Histogram['sumBk'] = G2pwd.getBackground(hfx,parmDict,bakType,calcControls[hfx+'histType'],x,fixedBkg)
3153    yc = np.zeros_like(yb)
3154    cw = np.diff(ma.getdata(x))
3155    cw = np.append(cw,cw[-1])
3156       
3157    if 'C' in calcControls[hfx+'histType']:   
3158        shl = max(parmDict[hfx+'SH/L'],0.002)
3159        Ka2 = False
3160        if hfx+'Lam1' in (parmDict.keys()):
3161            wave = parmDict[hfx+'Lam1']
3162            Ka2 = True
3163            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3164            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3165        else:
3166            wave = parmDict[hfx+'Lam']
3167    else:
3168        shl = 0.
3169    for phase in Histogram['Reflection Lists']:
3170        refDict = Histogram['Reflection Lists'][phase]
3171        if phase not in Phases:     #skips deleted or renamed phases silently!
3172            continue
3173        Phase = Phases[phase]
3174        if histogram and not histogram in Phase['Histograms']:
3175            continue
3176        pId = Phase['pId']
3177        pfx = '%d::'%(pId)
3178        phfx = '%d:%d:'%(pId,hId)
3179        hfx = ':%d:'%(hId)
3180        SGData = Phase['General']['SGData']
3181        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3182        im = 0
3183        if Phase['General'].get('Modulated',False):
3184            SSGData = Phase['General']['SSGData']
3185            im = 1  #offset in SS reflection list
3186            #??
3187        Dij = GetDij(phfx,SGData,parmDict)
3188        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]  #TODO: need to do someting if Dij << 0.
3189        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3190        if np.any(np.diag(G)<0.) or np.any(np.isnan(A)):
3191            raise G2obj.G2Exception('invalid metric tensor \n cell/Dij refinement not advised')
3192        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3193        Vst = np.sqrt(nl.det(G))    #V*
3194        if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
3195            if im:
3196                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
3197            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
3198                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3199            else:
3200                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
3201        badPeak = False
3202        # test to see if we are using multiprocessing here
3203        useMP,ncores = G2mp.InitMP()
3204        if len(refDict['RefList']) < 100: useMP = False       
3205        if useMP: # multiprocessing: create a set of initialized Python processes
3206            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
3207            profArgs = [[] for i in range(ncores)]
3208        if 'C' in calcControls[hfx+'histType']:
3209            for iref,refl in enumerate(refDict['RefList']):
3210                if im:
3211                    h,k,l,m = refl[:4]
3212                else:
3213                    h,k,l = refl[:3]
3214                Uniq = np.inner(refl[:3],SGMT)
3215                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
3216                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
3217                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3218                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3219                refl[11+im] *= Vst*Lorenz
3220                 
3221                if Phase['General'].get('doPawley'):
3222                    try:
3223                        if im:
3224                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3225                        else:
3226                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3227                        refl[9+im] = parmDict[pInd]
3228                    except KeyError:
3229#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3230                        continue
3231                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
3232                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3233                iFin = np.searchsorted(x,refl[5+im]+fmax)
3234                if not iBeg+iFin:       #peak below low limit - skip peak
3235                    continue
3236                elif not iBeg-iFin:     #peak above high limit - done
3237                    break
3238                elif iBeg > iFin:   #bad peak coeff - skip
3239                    badPeak = True
3240                    continue
3241                if useMP:
3242                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
3243                else:
3244                    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
3245                if Ka2:
3246                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3247                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
3248                    iBeg = np.searchsorted(x,pos2-fmin)
3249                    iFin = np.searchsorted(x,pos2+fmax)
3250                    if not iBeg+iFin:       #peak below low limit - skip peak
3251                        continue
3252                    elif not iBeg-iFin:     #peak above high limit - done
3253                        return yc,yb
3254                    elif iBeg > iFin:   #bad peak coeff - skip
3255                        continue
3256                    if useMP:
3257                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
3258                    else:
3259                        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
3260        elif 'T' in calcControls[hfx+'histType']:
3261            for iref,refl in enumerate(refDict['RefList']):
3262                if im:
3263                    h,k,l,m = refl[:4]
3264                else:
3265                    h,k,l = refl[:3]
3266                Uniq = np.inner(refl[:3],SGMT)
3267                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
3268                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
3269#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
3270                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
3271                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
3272                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3273                refl[11+im] *= Vst*Lorenz
3274                if Phase['General'].get('doPawley'):
3275                    try:
3276                        if im:
3277                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
3278                        else:
3279                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3280                        refl[9+im] = parmDict[pInd]
3281                    except KeyError:
3282#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3283                        continue
3284                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
3285                iBeg = np.searchsorted(x,refl[5+im]-fmin)
3286                iFin = np.searchsorted(x,refl[5+im]+fmax)
3287                if not iBeg+iFin:       #peak below low limit - skip peak
3288                    continue
3289                elif not iBeg-iFin:     #peak above high limit - done
3290                    break
3291                elif iBeg > iFin:   #bad peak coeff - skip
3292                    badPeak = True
3293                    continue
3294                if useMP:
3295                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin))
3296                else:
3297                    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]
3298#        print 'profile calc time: %.3fs'%(time.time()-time0)
3299        if useMP and 'C' in calcControls[hfx+'histType']:
3300            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
3301                yc += y
3302            MPpool.terminate()
3303        elif useMP:
3304            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
3305                yc += y
3306            MPpool.terminate()
3307    if badPeak:
3308        print ('ouch #4 bad profile coefficients yield negative peak width; some reflections skipped')
3309    if GSASIIpath.GetConfigValue('Show_timing',False):
3310        print ('getPowderProfile t=%.3f'%(time.time()-starttime))
3311    return yc,yb
3312   
3313def getPowderProfileDervMP(args):
3314    '''Computes the derivatives of the computed powder pattern with respect to all
3315    refined parameters.
3316    Multiprocessing version.
3317    '''
3318    import pytexture as ptx
3319    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics for each processor
3320    parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars = args[:9]
3321    prc,tprc,histogram = 0,1,None
3322    if len(args) >= 10: prc=args[9]
3323    if len(args) >= 11: tprc=args[10]
3324    if len(args) >= 12: histogram=args[11]
3325    def cellVaryDerv(pfx,SGData,dpdA): 
3326        if SGData['SGLaue'] in ['-1',]:
3327            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
3328                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
3329        elif SGData['SGLaue'] in ['2/m',]:
3330            if SGData['SGUniq'] == 'a':
3331                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA