source: trunk/GSASIIstrMath.py @ 3848

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

since id is a python routine which could be useful in debugging, replaced id with Id (or pid in one place)

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