source: trunk/GSASIIstrMath.py

Last change on this file was 5586, checked in by vondreele, 7 days ago

changes to allow generalized atomic form factors (beginning)
Some fixes to spin RB stuff

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