source: branch/2frame/GSASIImath.py @ 2967

Last change on this file since 2967 was 2967, checked in by toby, 6 years ago

misc changes for sphinx docs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 181.4 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2017-08-07 19:17:51 +0000 (Mon, 07 Aug 2017) $
5# $Author: toby $
6# $Revision: 2967 $
7# $URL: branch/2frame/GSASIImath.py $
8# $Id: GSASIImath.py 2967 2017-08-07 19:17:51Z toby $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import random as rn
18import numpy as np
19import numpy.linalg as nl
20import numpy.ma as ma
21import time
22import math
23import copy
24import GSASIIpath
25GSASIIpath.SetVersionNumber("$Revision: 2967 $")
26import GSASIIElem as G2el
27import GSASIIlattice as G2lat
28import GSASIIspc as G2spc
29import GSASIIpwd as G2pwd
30import numpy.fft as fft
31import scipy.optimize as so
32import pypowder as pwd
33
34sind = lambda x: np.sin(x*np.pi/180.)
35cosd = lambda x: np.cos(x*np.pi/180.)
36tand = lambda x: np.tan(x*np.pi/180.)
37asind = lambda x: 180.*np.arcsin(x)/np.pi
38acosd = lambda x: 180.*np.arccos(x)/np.pi
39atand = lambda x: 180.*np.arctan(x)/np.pi
40atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
41twopi = 2.0*np.pi
42twopisq = 2.0*np.pi**2
43nxs = np.newaxis
44   
45################################################################################
46##### Hessian least-squares Levenberg-Marquardt routine
47################################################################################
48
49def pinv(a, rcond=1e-15 ):
50    """
51    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
52    Modified from numpy.linalg.pinv; assumes a is Hessian & returns no. zeros found
53    Calculate the generalized inverse of a matrix using its
54    singular-value decomposition (SVD) and including all
55    *large* singular values.
56
57    :param array a: (M, M) array_like - here assumed to be LS Hessian
58      Matrix to be pseudo-inverted.
59    :param float rcond: Cutoff for small singular values.
60      Singular values smaller (in modulus) than
61      `rcond` * largest_singular_value (again, in modulus)
62      are set to zero.
63
64    :returns: B : (M, M) ndarray
65      The pseudo-inverse of `a`
66
67    Raises: LinAlgError
68      If the SVD computation does not converge.
69
70    Notes:
71      The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
72      defined as: "the matrix that 'solves' [the least-squares problem]
73      :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
74      :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
75
76    It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
77    value decomposition of A, then
78    :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
79    orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
80    of A's so-called singular values, (followed, typically, by
81    zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
82    consisting of the reciprocals of A's singular values
83    (again, followed by zeros). [1]
84
85    References:
86    .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
87           FL, Academic Press, Inc., 1980, pp. 139-142.
88
89    """
90    u, s, vt = nl.svd(a, 0)
91    cutoff = rcond*np.maximum.reduce(s)
92    s = np.where(s>cutoff,1./s,0.)
93    nzero = s.shape[0]-np.count_nonzero(s)
94#    res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))
95    res = np.dot(vt.T,s[:,nxs]*u.T)
96    return res,nzero
97
98def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
99   
100    """
101    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
102    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`   
103    where :math:`x = arg min(\sum_{y=0}^{N_{obs}} (func(y)^2,axis=0))`
104
105    :param function func: callable method or function
106        should take at least one (possibly length N vector) argument and
107        returns M floating point numbers.
108    :param np.ndarray x0: The starting estimate for the minimization of length N
109    :param function Hess: callable method or function
110        A required function or method to compute the weighted vector and Hessian for func.
111        It must be a symmetric NxN array
112    :param tuple args: Any extra arguments to func are placed in this tuple.
113    :param float ftol: Relative error desired in the sum of squares.
114    :param float xtol: Relative tolerance of zeros in the SVD solution in nl.pinv.
115    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
116        until other limits are met (ftol, xtol)
117    :param int lamda: initial Marquardt lambda=10**lamda
118    :param bool Print: True for printing results (residuals & times) by cycle
119
120    :returns: (x,cov_x,infodict) where
121
122      * x : ndarray
123        The solution (or the result of the last iteration for an unsuccessful
124        call).
125      * cov_x : ndarray
126        Uses the fjac and ipvt optional outputs to construct an
127        estimate of the jacobian around the solution.  ``None`` if a
128        singular matrix encountered (indicates very flat curvature in
129        some direction).  This matrix must be multiplied by the
130        residual standard deviation to get the covariance of the
131        parameter estimates -- see curve_fit.
132      * infodict : dict
133        a dictionary of optional outputs with the keys:
134
135         * 'fvec' : the function evaluated at the output
136         * 'num cyc':
137         * 'nfev':
138         * 'lamMax':
139         * 'psing':
140         * 'SVD0':
141           
142    """
143               
144    ifConverged = False
145    deltaChi2 = -10.
146    x0 = np.array(x0, ndmin=1)      #might be redundant?
147    n = len(x0)
148    if type(args) != type(()):
149        args = (args,)
150       
151    icycle = 0
152    One = np.ones((n,n))
153    lam = 10.**lamda
154    lamMax = lam
155    nfev = 0
156    if Print:
157        print ' Hessian Levenburg-Marquardt SVD refinement on %d variables:'%(n)
158    Lam = np.zeros((n,n))
159    while icycle < maxcyc:
160        time0 = time.time()
161        M = func(x0,*args)
162        nfev += 1
163        chisq0 = np.sum(M**2)
164        Yvec,Amat = Hess(x0,*args)
165        Adiag = np.sqrt(np.diag(Amat))
166        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
167        if np.any(psing):                #hard singularity in matrix
168            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
169        Anorm = np.outer(Adiag,Adiag)
170        Yvec /= Adiag
171        Amat /= Anorm
172        if Print:
173            print 'initial chi^2 %.5g'%(chisq0)
174        chitol = ftol
175        while True:
176            Lam = np.eye(Amat.shape[0])*lam
177            Amatlam = Amat*(One+Lam)
178            try:
179                Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
180            except nl.LinAlgError:
181                print 'ouch #1 bad SVD inversion; change parameterization'
182                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
183                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
184            Xvec = np.inner(Ainv,Yvec)      #solve
185            Xvec /= Adiag
186            M2 = func(x0+Xvec,*args)
187            nfev += 1
188            chisq1 = np.sum(M2**2)
189            if chisq1 > chisq0*(1.+chitol):
190                lam *= 10.
191                if Print:
192                    print 'new chi^2 %.5g, %d SVD zeros ; matrix modification needed; lambda now %.1e'%(chisq1,Nzeros,lam)
193            else:
194                x0 += Xvec
195                lam /= 10.
196                break
197            if lam > 10.e3:
198                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
199                break
200            chitol *= 2
201        lamMax = max(lamMax,lam)
202        deltaChi2 = (chisq0-chisq1)/chisq0
203        if Print:
204            print ' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Lambda: %.3g,  Delta: %.3g'%(
205                icycle,time.time()-time0,chisq1,lamMax,deltaChi2)
206        if deltaChi2 < ftol:
207            ifConverged = True
208            if Print: print "converged"
209            break
210        icycle += 1
211    M = func(x0,*args)
212    nfev += 1
213    Yvec,Amat = Hess(x0,*args)
214    Adiag = np.sqrt(np.diag(Amat))
215    Anorm = np.outer(Adiag,Adiag)
216    Lam = np.eye(Amat.shape[0])*lam
217    Amatlam = Amat/Anorm       
218    try:
219        Bmat,Nzero = pinv(Amatlam,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
220        if Print:
221            print 'Found %d SVD zeros'%(Nzero)
222#        Bmat = nl.inv(Amatlam); Nzeros = 0
223        Bmat = Bmat/Anorm
224        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
225    except nl.LinAlgError:
226        print 'ouch #2 linear algebra error in making v-cov matrix'
227        psing = []
228        if maxcyc:
229            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
230        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]         
231           
232def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
233   
234    """
235    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
236    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`   
237    where :math:`x = arg min(\sum_{y=0}^{N_{obs}} (func(y)^2,axis=0))`
238
239    :param function func: callable method or function
240        should take at least one (possibly length N vector) argument and
241        returns M floating point numbers.
242    :param np.ndarray x0: The starting estimate for the minimization of length N
243    :param function Hess: callable method or function
244        A required function or method to compute the weighted vector and Hessian for func.
245        It must be a symmetric NxN array
246    :param tuple args: Any extra arguments to func are placed in this tuple.
247    :param float ftol: Relative error desired in the sum of squares.
248    :param float xtol: Relative tolerance of zeros in the SVD solution in nl.pinv.
249    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
250        until other limits are met (ftol, xtol)
251    :param bool Print: True for printing results (residuals & times) by cycle
252
253    :returns: (x,cov_x,infodict) where
254
255      * x : ndarray
256        The solution (or the result of the last iteration for an unsuccessful
257        call).
258      * cov_x : ndarray
259        Uses the fjac and ipvt optional outputs to construct an
260        estimate of the jacobian around the solution.  ``None`` if a
261        singular matrix encountered (indicates very flat curvature in
262        some direction).  This matrix must be multiplied by the
263        residual standard deviation to get the covariance of the
264        parameter estimates -- see curve_fit.
265      * infodict : dict
266        a dictionary of optional outputs with the keys:
267
268         * 'fvec' : the function evaluated at the output
269         * 'num cyc':
270         * 'nfev':
271         * 'lamMax':0.
272         * 'psing':
273         * 'SVD0':
274           
275    """
276               
277    ifConverged = False
278    deltaChi2 = -10.
279    x0 = np.array(x0, ndmin=1)      #might be redundant?
280    n = len(x0)
281    if type(args) != type(()):
282        args = (args,)
283       
284    icycle = 0
285    nfev = 0
286    if Print:
287        print ' Hessian SVD refinement on %d variables:'%(n)
288    while icycle < maxcyc:
289        time0 = time.time()
290        M = func(x0,*args)
291        nfev += 1
292        chisq0 = np.sum(M**2)
293        Yvec,Amat = Hess(x0,*args)
294        Adiag = np.sqrt(np.diag(Amat))
295        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
296        if np.any(psing):                #hard singularity in matrix
297            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
298        Anorm = np.outer(Adiag,Adiag)
299        Yvec /= Adiag
300        Amat /= Anorm
301        if Print:
302            print 'initial chi^2 %.5g'%(chisq0)
303        try:
304            Ainv,Nzeros = pinv(Amat,xtol)    #do Moore-Penrose inversion (via SVD)
305        except nl.LinAlgError:
306            print 'ouch #1 bad SVD inversion; change parameterization'
307            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
308            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
309        Xvec = np.inner(Ainv,Yvec)      #solve
310        Xvec /= Adiag
311        M2 = func(x0+Xvec,*args)
312        nfev += 1
313        chisq1 = np.sum(M2**2)
314        deltaChi2 = (chisq0-chisq1)/chisq0
315        if Print:
316            print ' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Delta: %.3g'%(
317                icycle,time.time()-time0,chisq1,deltaChi2)
318        if deltaChi2 < ftol:
319            ifConverged = True
320            if Print: print "converged"
321            break
322        icycle += 1
323    M = func(x0,*args)
324    nfev += 1
325    Yvec,Amat = Hess(x0,*args)
326    Adiag = np.sqrt(np.diag(Amat))
327    Anorm = np.outer(Adiag,Adiag)
328    Amat = Amat/Anorm       
329    try:
330        Bmat,Nzero = pinv(Amat,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
331        print 'Found %d SVD zeros'%(Nzero)
332#        Bmat = nl.inv(Amatlam); Nzeros = 0
333        Bmat = Bmat/Anorm
334        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':[],
335            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
336    except nl.LinAlgError:
337        print 'ouch #2 linear algebra error in making v-cov matrix'
338        psing = []
339        if maxcyc:
340            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
341        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]         
342           
343def getVCov(varyNames,varyList,covMatrix):
344    '''obtain variance-covariance terms for a set of variables. NB: the varyList
345    and covMatrix were saved by the last least squares refinement so they must match.
346   
347    :param list varyNames: variable names to find v-cov matric for
348    :param list varyList: full list of all variables in v-cov matrix
349    :param nparray covMatrix: full variance-covariance matrix from the last
350     least squares refinement
351   
352    :returns: nparray vcov: variance-covariance matrix for the variables given
353     in varyNames
354   
355    '''
356    vcov = np.zeros((len(varyNames),len(varyNames)))
357    for i1,name1 in enumerate(varyNames):
358        for i2,name2 in enumerate(varyNames):
359            try:
360                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
361            except ValueError:
362                vcov[i1][i2] = 0.0
363#                if i1 == i2:
364#                    vcov[i1][i2] = 1e-20
365#                else:
366#                    vcov[i1][i2] = 0.0
367    return vcov
368   
369################################################################################
370##### Atom manipulations
371################################################################################
372
373def FindMolecule(ind,generalData,atomData):                    #uses numpy & masks - very fast even for proteins!
374
375    def getNeighbors(atom,radius):
376        Dx = UAtoms-np.array(atom[cx:cx+3])
377        dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of disorder "bonds" < 0.5A
378        sumR = Radii+radius
379        return set(ma.nonzero(ma.masked_greater(dist-factor*sumR,0.))[0])                #get indices of bonded atoms
380
381    import numpy.ma as ma
382    indices = (-1,0,1)
383    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices],dtype='f')
384    cx,ct,cs,ci = generalData['AtomPtrs']
385    DisAglCtls = generalData['DisAglCtls']
386    SGData = generalData['SGData']
387    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
388    radii = DisAglCtls['BondRadii']
389    atomTypes = DisAglCtls['AtomTypes']
390    factor = DisAglCtls['Factors'][0]
391    unit = np.zeros(3)
392    try:
393        indH = atomTypes.index('H')
394        radii[indH] = 0.5
395    except:
396        pass           
397    nAtom = len(atomData)
398    Indx = range(nAtom)
399    UAtoms = []
400    Radii = []
401    for atom in atomData:
402        UAtoms.append(np.array(atom[cx:cx+3]))
403        Radii.append(radii[atomTypes.index(atom[ct])])
404    UAtoms = np.array(UAtoms)
405    Radii = np.array(Radii)
406    for nOp,Op in enumerate(SGData['SGOps'][1:]):
407        UAtoms = np.concatenate((UAtoms,(np.inner(Op[0],UAtoms[:nAtom]).T+Op[1])))
408        Radii = np.concatenate((Radii,Radii[:nAtom]))
409        Indx += Indx[:nAtom]
410    for icen,cen in enumerate(SGData['SGCen'][1:]):
411        UAtoms = np.concatenate((UAtoms,(UAtoms+cen)))
412        Radii = np.concatenate((Radii,Radii))
413        Indx += Indx[:nAtom]
414    if SGData['SGInv']:
415        UAtoms = np.concatenate((UAtoms,-UAtoms))
416        Radii = np.concatenate((Radii,Radii))
417        Indx += Indx
418    UAtoms %= 1.
419    mAtoms = len(UAtoms)
420    for unit in Units:
421        if np.any(unit):    #skip origin cell
422            UAtoms = np.concatenate((UAtoms,UAtoms[:mAtoms]+unit))
423            Radii = np.concatenate((Radii,Radii[:mAtoms]))
424            Indx += Indx[:mAtoms]
425    UAtoms = np.array(UAtoms)
426    Radii = np.array(Radii)
427    newAtoms = [atomData[ind],]
428    atomData[ind] = None
429    radius = Radii[ind]
430    IndB = getNeighbors(newAtoms[-1],radius)
431    while True:
432        if not len(IndB):
433            break
434        indb = IndB.pop()
435        if atomData[Indx[indb]] == None:
436            continue
437        while True:
438            try:
439                jndb = IndB.index(indb)
440                IndB.remove(jndb)
441            except:
442                break
443        newAtom = copy.copy(atomData[Indx[indb]])
444        newAtom[cx:cx+3] = UAtoms[indb]     #NB: thermal Uij, etc. not transformed!
445        newAtoms.append(newAtom)
446        atomData[Indx[indb]] = None
447        IndB = set(list(IndB)+list(getNeighbors(newAtoms[-1],radius)))
448        if len(IndB) > nAtom:
449            return 'Assemble molecule cannot be used on extended structures'
450    for atom in atomData:
451        if atom != None:
452            newAtoms.append(atom)
453    return newAtoms
454       
455def FindAtomIndexByIDs(atomData,loc,IDs,Draw=True):
456    '''finds the set of atom array indices for a list of atom IDs. Will search
457    either the Atom table or the drawAtom table.
458   
459    :param list atomData: Atom or drawAtom table containting coordinates, etc.
460    :param int loc: location of atom id in atomData record
461    :param list IDs: atom IDs to be found
462    :param bool Draw: True if drawAtom table to be searched; False if Atom table
463      is searched
464   
465    :returns: list indx: atom (or drawAtom) indices
466   
467    '''
468    indx = []
469    for i,atom in enumerate(atomData):
470        if Draw and atom[loc] in IDs:
471            indx.append(i)
472        elif atom[loc] in IDs:
473            indx.append(i)
474    return indx
475
476def FillAtomLookUp(atomData,indx):
477    '''create a dictionary of atom indexes with atom IDs as keys
478   
479    :param list atomData: Atom table to be used
480   
481    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
482   
483    '''
484    atomLookUp = {}
485    for iatm,atom in enumerate(atomData):
486        atomLookUp[atom[indx]] = iatm
487    return atomLookUp
488
489def GetAtomsById(atomData,atomLookUp,IdList):
490    '''gets a list of atoms from Atom table that match a set of atom IDs
491   
492    :param list atomData: Atom table to be used
493    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
494    :param list IdList: atom IDs to be found
495   
496    :returns: list atoms: list of atoms found
497   
498    '''
499    atoms = []
500    for id in IdList:
501        atoms.append(atomData[atomLookUp[id]])
502    return atoms
503   
504def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
505    '''gets atom parameters for atoms using atom IDs
506   
507    :param list atomData: Atom table to be used
508    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
509    :param list IdList: atom IDs to be found
510    :param int itemLoc: pointer to desired 1st item in an atom table entry
511    :param int numItems: number of items to be retrieved
512   
513    :returns: type name: description
514   
515    '''
516    Items = []
517    if not isinstance(IdList,list):
518        IdList = [IdList,]
519    for id in IdList:
520        if numItems == 1:
521            Items.append(atomData[atomLookUp[id]][itemLoc])
522        else:
523            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
524    return Items
525   
526def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
527    '''default doc string
528   
529    :param type name: description
530   
531    :returns: type name: description
532   
533    '''
534    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
535    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
536    XYZ = []
537    for ind in indx:
538        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
539        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
540        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
541    return XYZ
542   
543#def TransformAtoms(Atoms,cx,cia,Trans,Vec):
544#    for Atom in Atoms:
545#        XYZ = Atom[cx:cx+3]
546#        if 'A' in Atom[cia]:
547#            U6 = Atom[cia+2:cia+8]
548   
549
550def FindNeighbors(phase,FrstName,AtNames,notName=''):
551    General = phase['General']
552    cx,ct,cs,cia = General['AtomPtrs']
553    Atoms = phase['Atoms']
554    atNames = [atom[ct-1] for atom in Atoms]
555    Cell = General['Cell'][1:7]
556    Amat,Bmat = G2lat.cell2AB(Cell)
557    atTypes = General['AtomTypes']
558    Radii = np.array(General['BondRadii'])
559    DisAglCtls = General['DisAglCtls']   
560    radiusFactor = DisAglCtls['Factors'][0]
561    AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
562    Orig = atNames.index(FrstName)
563    OId = Atoms[Orig][cia+8]
564    OType = Atoms[Orig][ct]
565    XYZ = getAtomXYZ(Atoms,cx)       
566    Neigh = []
567    Ids = []
568    Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
569    dist = np.sqrt(np.sum(Dx**2,axis=1))
570    sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
571    IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))
572    for j in IndB[0]:
573        if j != Orig:
574            if AtNames[j] != notName:
575                Neigh.append([AtNames[j],dist[j],True])
576                Ids.append(Atoms[j][cia+8])
577    return Neigh,[OId,Ids]
578   
579def FindAllNeighbors(phase,FrstName,AtNames,notName=''):
580    General = phase['General']
581    cx,ct,cs,cia = General['AtomPtrs']
582    Atoms = phase['Atoms']
583    atNames = [atom[ct-1] for atom in Atoms]
584    Cell = General['Cell'][1:7]
585    Amat,Bmat = G2lat.cell2AB(Cell)
586    SGData = General['SGData']
587    indices = (-1,0,1)
588    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
589    atTypes = General['AtomTypes']
590    Radii = np.array(General['BondRadii'])
591    DisAglCtls = General['DisAglCtls']   
592    radiusFactor = DisAglCtls['Factors'][0]
593    AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
594    Orig = atNames.index(FrstName)
595    OId = Atoms[Orig][cia+8]
596    OType = Atoms[Orig][ct]
597    XYZ = getAtomXYZ(Atoms,cx)       
598    Oxyz = XYZ[Orig]
599    Neigh = []
600    Ids = []
601    sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
602    sumR = np.reshape(np.tile(sumR,27),(27,-1))
603    results = []
604    for xyz in XYZ:
605        results.append(G2spc.GenAtom(xyz,SGData,False,Move=False))
606    for iA,result in enumerate(results):
607        if iA != Orig:               
608            for [Txyz,Top,Tunit] in result:
609                Dx = np.array([Txyz-Oxyz+unit for unit in Units])
610                dx = np.inner(Dx,Amat)
611                dist = np.sqrt(np.sum(dx**2,axis=1))
612                IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR[:,iA],0.))
613        #        GSASIIpath.IPyBreak()
614                for iU in IndB[0]:
615                    if AtNames[iA] != notName:
616                        unit = Units[iU]
617                        if np.any(unit):
618                            Topstr = ' +(%4d)[%2d,%2d,%2d]'%(Top,unit[0],unit[1],unit[2])
619                        else:
620                            Topstr = ' +(%4d)'%(Top)
621                        Neigh.append([AtNames[iA]+Topstr,dist[iU]])
622                        Ids.append(Atoms[iA][cia+8])
623    return Neigh,[OId,Ids]
624   
625def calcBond(A,Ax,Bx,MTCU):
626    cell = G2lat.A2cell(A)
627    Amat,Bmat = G2lat.cell2AB(cell)
628    M,T,C,U = MTCU
629    Btx = np.inner(M,Bx)+T+C+U
630    Dx = Btx-Ax
631    dist = np.sqrt(np.inner(Amat,Dx))
632    return dist
633   
634def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
635   
636    def getTransMat(RXYZ,OXYZ,TXYZ,Amat):
637        Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T           
638        Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,nxs]
639        Mat2 = np.cross(Vec[0],Vec[1])      #UxV
640        Mat2 /= np.sqrt(np.sum(Mat2**2))
641        Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
642        return nl.inv(np.array([Vec[0],Mat2,Mat3]))       
643   
644    cx,ct,cs,cia = General['AtomPtrs']
645    Cell = General['Cell'][1:7]
646    Amat,Bmat = G2lat.cell2AB(Cell)
647    nBonds = AddHydId[-1]+len(AddHydId[1])
648    Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0]
649    OXYZ = np.array(Oatom[cx:cx+3])
650    if 'I' in Oatom[cia]:
651        Uiso = Oatom[cia+1]
652    else:
653        Uiso = (Oatom[cia+2]+Oatom[cia+3]+Oatom[cia+4])/3.0       #simple average
654    Uiso = max(Uiso,0.005)                      #set floor!
655    Tatoms = GetAtomsById(Atoms,AtLookUp,AddHydId[1])
656    TXYZ = np.array([tatom[cx:cx+3] for tatom in Tatoms]) #3 x xyz
657    if nBonds == 4:
658        if AddHydId[-1] == 1:
659            Vec = TXYZ-OXYZ
660            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0))
661            Vec = np.sum(Vec/Len,axis=0)
662            Len = np.sqrt(np.sum(Vec**2))
663            Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len
664            HU = 1.1*Uiso
665            return [Hpos,],[HU,]
666        elif AddHydId[-1] == 2:
667            Vec = np.inner(Amat,TXYZ-OXYZ).T
668            Vec[0] += Vec[1]            #U - along bisector
669            Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,nxs]
670            Mat2 = np.cross(Vec[0],Vec[1])      #UxV
671            Mat2 /= np.sqrt(np.sum(Mat2**2))
672            Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
673            iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
674            Hpos = np.array([[-0.97*cosd(54.75),0.97*sind(54.75),0.],
675                [-0.97*cosd(54.75),-0.97*sind(54.75),0.]])
676            HU = 1.2*Uiso*np.ones(2)
677            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
678            return Hpos,HU
679        else:
680            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
681            RXYZ = np.array(Ratom[cx:cx+3])
682            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
683            a = 0.96*cosd(70.5)
684            b = 0.96*sind(70.5)
685            Hpos = np.array([[a,0.,-b],[a,-b*cosd(30.),0.5*b],[a,b*cosd(30.),0.5*b]])
686            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
687            HU = 1.5*Uiso*np.ones(3)
688            return Hpos,HU         
689    elif nBonds == 3:
690        if AddHydId[-1] == 1:
691            Vec = np.sum(TXYZ-OXYZ,axis=0)               
692            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
693            Vec = -0.93*Vec/Len
694            Hpos = OXYZ+Vec
695            HU = 1.1*Uiso
696            return [Hpos,],[HU,]
697        elif AddHydId[-1] == 2:
698            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
699            RXYZ = np.array(Ratom[cx:cx+3])
700            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
701            a = 0.93*cosd(60.)
702            b = 0.93*sind(60.)
703            Hpos = [[a,b,0],[a,-b,0]]
704            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
705            HU = 1.2*Uiso*np.ones(2)
706            return Hpos,HU
707    else:   #2 bonds
708        if 'C' in Oatom[ct]:
709            Vec = TXYZ[0]-OXYZ
710            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
711            Vec = -0.93*Vec/Len
712            Hpos = OXYZ+Vec
713            HU = 1.1*Uiso
714            return [Hpos,],[HU,]
715        elif 'O' in Oatom[ct]:
716            mapData = General['Map']
717            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
718            RXYZ = np.array(Ratom[cx:cx+3])
719            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
720            a = 0.82*cosd(70.5)
721            b = 0.82*sind(70.5)
722            azm = np.arange(0.,360.,5.)
723            Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm])
724            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
725            Rhos = np.array([getRho(pos,mapData) for pos in Hpos])
726            imax = np.argmax(Rhos)
727            HU = 1.5*Uiso
728            return [Hpos[imax],],[HU,]
729    return [],[]
730       
731#def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
732#    '''default doc string
733#   
734#    :param type name: description
735#   
736#    :returns: type name: description
737#   
738#    '''
739#    for atom in atomData:
740#        XYZ = np.inner(Amat,atom[cx:cx+3])
741#        if atom[cia] == 'A':
742#            UIJ = atom[cia+2:cia+8]
743               
744def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
745    '''default doc string
746   
747    :param type name: description
748   
749    :returns: type name: description
750   
751    '''
752    TLStype,TLS = rbObj['ThermalMotion'][:2]
753    Tmat = np.zeros((3,3))
754    Lmat = np.zeros((3,3))
755    Smat = np.zeros((3,3))
756    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
757        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
758    if 'T' in TLStype:
759        Tmat = G2lat.U6toUij(TLS[:6])
760    if 'L' in TLStype:
761        Lmat = G2lat.U6toUij(TLS[6:12])
762    if 'S' in TLStype:
763        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
764    XYZ = np.inner(Amat,xyz)
765    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
766    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
767    beta = np.inner(np.inner(g,Umat),g)
768    return G2lat.UijtoU6(beta)*gvec   
769       
770def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
771    '''default doc string
772   
773    :param type name: description
774   
775    :returns: type name: description
776   
777    '''
778    cx,ct,cs,cia = atPtrs
779    TLStype,TLS = rbObj['ThermalMotion'][:2]
780    Tmat = np.zeros((3,3))
781    Lmat = np.zeros((3,3))
782    Smat = np.zeros((3,3))
783    G,g = G2lat.A2Gmat(Amat)
784    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
785    if 'T' in TLStype:
786        Tmat = G2lat.U6toUij(TLS[:6])
787    if 'L' in TLStype:
788        Lmat = G2lat.U6toUij(TLS[6:12])
789    if 'S' in TLStype:
790        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
791    for atom in atomData:
792        XYZ = np.inner(Amat,atom[cx:cx+3])
793        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
794        if 'U' in TLStype:
795            atom[cia+1] = TLS[0]
796            atom[cia] = 'I'
797        else:
798            atom[cia] = 'A'
799            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
800            beta = np.inner(np.inner(g,Umat),g)
801            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
802
803def GetXYZDist(xyz,XYZ,Amat):
804    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
805        and are in crystal coordinates; Amat is crystal to Cart matrix
806   
807    :param type name: description
808   
809    :returns: type name: description
810   
811    '''
812    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
813
814def getAtomXYZ(atoms,cx):
815    '''default doc string
816   
817    :param type name: description
818   
819    :returns: type name: description
820   
821    '''
822    XYZ = []
823    for atom in atoms:
824        XYZ.append(atom[cx:cx+3])
825    return np.array(XYZ)
826
827def RotateRBXYZ(Bmat,Cart,oriQ):
828    '''rotate & transform cartesian coordinates to crystallographic ones
829    no translation applied. To be used for numerical derivatives
830   
831    :param type name: description
832   
833    :returns: type name: description
834   
835    '''
836    ''' returns crystal coordinates for atoms described by RBObj
837    '''
838    XYZ = np.zeros_like(Cart)
839    for i,xyz in enumerate(Cart):
840        XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
841    return XYZ
842
843def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
844    '''default doc string
845   
846    :param type name: description
847   
848    :returns: type name: description
849   
850    '''
851    ''' returns crystal coordinates for atoms described by RBObj
852    '''
853    RBRes = RBData[RBType][RBObj['RBId']]
854    if RBType == 'Vector':
855        vecs = RBRes['rbVect']
856        mags = RBRes['VectMag']
857        Cart = np.zeros_like(vecs[0])
858        for vec,mag in zip(vecs,mags):
859            Cart += vec*mag
860    elif RBType == 'Residue':
861        Cart = np.array(RBRes['rbXYZ'])
862        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
863            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
864            Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
865    XYZ = np.zeros_like(Cart)
866    for i,xyz in enumerate(Cart):
867        XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
868    return XYZ,Cart
869
870def UpdateMCSAxyz(Bmat,MCSA):
871    '''default doc string
872   
873    :param type name: description
874   
875    :returns: type name: description
876   
877    '''
878    xyz = []
879    atTypes = []
880    iatm = 0
881    for model in MCSA['Models'][1:]:        #skip the MD model
882        if model['Type'] == 'Atom':
883            xyz.append(model['Pos'][0])
884            atTypes.append(model['atType'])
885            iatm += 1
886        else:
887            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
888            Pos = np.array(model['Pos'][0])
889            Ori = np.array(model['Ori'][0])
890            Qori = AVdeg2Q(Ori[0],Ori[1:])
891            if model['Type'] == 'Vector':
892                vecs = RBRes['rbVect']
893                mags = RBRes['VectMag']
894                Cart = np.zeros_like(vecs[0])
895                for vec,mag in zip(vecs,mags):
896                    Cart += vec*mag
897            elif model['Type'] == 'Residue':
898                Cart = np.array(RBRes['rbXYZ'])
899                for itor,seq in enumerate(RBRes['rbSeq']):
900                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
901                    Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
902            if model['MolCent'][1]:
903                Cart -= model['MolCent'][0]
904            for i,x in enumerate(Cart):
905                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
906                atType = RBRes['rbTypes'][i]
907                atTypes.append(atType)
908                iatm += 1
909    return np.array(xyz),atTypes
910   
911def SetMolCent(model,RBData):
912    '''default doc string
913   
914    :param type name: description
915   
916    :returns: type name: description
917   
918    '''
919    rideList = []
920    RBRes = RBData[model['Type']][model['RBId']]
921    if model['Type'] == 'Vector':
922        vecs = RBRes['rbVect']
923        mags = RBRes['VectMag']
924        Cart = np.zeros_like(vecs[0])
925        for vec,mag in zip(vecs,mags):
926            Cart += vec*mag
927    elif model['Type'] == 'Residue':
928        Cart = np.array(RBRes['rbXYZ'])
929        for seq in RBRes['rbSeq']:
930            rideList += seq[3]
931    centList = set(range(len(Cart)))-set(rideList)
932    cent = np.zeros(3)
933    for i in centList:
934        cent += Cart[i]
935    model['MolCent'][0] = cent/len(centList)   
936   
937def UpdateRBUIJ(Bmat,Cart,RBObj):
938    '''default doc string
939   
940    :param type name: description
941   
942    :returns: type name: description
943   
944    '''
945    ''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
946    '''
947    TLStype,TLS = RBObj['ThermalMotion'][:2]
948    T = np.zeros(6)
949    L = np.zeros(6)
950    S = np.zeros(8)
951    if 'T' in TLStype:
952        T = TLS[:6]
953    if 'L' in TLStype:
954        L = np.array(TLS[6:12])*(np.pi/180.)**2
955    if 'S' in TLStype:
956        S = np.array(TLS[12:])*(np.pi/180.)
957    g = nl.inv(np.inner(Bmat,Bmat))
958    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
959        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
960    Uout = []
961    Q = RBObj['Orient'][0]
962    for X in Cart:
963        X = prodQVQ(Q,X)
964        if 'U' in TLStype:
965            Uout.append(['I',TLS[0],0,0,0,0,0,0])
966        elif not 'N' in TLStype:
967            U = [0,0,0,0,0,0]
968            U[0] = T[0]+L[1]*X[2]**2+L[2]*X[1]**2-2.0*L[5]*X[1]*X[2]+2.0*(S[2]*X[2]-S[4]*X[1])
969            U[1] = T[1]+L[0]*X[2]**2+L[2]*X[0]**2-2.0*L[4]*X[0]*X[2]+2.0*(S[5]*X[0]-S[0]*X[2])
970            U[2] = T[2]+L[1]*X[0]**2+L[0]*X[1]**2-2.0*L[3]*X[1]*X[0]+2.0*(S[1]*X[1]-S[3]*X[0])
971            U[3] = T[3]+L[4]*X[1]*X[2]+L[5]*X[0]*X[2]-L[3]*X[2]**2-L[2]*X[0]*X[1]+  \
972                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
973            U[4] = T[4]+L[3]*X[1]*X[2]+L[5]*X[0]*X[1]-L[4]*X[1]**2-L[1]*X[0]*X[2]+  \
974                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
975            U[5] = T[5]+L[3]*X[0]*X[2]+L[4]*X[0]*X[1]-L[5]*X[0]**2-L[0]*X[2]*X[1]+  \
976                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
977            Umat = G2lat.U6toUij(U)
978            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
979            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
980        else:
981            Uout.append(['N',])
982    return Uout
983   
984def GetSHCoeff(pId,parmDict,SHkeys):
985    '''default doc string
986   
987    :param type name: description
988   
989    :returns: type name: description
990   
991    '''
992    SHCoeff = {}
993    for shkey in SHkeys:
994        shname = str(pId)+'::'+shkey
995        SHCoeff[shkey] = parmDict[shname]
996    return SHCoeff
997       
998def getMass(generalData):
999    '''Computes mass of unit cell contents
1000   
1001    :param dict generalData: The General dictionary in Phase
1002   
1003    :returns: float mass: Crystal unit cell mass in AMU.
1004   
1005    '''
1006    mass = 0.
1007    for i,elem in enumerate(generalData['AtomTypes']):
1008        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
1009    return max(mass,1.0)   
1010
1011def getDensity(generalData):
1012    '''calculate crystal structure density
1013   
1014    :param dict generalData: The General dictionary in Phase
1015   
1016    :returns: float density: crystal density in gm/cm^3
1017   
1018    '''
1019    mass = getMass(generalData)
1020    Volume = generalData['Cell'][7]
1021    density = mass/(0.6022137*Volume)
1022    return density,Volume/mass
1023   
1024def getWave(Parms):
1025    '''returns wavelength from Instrument parameters dictionary
1026   
1027    :param dict Parms: Instrument parameters;
1028        must contain:
1029        Lam: single wavelength
1030        or
1031        Lam1: Ka1 radiation wavelength
1032   
1033    :returns: float wave: wavelength
1034   
1035    '''
1036    try:
1037        return Parms['Lam'][1]
1038    except KeyError:
1039        return Parms['Lam1'][1]
1040       
1041def getMeanWave(Parms):
1042    '''returns mean wavelength from Instrument parameters dictionary
1043   
1044    :param dict Parms: Instrument parameters;
1045        must contain:
1046        Lam: single wavelength
1047        or
1048        Lam1,Lam2: Ka1,Ka2 radiation wavelength
1049        I(L2)/I(L1): Ka2/Ka1 ratio
1050   
1051    :returns: float wave: mean wavelength
1052   
1053    '''
1054    try:
1055        return Parms['Lam'][1]
1056    except KeyError:
1057        meanLam = (Parms['Lam1'][1]+Parms['I(L2)/I(L1)'][1]*Parms['Lam2'][1])/   \
1058            (1.+Parms['I(L2)/I(L1)'][1])
1059        return meanLam
1060   
1061       
1062def El2Mass(Elements):
1063    '''compute molecular weight from Elements
1064   
1065    :param dict Elements: elements in molecular formula;
1066        each must contain
1067        Num: number of atoms in formula
1068        Mass: at. wt.
1069   
1070    :returns: float mass: molecular weight.
1071   
1072    '''
1073    mass = 0
1074    for El in Elements:
1075        mass += Elements[El]['Num']*Elements[El]['Mass']
1076    return mass
1077       
1078def Den2Vol(Elements,density):
1079    '''converts density to molecular volume
1080   
1081    :param dict Elements: elements in molecular formula;
1082        each must contain
1083        Num: number of atoms in formula
1084        Mass: at. wt.
1085    :param float density: material density in gm/cm^3
1086   
1087    :returns: float volume: molecular volume in A^3
1088   
1089    '''
1090    return El2Mass(Elements)/(density*0.6022137)
1091   
1092def Vol2Den(Elements,volume):
1093    '''converts volume to density
1094   
1095    :param dict Elements: elements in molecular formula;
1096        each must contain
1097        Num: number of atoms in formula
1098        Mass: at. wt.
1099    :param float volume: molecular volume in A^3
1100   
1101    :returns: float density: material density in gm/cm^3
1102   
1103    '''
1104    return El2Mass(Elements)/(volume*0.6022137)
1105   
1106def El2EstVol(Elements):
1107    '''Estimate volume from molecular formula; assumes atom volume = 10A^3
1108   
1109    :param dict Elements: elements in molecular formula;
1110        each must contain
1111        Num: number of atoms in formula
1112   
1113    :returns: float volume: estimate of molecular volume in A^3
1114   
1115    '''
1116    vol = 0
1117    for El in Elements:
1118        vol += 10.*Elements[El]['Num']
1119    return vol
1120   
1121def XScattDen(Elements,vol,wave=0.):
1122    '''Estimate X-ray scattering density from molecular formula & volume;
1123    ignores valence, but includes anomalous effects
1124   
1125    :param dict Elements: elements in molecular formula;
1126        each element must contain
1127        Num: number of atoms in formula
1128        Z: atomic number
1129    :param float vol: molecular volume in A^3
1130    :param float wave: optional wavelength in A
1131   
1132    :returns: float rho: scattering density in 10^10cm^-2;
1133        if wave > 0 the includes f' contribution
1134    :returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
1135    :returns: float fpp: if wave>0 f" in 10^10cm^-2; otherwise 0
1136   
1137    '''
1138    rho = 0
1139    mu = 0
1140    fpp = 0
1141    if wave: 
1142        Xanom = XAnomAbs(Elements,wave)
1143    for El in Elements:
1144        f0 = Elements[El]['Z']
1145        if wave:
1146            f0 += Xanom[El][0]
1147            fpp += Xanom[El][1]*Elements[El]['Num']
1148            mu += Xanom[El][2]*Elements[El]['Num']
1149        rho += Elements[El]['Num']*f0
1150    return 28.179*rho/vol,0.1*mu/vol,28.179*fpp/vol
1151   
1152def NCScattDen(Elements,vol,wave=0.):
1153    '''Estimate neutron scattering density from molecular formula & volume;
1154    ignores valence, but includes anomalous effects
1155   
1156    :param dict Elements: elements in molecular formula;
1157        each element must contain
1158        Num: number of atoms in formula
1159        Z: atomic number
1160    :param float vol: molecular volume in A^3
1161    :param float wave: optional wavelength in A
1162   
1163    :returns: float rho: scattering density in 10^10cm^-2;
1164        if wave > 0 the includes f' contribution
1165    :returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
1166    :returns: float fpp: if wave>0 f" in 10^10cm^-2; otherwise 0
1167   
1168    '''
1169    rho = 0
1170    mu = 0
1171    bpp = 0
1172    for El in Elements:
1173        isotope = Elements[El]['Isotope']
1174        b0 = Elements[El]['Isotopes'][isotope]['SL'][0]
1175        mu += Elements[El]['Isotopes'][isotope].get('SA',0.)*Elements[El]['Num']
1176        if wave and 'BW-LS' in Elements[El]['Isotopes'][isotope]:
1177            Re,Im,E0,gam,A,E1,B,E2 = Elements[El]['Isotopes'][isotope]['BW-LS'][1:]
1178            Emev = 81.80703/wave**2
1179            T0 = Emev-E0
1180            T1 = Emev-E1
1181            T2 = Emev-E2
1182            D0 = T0**2+gam**2
1183            D1 = T1**2+gam**2
1184            D2 = T2**2+gam**2
1185            b0 += Re*(T0/D0+A*T1/D1+B*T2/D2)
1186            bpp += Im*(1/D0+A/D1+B/D2)
1187        else:
1188            bpp += Elements[El]['Isotopes'][isotope]['SL'][1]
1189        rho += Elements[El]['Num']*b0
1190    if wave: mu *= wave
1191    return 100.*rho/vol,mu/vol,100.*bpp/vol
1192   
1193def wavekE(wavekE):
1194    '''Convert wavelength to energy & vise versa
1195   
1196    :param float waveKe:wavelength in A or energy in kE
1197   
1198    :returns float waveKe:the other one
1199   
1200    '''
1201    return 12.397639/wavekE
1202       
1203def XAnomAbs(Elements,wave):
1204    kE = wavekE(wave)
1205    Xanom = {}
1206    for El in Elements:
1207        Orbs = G2el.GetXsectionCoeff(El)
1208        Xanom[El] = G2el.FPcalc(Orbs, kE)
1209    return Xanom        #f',f", mu
1210   
1211################################################################################
1212#### Modulation math
1213################################################################################
1214
1215def makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast):
1216    '''
1217    waveTypes: array nAtoms: 'Fourier','ZigZag' or 'Block'
1218    FSSdata: array 2 x atoms x waves    (sin,cos terms)
1219    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1220    USSdata: array 2x6 x atoms X waves (sin,cos terms)
1221    Mast: array orthogonalization matrix for Uij
1222    '''
1223    ngl = 32
1224    glTau,glWt = pwd.pygauleg(0.,1.,ngl)         #get Gauss-Legendre intervals & weights
1225    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
1226    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
1227    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
1228    Bf = np.array(FSSdata[1]).T    #cos frac mods...
1229    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods as betaij
1230    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods as betaij
1231    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] 
1232    if nWaves[0]:
1233        tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau  #Fwaves x ngl
1234        FmodA = Af[:,:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X ngl
1235        FmodB = Bf[:,:,nxs]*np.cos(twopi*tauF[nxs,:,:])
1236        Fmod = np.sum(1.0+FmodA+FmodB,axis=1)             #atoms X ngl; sum waves
1237    else:
1238        Fmod = 1.0           
1239    XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1240    XmodA = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1241    XmodB = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1242    for iatm in range(Ax.shape[0]):
1243        nx = 0
1244        if 'ZigZag' in waveTypes[iatm]:
1245            nx = 1
1246            Tmm = Ax[iatm][0][:2]                       
1247            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
1248            XmodZ[iatm][0] += posZigZag(glTau,Tmm,XYZmax).T
1249        elif 'Block' in waveTypes[iatm]:
1250            nx = 1
1251            Tmm = Ax[iatm][0][:2]                       
1252            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
1253            XmodZ[iatm][0] += posBlock(glTau,Tmm,XYZmax).T
1254        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
1255        if nx:   
1256            XmodA[iatm][:-nx] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1257            XmodB[iatm][:-nx] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
1258        else:
1259            XmodA[iatm] = Ax[iatm,:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1260            XmodB[iatm] = Bx[iatm,:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
1261    Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X 3 X ngl; sum waves
1262    Xmod = np.swapaxes(Xmod,1,2)
1263    if nWaves[2]:
1264        tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
1265        UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x ngl
1266        UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
1267        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x ngl; sum waves
1268    else:
1269        Umod = 1.0
1270#    GSASIIpath.IPyBreak()
1271    return ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt
1272       
1273def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1274    '''
1275    H: array nRefBlk x ops X hklt
1276    HP: array nRefBlk x ops X hklt proj to hkl
1277    Fmod: array 2 x atoms x waves    (sin,cos terms)
1278    Xmod: array atoms X 3 X ngl
1279    Umod: array atoms x 3x3 x ngl
1280    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1281    '''
1282   
1283    if nWaves[2]:
1284        if len(HP.shape) > 2:
1285            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1286        else:
1287            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1288    else:
1289        HbH = 1.0
1290    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
1291    if len(H.shape) > 2:
1292        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
1293        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1294    else:
1295        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1296        HdotXD = twopi*(HdotX+D[:,nxs,:])
1297    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1298    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
1299    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1300   
1301def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1302    '''
1303    H: array nRefBlk x tw x ops X hklt
1304    HP: array nRefBlk x tw x ops X hklt proj to hkl
1305    Fmod: array 2 x atoms x waves    (sin,cos terms)
1306    Xmod: array atoms X ngl X 3
1307    Umod: array atoms x ngl x 3x3
1308    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1309    '''
1310   
1311    if nWaves[2]:
1312        if len(HP.shape) > 3:   #Blocks of reflections
1313            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1314        else:   #single reflections
1315            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1316    else:
1317        HbH = 1.0
1318    HdotX = np.inner(HP,Xmod)                   #refBlk x tw x ops x atoms X ngl
1319    if len(H.shape) > 3:
1320        D = glTau*H[:,:,:,3:,nxs]              #m*e*tau; refBlk x tw x ops X ngl
1321        HdotXD = twopi*(HdotX+D[:,:,:,nxs,:])
1322    else:
1323        D = H*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1324        HdotXD = twopi*(HdotX+D[:,nxs,:])
1325    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1326    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
1327    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1328   
1329def makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast):
1330    '''
1331    FSSdata: array 2 x atoms x waves    (sin,cos terms)
1332    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1333    USSdata: array 2x6 x atoms X waves (sin,cos terms)
1334    Mast: array orthogonalization matrix for Uij
1335    '''
1336    glTau,glWt = pwd.pygauleg(0.,1.,ngl)         #get Gauss-Legendre intervals & weights
1337    waveShapes = [FSSdata.T.shape,XSSdata.T.shape,USSdata.T.shape]
1338    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
1339    Bf = np.array(FSSdata[1]).T    #cos frac mods...
1340    Ax = np.array(XSSdata[:3]).T   #...cos pos mods x waves x atoms
1341    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
1342    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
1343    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
1344    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] 
1345    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl
1346    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1347    ZtauXt = np.zeros((Ax.shape[0],2,3,ngl))               #atoms x Tminmax x 3 x ngl
1348    ZtauXx = np.zeros((Ax.shape[0],3,ngl))               #atoms x XYZmax x ngl
1349    for iatm in range(Ax.shape[0]):
1350        nx = 0
1351        if 'ZigZag' in waveTypes[iatm]:
1352            nx = 1
1353            Tmm = Ax[iatm][0][:2]                       
1354            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
1355            ZtauXt[iatm],ZtauXx[iatm] = posZigZagDerv(glTau,Tmm,XYZmax)
1356        elif 'Block' in waveTypes[iatm]:
1357            nx = 1
1358            Tmm = Ax[iatm][0][:2]                       
1359            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
1360            ZtauXt[iatm],ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax)
1361        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
1362        if nx:   
1363            StauX[iatm][:-nx] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
1364            CtauX[iatm][:-nx] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
1365        else:
1366            StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
1367            CtauX[iatm] = np.ones_like(Bx)[iatm,:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
1368#    GSASIIpath.IPyBreak()
1369    if nWaves[0]:
1370        tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau  #Fwaves x ngl
1371        StauF = np.ones_like(Af)[:,:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
1372        CtauF = np.ones_like(Bf)[:,:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
1373    else:
1374        StauF = 1.0
1375        CtauF = 1.0
1376    if nWaves[2]:
1377        tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
1378        StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodA/dAu
1379        CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodB/dBu
1380        UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x ngl
1381        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
1382#derivs need to be ops x atoms x waves x 6uij; ops x atoms x waves x ngl x 6uij before sum
1383        StauU = np.rollaxis(np.rollaxis(np.swapaxes(StauU,2,4),-1),-1)
1384        CtauU = np.rollaxis(np.rollaxis(np.swapaxes(CtauU,2,4),-1),-1)
1385    else:
1386        StauU = 1.0
1387        CtauU = 1.0
1388        UmodA = 0.
1389        UmodB = 0.
1390    return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
1391   
1392def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
1393    '''
1394    H: array ops X hklt proj to hkl
1395    HP: array ops X hklt proj to hkl
1396    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
1397    '''
1398   
1399    Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
1400    dGdMfC = np.zeros(Mf)
1401    dGdMfS = np.zeros(Mf)
1402    Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
1403    dGdMxC = np.zeros(Mx)
1404    dGdMxS = np.zeros(Mx)
1405    Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
1406    dGdMuC = np.zeros(Mu)
1407    dGdMuS = np.zeros(Mu)
1408   
1409    D = twopi*H[:,3][:,nxs]*glTau[nxs,:]              #m*e*tau; ops X ngl
1410    HdotX = twopi*np.inner(HP,Xmod)        #ops x atoms X ngl
1411    HdotXD = HdotX+D[:,nxs,:]
1412    if nWaves[2]:
1413        Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
1414        HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #ops x atoms x waves x ngl
1415        HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #ops x atoms x waves x ngl
1416        HbH = np.exp(-np.sum(HuH,axis=-2)) # ops x atoms x ngl; sum waves - OK vs Modulation version
1417        part1 = -np.exp(-HuH)*Fmod    #ops x atoms x waves x ngl
1418        dUdAu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
1419        dUdBu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
1420        dGdMuCa = np.sum(part1[:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1421        dGdMuCb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1422        dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
1423        dGdMuSa = np.sum(part1[:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1424        dGdMuSb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1425        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
1426    else:
1427        HbH = np.ones_like(HdotX)
1428    dHdXA = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
1429    dHdXB = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,:,:,:,:]    #ditto - cos waves
1430# ops x atoms x waves x 2xyz - real part - good
1431    dGdMxCa = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1432    dGdMxCb = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1433    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
1434# ops x atoms x waves x 2xyz - imag part - good
1435    dGdMxSa = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1436    dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1437    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
1438# ZigZag/Block waves - problems here?
1439    dHdXZt = -twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
1440    dHdXZx = twopi*HP[:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
1441    dGdMzCt = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1442    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
1443    dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
1444    dGdMzSt = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1445    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
1446    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
1447#    GSASIIpath.IPyBreak()
1448    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
1449   
1450def ModulationDerv2(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
1451    '''
1452    H: array refBlk x ops X hklt proj to hkl
1453    HP: array refBlk x ops X hklt proj to hkl
1454    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
1455    '''
1456   
1457    Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
1458    dGdMfC = np.zeros(Mf)
1459    dGdMfS = np.zeros(Mf)
1460    Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
1461    dGdMxC = np.zeros(Mx)
1462    dGdMxS = np.zeros(Mx)
1463    Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
1464    dGdMuC = np.zeros(Mu)
1465    dGdMuS = np.zeros(Mu)
1466   
1467    D = twopi*H[:,:,3,nxs]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
1468    HdotX = twopi*np.inner(HP,Xmod)        #refBlk x ops x atoms X ngl
1469    HdotXD = HdotX+D[:,:,nxs,:]
1470    if nWaves[2]:
1471        Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
1472        HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
1473        HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #refBlk x ops x atoms x waves x ngl
1474        HbH = np.exp(-np.sum(HuH,axis=-2)) #refBlk x ops x atoms x ngl; sum waves - OK vs Modulation version
1475        part1 = -np.exp(-HuH)*Fmod    #refBlk x ops x atoms x waves x ngl
1476        dUdAu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
1477        dUdBu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
1478        dGdMuCa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1479        dGdMuCb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1480        dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
1481        dGdMuSa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1482        dGdMuSb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1483        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
1484    else:
1485        HbH = np.ones_like(HdotX)
1486    dHdXA = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
1487    dHdXB = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,nxs,:,:,:,:]    #ditto - cos waves
1488# ops x atoms x waves x 2xyz - real part - good
1489    dGdMxCa = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1490    dGdMxCb = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1491    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
1492# ops x atoms x waves x 2xyz - imag part - good
1493    dGdMxSa = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
1494    dGdMxSb = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)   
1495    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
1496# ZigZag/Block waves - problems here?
1497    dHdXZt = -twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
1498    dHdXZx = twopi*HP[:,:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
1499    dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1500    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1501    dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
1502    dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2)
1503    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1504    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
1505#    GSASIIpath.IPyBreak()
1506    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
1507   
1508def posFourier(tau,psin,pcos):
1509    A = np.array([ps[:,nxs]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
1510    B = np.array([pc[:,nxs]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
1511    return np.sum(A,axis=0)+np.sum(B,axis=0)
1512   
1513def posZigZag(T,Tmm,Xmax):
1514    DT = Tmm[1]-Tmm[0]
1515    Su = 2.*Xmax/DT
1516    Sd = 2.*Xmax/(1.-DT)
1517    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax+Su*((t-Tmm[0])%1.),Xmax-Sd*((t-Tmm[1])%1.)) for t in T])
1518    return A
1519   
1520def posZigZagDerv(T,Tmm,Xmax):
1521    DT = Tmm[1]-Tmm[0]
1522    Su = 2.*Xmax/DT
1523    Sd = 2.*Xmax/(1.-DT)
1524    dAdT = np.zeros((2,3,len(T)))
1525    dAdT[0] = np.array([np.where(Tmm[0] < t <= Tmm[1],Su*(t-Tmm[0]-1)/DT,-Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
1526    dAdT[1] = np.array([np.where(Tmm[0] < t <= Tmm[1],-Su*(t-Tmm[0])/DT,Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
1527    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.+2.*(t-Tmm[0])/DT,1.-2.*(t-Tmm[1])%1./DT) for t in T])
1528    return dAdT,dAdX
1529
1530def posBlock(T,Tmm,Xmax):
1531    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax,Xmax) for t in T])
1532    return A
1533   
1534def posBlockDerv(T,Tmm,Xmax):
1535    dAdT = np.zeros((2,3,len(T)))
1536    ind = np.searchsorted(T,Tmm)
1537    dAdT[0,:,ind[0]] = -Xmax/len(T)
1538    dAdT[1,:,ind[1]] = Xmax/len(T)
1539    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t <= Tmm[1],-1.,1.) for t in T])  #OK
1540    return dAdT,dAdX
1541   
1542def fracCrenel(tau,Toff,Twid):
1543    Tau = (tau-Toff)%1.
1544    A = np.where(Tau<Twid,1.,0.)
1545    return A
1546   
1547def fracFourier(tau,fsin,fcos):
1548    A = np.array([fs[:,nxs]*np.sin(2.*np.pi*(i+1)*tau) for i,fs in enumerate(fsin)])
1549    B = np.array([fc[:,nxs]*np.cos(2.*np.pi*(i+1)*tau) for i,fc in enumerate(fcos)])
1550    return np.sum(A,axis=0)+np.sum(B,axis=0)
1551   
1552def ApplyModulation(data,tau):
1553    '''Applies modulation to drawing atom positions & Uijs for given tau
1554    '''
1555    generalData = data['General']
1556    SGData = generalData['SGData']
1557    SSGData = generalData['SSGData']
1558    cx,ct,cs,cia = generalData['AtomPtrs']
1559    drawingData = data['Drawing']
1560    dcx,dct,dcs,dci = drawingData['atomPtrs']
1561    atoms = data['Atoms']
1562    drawAtoms = drawingData['Atoms']
1563    Fade = np.zeros(len(drawAtoms))
1564    for atom in atoms:   
1565        atxyz = G2spc.MoveToUnitCell(np.array(atom[cx:cx+3]))[0]
1566        atuij = np.array(atom[cia+2:cia+8])
1567        waveType = atom[-1]['SS1']['waveType']
1568        Sfrac = atom[-1]['SS1']['Sfrac']
1569        Spos = atom[-1]['SS1']['Spos']
1570        Sadp = atom[-1]['SS1']['Sadp']
1571        indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
1572        for ind in indx:
1573            drawatom = drawAtoms[ind]
1574            opr = drawatom[dcs-1]
1575            sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData)
1576            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz)
1577            tauT *= icent       #invert wave on -1
1578            wave = np.zeros(3)
1579            uwave = np.zeros(6)
1580            #how do I handle Sfrac? - fade the atoms?
1581            if len(Sfrac):
1582                scof = []
1583                ccof = []
1584                for i,sfrac in enumerate(Sfrac):
1585                    if not i and 'Crenel' in waveType:
1586                        Fade[ind] += fracCrenel(tauT,sfrac[0][0],sfrac[0][1])
1587                    else:
1588                        scof.append(sfrac[0][0])
1589                        ccof.append(sfrac[0][1])
1590                    if len(scof):
1591                        Fade[ind] += np.sum(fracFourier(tauT,scof,ccof))                           
1592            if len(Spos):
1593                scof = []
1594                ccof = []
1595                for i,spos in enumerate(Spos):
1596                    if waveType in ['ZigZag','Block'] and not i:
1597                        Tminmax = spos[0][:2]
1598                        XYZmax = np.array(spos[0][2:])
1599                        if waveType == 'Block':
1600                            wave = np.array(posBlock([tauT,],Tminmax,XYZmax))[0]
1601                        elif waveType == 'ZigZag':
1602                            wave = np.array(posZigZag([tauT,],Tminmax,XYZmax))[0]
1603                    else:
1604                        scof.append(spos[0][:3])
1605                        ccof.append(spos[0][3:])
1606                if len(scof):
1607                    wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1608            if len(Sadp):
1609                scof = []
1610                ccof = []
1611                for i,sadp in enumerate(Sadp):
1612                    scof.append(sadp[0][:6])
1613                    ccof.append(sadp[0][6:])
1614                uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1615            if atom[cia] == 'A':                   
1616                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
1617                drawatom[dcx:dcx+3] = X
1618                drawatom[dci-6:dci] = U
1619            else:
1620                X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
1621                drawatom[dcx:dcx+3] = X
1622    return drawAtoms,Fade
1623   
1624# gauleg.py Gauss Legendre numerical quadrature, x and w computation
1625# integrate from a to b using n evaluations of the function f(x) 
1626# usage: from gauleg import gaulegf         
1627#        x,w = gaulegf( a, b, n)                               
1628#        area = 0.0                                           
1629#        for i in range(1,n+1):          #  yes, 1..n                   
1630#          area += w[i]*f(x[i])                                   
1631
1632def gaulegf(a, b, n):
1633  x = range(n+1) # x[0] unused
1634  w = range(n+1) # w[0] unused
1635  eps = 3.0E-14
1636  m = (n+1)/2
1637  xm = 0.5*(b+a)
1638  xl = 0.5*(b-a)
1639  for i in range(1,m+1):
1640    z = math.cos(3.141592654*(i-0.25)/(n+0.5))
1641    while True:
1642      p1 = 1.0
1643      p2 = 0.0
1644      for j in range(1,n+1):
1645        p3 = p2
1646        p2 = p1
1647        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
1648
1649      pp = n*(z*p1-p2)/(z*z-1.0)
1650      z1 = z
1651      z = z1 - p1/pp
1652      if abs(z-z1) <= eps:
1653            break
1654
1655    x[i] = xm - xl*z
1656    x[n+1-i] = xm + xl*z
1657    w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
1658    w[n+1-i] = w[i]
1659  return np.array(x), np.array(w)
1660# end gaulegf
1661   
1662   
1663def BessJn(nmax,x):
1664    ''' compute Bessel function J(n,x) from scipy routine & recurrance relation
1665    returns sequence of J(n,x) for n in range [-nmax...0...nmax]
1666   
1667    :param  integer nmax: maximul order for Jn(x)
1668    :param  float x: argument for Jn(x)
1669   
1670    :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
1671   
1672    '''
1673    import scipy.special as sp
1674    bessJn = np.zeros(2*nmax+1)
1675    bessJn[nmax] = sp.j0(x)
1676    bessJn[nmax+1] = sp.j1(x)
1677    bessJn[nmax-1] = -bessJn[nmax+1]
1678    for i in range(2,nmax+1):
1679        bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
1680        bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i
1681    return bessJn
1682   
1683def BessIn(nmax,x):
1684    ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
1685    returns sequence of I(n,x) for n in range [-nmax...0...nmax]
1686   
1687    :param  integer nmax: maximul order for In(x)
1688    :param  float x: argument for In(x)
1689   
1690    :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
1691   
1692    '''
1693    import scipy.special as sp
1694    bessIn = np.zeros(2*nmax+1)
1695    bessIn[nmax] = sp.i0(x)
1696    bessIn[nmax+1] = sp.i1(x)
1697    bessIn[nmax-1] = bessIn[nmax+1]
1698    for i in range(2,nmax+1):
1699        bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x
1700        bessIn[nmax-i] = bessIn[i+nmax]
1701    return bessIn
1702       
1703   
1704################################################################################
1705##### distance, angle, planes, torsion stuff
1706################################################################################
1707
1708def CalcDist(distance_dict, distance_atoms, parmDict):
1709    if not len(parmDict):
1710        return 0.
1711    pId = distance_dict['pId']
1712    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1713    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1714    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1715    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1716    inv = 1
1717    symNo = distance_dict['symNo']
1718    if symNo < 0:
1719        inv = -1
1720        symNo *= -1
1721    cen = symNo/100
1722    op = symNo%100-1
1723    M,T = distance_dict['SGData']['SGOps'][op]
1724    D = T*inv+distance_dict['SGData']['SGCen'][cen]
1725    D += distance_dict['cellNo']
1726    Txyz = np.inner(M*inv,Txyz)+D
1727    dist = np.sqrt(np.sum(np.inner(Amat,(Txyz-Oxyz))**2))
1728#    GSASIIpath.IPyBreak()
1729    return dist   
1730   
1731def CalcDistDeriv(distance_dict, distance_atoms, parmDict):
1732    if not len(parmDict):
1733        return None
1734    pId = distance_dict['pId']
1735    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1736    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1737    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1738    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1739    symNo = distance_dict['symNo']
1740    Tunit = distance_dict['cellNo']
1741    SGData = distance_dict['SGData']   
1742    deriv = getDistDerv(Oxyz,Txyz,Amat,Tunit,symNo,SGData)
1743    return deriv
1744   
1745def CalcAngle(angle_dict, angle_atoms, parmDict):
1746    if not len(parmDict):
1747        return 0.
1748    pId = angle_dict['pId']
1749    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1750    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1751    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1752    Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1753    Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1754    ABxyz = [Axyz,Bxyz]
1755    symNo = angle_dict['symNo']
1756    vec = np.zeros((2,3))
1757    for i in range(2):
1758        inv = 1
1759        if symNo[i] < 0:
1760            inv = -1
1761        cen = inv*symNo[i]/100
1762        op = inv*symNo[i]%100-1
1763        M,T = angle_dict['SGData']['SGOps'][op]
1764        D = T*inv+angle_dict['SGData']['SGCen'][cen]
1765        D += angle_dict['cellNo'][i]
1766        ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
1767        vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
1768        dist = np.sqrt(np.sum(vec[i]**2))
1769        if not dist:
1770            return 0.
1771        vec[i] /= dist
1772    angle = acosd(np.sum(vec[0]*vec[1]))
1773#    GSASIIpath.IPyBreak()
1774    return angle
1775
1776def CalcAngleDeriv(angle_dict, angle_atoms, parmDict):
1777    if not len(parmDict):
1778        return None
1779    pId = angle_dict['pId']
1780    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1781    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1782    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1783    Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1784    Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1785    symNo = angle_dict['symNo']
1786    Tunit = angle_dict['cellNo']
1787    SGData = angle_dict['SGData']   
1788    deriv = getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData)
1789    return deriv
1790
1791def getSyXYZ(XYZ,ops,SGData):
1792    '''default doc
1793   
1794    :param type name: description
1795   
1796    :returns: type name: description
1797   
1798    '''
1799    XYZout = np.zeros_like(XYZ)
1800    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
1801        if op == '1':
1802            XYZout[i] = xyz
1803        else:
1804            oprs = op.split('+')
1805            unit = [0,0,0]
1806            if len(oprs)>1:
1807                unit = np.array(list(eval(oprs[1])))
1808            syop =int(oprs[0])
1809            inv = syop/abs(syop)
1810            syop *= inv
1811            cent = syop/100
1812            syop %= 100
1813            syop -= 1
1814            M,T = SGData['SGOps'][syop]
1815            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
1816    return XYZout
1817   
1818def getRestDist(XYZ,Amat):
1819    '''default doc string
1820   
1821    :param type name: description
1822   
1823    :returns: type name: description
1824   
1825    '''
1826    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1827   
1828def getRestDeriv(Func,XYZ,Amat,ops,SGData):
1829    '''default doc string
1830   
1831    :param type name: description
1832   
1833    :returns: type name: description
1834   
1835    '''
1836    deriv = np.zeros((len(XYZ),3))
1837    dx = 0.00001
1838    for j,xyz in enumerate(XYZ):
1839        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1840            XYZ[j] -= x
1841            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1842            XYZ[j] += 2*x
1843            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1844            XYZ[j] -= x
1845            deriv[j][i] = (d1-d2)/(2*dx)
1846    return deriv.flatten()
1847
1848def getRestAngle(XYZ,Amat):
1849    '''default doc string
1850   
1851    :param type name: description
1852   
1853    :returns: type name: description
1854   
1855    '''
1856   
1857    def calcVec(Ox,Tx,Amat):
1858        return np.inner(Amat,(Tx-Ox))
1859
1860    VecA = calcVec(XYZ[1],XYZ[0],Amat)
1861    VecA /= np.sqrt(np.sum(VecA**2))
1862    VecB = calcVec(XYZ[1],XYZ[2],Amat)
1863    VecB /= np.sqrt(np.sum(VecB**2))
1864    edge = VecB-VecA
1865    edge = np.sum(edge**2)
1866    angle = (2.-edge)/2.
1867    angle = max(angle,-1.)
1868    return acosd(angle)
1869   
1870def getRestPlane(XYZ,Amat):
1871    '''default doc string
1872   
1873    :param type name: description
1874   
1875    :returns: type name: description
1876   
1877    '''
1878    sumXYZ = np.zeros(3)
1879    for xyz in XYZ:
1880        sumXYZ += xyz
1881    sumXYZ /= len(XYZ)
1882    XYZ = np.array(XYZ)-sumXYZ
1883    XYZ = np.inner(Amat,XYZ).T
1884    Zmat = np.zeros((3,3))
1885    for i,xyz in enumerate(XYZ):
1886        Zmat += np.outer(xyz.T,xyz)
1887    Evec,Emat = nl.eig(Zmat)
1888    Evec = np.sqrt(Evec)/(len(XYZ)-3)
1889    Order = np.argsort(Evec)
1890    return Evec[Order[0]]
1891   
1892def getRestChiral(XYZ,Amat):   
1893    '''default doc string
1894   
1895    :param type name: description
1896   
1897    :returns: type name: description
1898   
1899    '''
1900    VecA = np.empty((3,3))   
1901    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1902    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
1903    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
1904    return nl.det(VecA)
1905   
1906def getRestTorsion(XYZ,Amat):
1907    '''default doc string
1908   
1909    :param type name: description
1910   
1911    :returns: type name: description
1912   
1913    '''
1914    VecA = np.empty((3,3))
1915    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1916    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
1917    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
1918    D = nl.det(VecA)
1919    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
1920    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
1921    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
1922    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
1923    Ang = 1.0
1924    if abs(P12) < 1.0 and abs(P23) < 1.0:
1925        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
1926    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
1927    return TOR
1928   
1929def calcTorsionEnergy(TOR,Coeff=[]):
1930    '''default doc string
1931   
1932    :param type name: description
1933   
1934    :returns: type name: description
1935   
1936    '''
1937    sum = 0.
1938    if len(Coeff):
1939        cof = np.reshape(Coeff,(3,3)).T
1940        delt = TOR-cof[1]
1941        delt = np.where(delt<-180.,delt+360.,delt)
1942        delt = np.where(delt>180.,delt-360.,delt)
1943        term = -cof[2]*delt**2
1944        val = cof[0]*np.exp(term/1000.0)
1945        pMax = cof[0][np.argmin(val)]
1946        Eval = np.sum(val)
1947        sum = Eval-pMax
1948    return sum,Eval
1949
1950def getTorsionDeriv(XYZ,Amat,Coeff):
1951    '''default doc string
1952   
1953    :param type name: description
1954   
1955    :returns: type name: description
1956   
1957    '''
1958    deriv = np.zeros((len(XYZ),3))
1959    dx = 0.00001
1960    for j,xyz in enumerate(XYZ):
1961        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1962            XYZ[j] -= x
1963            tor = getRestTorsion(XYZ,Amat)
1964            p1,d1 = calcTorsionEnergy(tor,Coeff)
1965            XYZ[j] += 2*x
1966            tor = getRestTorsion(XYZ,Amat)
1967            p2,d2 = calcTorsionEnergy(tor,Coeff)           
1968            XYZ[j] -= x
1969            deriv[j][i] = (p2-p1)/(2*dx)
1970    return deriv.flatten()
1971
1972def getRestRama(XYZ,Amat):
1973    '''Computes a pair of torsion angles in a 5 atom string
1974   
1975    :param nparray XYZ: crystallographic coordinates of 5 atoms
1976    :param nparray Amat: crystal to cartesian transformation matrix
1977   
1978    :returns: list (phi,psi) two torsion angles in degrees
1979   
1980    '''
1981    phi = getRestTorsion(XYZ[:5],Amat)
1982    psi = getRestTorsion(XYZ[1:],Amat)
1983    return phi,psi
1984   
1985def calcRamaEnergy(phi,psi,Coeff=[]):
1986    '''Computes pseudo potential energy from a pair of torsion angles and a
1987    numerical description of the potential energy surface. Used to create
1988    penalty function in LS refinement:     
1989    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
1990
1991    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
1992   
1993    :param float phi: first torsion angle (:math:`\\phi`)
1994    :param float psi: second torsion angle (:math:`\\psi`)
1995    :param list Coeff: pseudo potential coefficients
1996   
1997    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
1998      sum is used for penalty function.
1999   
2000    '''
2001    sum = 0.
2002    Eval = 0.
2003    if len(Coeff):
2004        cof = Coeff.T
2005        dPhi = phi-cof[1]
2006        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
2007        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
2008        dPsi = psi-cof[2]
2009        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
2010        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
2011        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
2012        val = cof[0]*np.exp(val/1000.)
2013        pMax = cof[0][np.argmin(val)]
2014        Eval = np.sum(val)
2015        sum = Eval-pMax
2016    return sum,Eval
2017
2018def getRamaDeriv(XYZ,Amat,Coeff):
2019    '''Computes numerical derivatives of torsion angle pair pseudo potential
2020    with respect of crystallographic atom coordinates of the 5 atom sequence
2021   
2022    :param nparray XYZ: crystallographic coordinates of 5 atoms
2023    :param nparray Amat: crystal to cartesian transformation matrix
2024    :param list Coeff: pseudo potential coefficients
2025   
2026    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
2027     crystallographic xyz coordinates.
2028   
2029    '''
2030    deriv = np.zeros((len(XYZ),3))
2031    dx = 0.00001
2032    for j,xyz in enumerate(XYZ):
2033        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2034            XYZ[j] -= x
2035            phi,psi = getRestRama(XYZ,Amat)
2036            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
2037            XYZ[j] += 2*x
2038            phi,psi = getRestRama(XYZ,Amat)
2039            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
2040            XYZ[j] -= x
2041            deriv[j][i] = (p2-p1)/(2*dx)
2042    return deriv.flatten()
2043
2044def getRestPolefig(ODFln,SamSym,Grid):
2045    '''default doc string
2046   
2047    :param type name: description
2048   
2049    :returns: type name: description
2050   
2051    '''
2052    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
2053    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
2054    R = np.where(R <= 1.,2.*atand(R),0.0)
2055    Z = np.zeros_like(R)
2056    Z = G2lat.polfcal(ODFln,SamSym,R,P)
2057    Z = np.reshape(Z,(Grid,Grid))
2058    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
2059
2060def getRestPolefigDerv(HKL,Grid,SHCoeff):
2061    '''default doc string
2062   
2063    :param type name: description
2064   
2065    :returns: type name: description
2066   
2067    '''
2068    pass
2069       
2070def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
2071    '''default doc string
2072   
2073    :param type name: description
2074   
2075    :returns: type name: description
2076   
2077    '''
2078    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
2079        TxT = inv*(np.inner(M,Tx)+T)+C+U
2080        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
2081       
2082    inv = Top/abs(Top)
2083    cent = abs(Top)/100
2084    op = abs(Top)%100-1
2085    M,T = SGData['SGOps'][op]
2086    C = SGData['SGCen'][cent]
2087    dx = .00001
2088    deriv = np.zeros(6)
2089    for i in [0,1,2]:
2090        Oxyz[i] -= dx
2091        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2092        Oxyz[i] += 2*dx
2093        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2094        Oxyz[i] -= dx
2095        Txyz[i] -= dx
2096        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2097        Txyz[i] += 2*dx
2098        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2099        Txyz[i] -= dx
2100    return deriv
2101   
2102def getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData):
2103   
2104    def calcAngle(Oxyz,ABxyz,Amat,Tunit,symNo,SGData):
2105        vec = np.zeros((2,3))
2106        for i in range(2):
2107            inv = 1
2108            if symNo[i] < 0:
2109                inv = -1
2110            cen = inv*symNo[i]/100
2111            op = inv*symNo[i]%100-1
2112            M,T = SGData['SGOps'][op]
2113            D = T*inv+SGData['SGCen'][cen]
2114            D += Tunit[i]
2115            ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
2116            vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
2117            dist = np.sqrt(np.sum(vec[i]**2))
2118            if not dist:
2119                return 0.
2120            vec[i] /= dist
2121        angle = acosd(np.sum(vec[0]*vec[1]))
2122    #    GSASIIpath.IPyBreak()
2123        return angle
2124       
2125    dx = .00001
2126    deriv = np.zeros(9)
2127    for i in [0,1,2]:
2128        Oxyz[i] -= dx
2129        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2130        Oxyz[i] += 2*dx
2131        deriv[i] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2132        Oxyz[i] -= dx
2133        Axyz[i] -= dx
2134        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2135        Axyz[i] += 2*dx
2136        deriv[i+3] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2137        Axyz[i] -= dx
2138        Bxyz[i] -= dx
2139        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2140        Bxyz[i] += 2*dx
2141        deriv[i+6] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2142        Bxyz[i] -= dx
2143    return deriv
2144   
2145def getAngSig(VA,VB,Amat,SGData,covData={}):
2146    '''default doc string
2147   
2148    :param type name: description
2149   
2150    :returns: type name: description
2151   
2152    '''
2153    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
2154        TxT = inv*(np.inner(M,Tx)+T)+C+U
2155        return np.inner(Amat,(TxT-Ox))
2156       
2157    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
2158        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
2159        VecA /= np.sqrt(np.sum(VecA**2))
2160        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
2161        VecB /= np.sqrt(np.sum(VecB**2))
2162        edge = VecB-VecA
2163        edge = np.sum(edge**2)
2164        angle = (2.-edge)/2.
2165        angle = max(angle,-1.)
2166        return acosd(angle)
2167       
2168    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
2169    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
2170    invA = invB = 1
2171    invA = TopA/abs(TopA)
2172    invB = TopB/abs(TopB)
2173    centA = abs(TopA)/100
2174    centB = abs(TopB)/100
2175    opA = abs(TopA)%100-1
2176    opB = abs(TopB)%100-1
2177    MA,TA = SGData['SGOps'][opA]
2178    MB,TB = SGData['SGOps'][opB]
2179    CA = SGData['SGCen'][centA]
2180    CB = SGData['SGCen'][centB]
2181    if 'covMatrix' in covData:
2182        covMatrix = covData['covMatrix']
2183        varyList = covData['varyList']
2184        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
2185        dx = .00001
2186        dadx = np.zeros(9)
2187        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2188        for i in [0,1,2]:
2189            OxA[i] -= dx
2190            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2191            OxA[i] += 2*dx
2192            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2193            OxA[i] -= dx
2194           
2195            TxA[i] -= dx
2196            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2197            TxA[i] += 2*dx
2198            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2199            TxA[i] -= dx
2200           
2201            TxB[i] -= dx
2202            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2203            TxB[i] += 2*dx
2204            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2205            TxB[i] -= dx
2206           
2207        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2208        if sigAng < 0.01:
2209            sigAng = 0.0
2210        return Ang,sigAng
2211    else:
2212        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
2213
2214def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
2215    '''default doc string
2216   
2217    :param type name: description
2218   
2219    :returns: type name: description
2220   
2221    '''
2222    def calcDist(Atoms,SyOps,Amat):
2223        XYZ = []
2224        for i,atom in enumerate(Atoms):
2225            Inv,M,T,C,U = SyOps[i]
2226            XYZ.append(np.array(atom[1:4]))
2227            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2228            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2229        V1 = XYZ[1]-XYZ[0]
2230        return np.sqrt(np.sum(V1**2))
2231       
2232    SyOps = []
2233    names = []
2234    for i,atom in enumerate(Oatoms):
2235        names += atom[-1]
2236        Op,unit = Atoms[i][-1]
2237        inv = Op/abs(Op)
2238        m,t = SGData['SGOps'][abs(Op)%100-1]
2239        c = SGData['SGCen'][abs(Op)/100]
2240        SyOps.append([inv,m,t,c,unit])
2241    Dist = calcDist(Oatoms,SyOps,Amat)
2242   
2243    sig = -0.001
2244    if 'covMatrix' in covData:
2245        dx = .00001
2246        dadx = np.zeros(6)
2247        for i in range(6):
2248            ia = i/3
2249            ix = i%3
2250            Oatoms[ia][ix+1] += dx
2251            a0 = calcDist(Oatoms,SyOps,Amat)
2252            Oatoms[ia][ix+1] -= 2*dx
2253            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2254        covMatrix = covData['covMatrix']
2255        varyList = covData['varyList']
2256        DistVcov = getVCov(names,varyList,covMatrix)
2257        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
2258        if sig < 0.001:
2259            sig = -0.001
2260   
2261    return Dist,sig
2262
2263def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
2264    '''default doc string
2265   
2266    :param type name: description
2267   
2268    :returns: type name: description
2269   
2270    '''
2271
2272    def calcAngle(Atoms,SyOps,Amat):
2273        XYZ = []
2274        for i,atom in enumerate(Atoms):
2275            Inv,M,T,C,U = SyOps[i]
2276            XYZ.append(np.array(atom[1:4]))
2277            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2278            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2279        V1 = XYZ[1]-XYZ[0]
2280        V1 /= np.sqrt(np.sum(V1**2))
2281        V2 = XYZ[1]-XYZ[2]
2282        V2 /= np.sqrt(np.sum(V2**2))
2283        V3 = V2-V1
2284        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2285        return acosd(cang)
2286
2287    SyOps = []
2288    names = []
2289    for i,atom in enumerate(Oatoms):
2290        names += atom[-1]
2291        Op,unit = Atoms[i][-1]
2292        inv = Op/abs(Op)
2293        m,t = SGData['SGOps'][abs(Op)%100-1]
2294        c = SGData['SGCen'][abs(Op)/100]
2295        SyOps.append([inv,m,t,c,unit])
2296    Angle = calcAngle(Oatoms,SyOps,Amat)
2297   
2298    sig = -0.01
2299    if 'covMatrix' in covData:
2300        dx = .00001
2301        dadx = np.zeros(9)
2302        for i in range(9):
2303            ia = i/3
2304            ix = i%3
2305            Oatoms[ia][ix+1] += dx
2306            a0 = calcAngle(Oatoms,SyOps,Amat)
2307            Oatoms[ia][ix+1] -= 2*dx
2308            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2309        covMatrix = covData['covMatrix']
2310        varyList = covData['varyList']
2311        AngVcov = getVCov(names,varyList,covMatrix)
2312        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2313        if sig < 0.01:
2314            sig = -0.01
2315   
2316    return Angle,sig
2317
2318def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
2319    '''default doc string
2320   
2321    :param type name: description
2322   
2323    :returns: type name: description
2324   
2325    '''
2326
2327    def calcTorsion(Atoms,SyOps,Amat):
2328       
2329        XYZ = []
2330        for i,atom in enumerate(Atoms):
2331            Inv,M,T,C,U = SyOps[i]
2332            XYZ.append(np.array(atom[1:4]))
2333            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2334            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2335        V1 = XYZ[1]-XYZ[0]
2336        V2 = XYZ[2]-XYZ[1]
2337        V3 = XYZ[3]-XYZ[2]
2338        V1 /= np.sqrt(np.sum(V1**2))
2339        V2 /= np.sqrt(np.sum(V2**2))
2340        V3 /= np.sqrt(np.sum(V3**2))
2341        M = np.array([V1,V2,V3])
2342        D = nl.det(M)
2343        P12 = np.dot(V1,V2)
2344        P13 = np.dot(V1,V3)
2345        P23 = np.dot(V2,V3)
2346        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2347        return Tors
2348           
2349    SyOps = []
2350    names = []
2351    for i,atom in enumerate(Oatoms):
2352        names += atom[-1]
2353        Op,unit = Atoms[i][-1]
2354        inv = Op/abs(Op)
2355        m,t = SGData['SGOps'][abs(Op)%100-1]
2356        c = SGData['SGCen'][abs(Op)/100]
2357        SyOps.append([inv,m,t,c,unit])
2358    Tors = calcTorsion(Oatoms,SyOps,Amat)
2359   
2360    sig = -0.01
2361    if 'covMatrix' in covData:
2362        dx = .00001
2363        dadx = np.zeros(12)
2364        for i in range(12):
2365            ia = i/3
2366            ix = i%3
2367            Oatoms[ia][ix+1] -= dx
2368            a0 = calcTorsion(Oatoms,SyOps,Amat)
2369            Oatoms[ia][ix+1] += 2*dx
2370            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2371            Oatoms[ia][ix+1] -= dx           
2372        covMatrix = covData['covMatrix']
2373        varyList = covData['varyList']
2374        TorVcov = getVCov(names,varyList,covMatrix)
2375        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
2376        if sig < 0.01:
2377            sig = -0.01
2378   
2379    return Tors,sig
2380       
2381def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
2382    '''default doc string
2383   
2384    :param type name: description
2385   
2386    :returns: type name: description
2387   
2388    '''
2389
2390    def calcDist(Atoms,SyOps,Amat):
2391        XYZ = []
2392        for i,atom in enumerate(Atoms):
2393            Inv,M,T,C,U = SyOps[i]
2394            XYZ.append(np.array(atom[1:4]))
2395            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2396            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2397        V1 = XYZ[1]-XYZ[0]
2398        return np.sqrt(np.sum(V1**2))
2399       
2400    def calcAngle(Atoms,SyOps,Amat):
2401        XYZ = []
2402        for i,atom in enumerate(Atoms):
2403            Inv,M,T,C,U = SyOps[i]
2404            XYZ.append(np.array(atom[1:4]))
2405            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2406            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2407        V1 = XYZ[1]-XYZ[0]
2408        V1 /= np.sqrt(np.sum(V1**2))
2409        V2 = XYZ[1]-XYZ[2]
2410        V2 /= np.sqrt(np.sum(V2**2))
2411        V3 = V2-V1
2412        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2413        return acosd(cang)
2414
2415    def calcTorsion(Atoms,SyOps,Amat):
2416       
2417        XYZ = []
2418        for i,atom in enumerate(Atoms):
2419            Inv,M,T,C,U = SyOps[i]
2420            XYZ.append(np.array(atom[1:4]))
2421            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2422            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2423        V1 = XYZ[1]-XYZ[0]
2424        V2 = XYZ[2]-XYZ[1]
2425        V3 = XYZ[3]-XYZ[2]
2426        V1 /= np.sqrt(np.sum(V1**2))
2427        V2 /= np.sqrt(np.sum(V2**2))
2428        V3 /= np.sqrt(np.sum(V3**2))
2429        M = np.array([V1,V2,V3])
2430        D = nl.det(M)
2431        P12 = np.dot(V1,V2)
2432        P13 = np.dot(V1,V3)
2433        P23 = np.dot(V2,V3)
2434        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2435        return Tors
2436           
2437    SyOps = []
2438    names = []
2439    for i,atom in enumerate(Oatoms):
2440        names += atom[-1]
2441        Op,unit = Atoms[i][-1]
2442        inv = Op/abs(Op)
2443        m,t = SGData['SGOps'][abs(Op)%100-1]
2444        c = SGData['SGCen'][abs(Op)/100]
2445        SyOps.append([inv,m,t,c,unit])
2446    M = len(Oatoms)
2447    if M == 2:
2448        Val = calcDist(Oatoms,SyOps,Amat)
2449    elif M == 3:
2450        Val = calcAngle(Oatoms,SyOps,Amat)
2451    else:
2452        Val = calcTorsion(Oatoms,SyOps,Amat)
2453   
2454    sigVals = [-0.001,-0.01,-0.01]
2455    sig = sigVals[M-3]
2456    if 'covMatrix' in covData:
2457        dx = .00001
2458        N = M*3
2459        dadx = np.zeros(N)
2460        for i in range(N):
2461            ia = i/3
2462            ix = i%3
2463            Oatoms[ia][ix+1] += dx
2464            if M == 2:
2465                a0 = calcDist(Oatoms,SyOps,Amat)
2466            elif M == 3:
2467                a0 = calcAngle(Oatoms,SyOps,Amat)
2468            else:
2469                a0 = calcTorsion(Oatoms,SyOps,Amat)
2470            Oatoms[ia][ix+1] -= 2*dx
2471            if M == 2:
2472                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2473            elif M == 3:
2474                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2475            else:
2476                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2477        covMatrix = covData['covMatrix']
2478        varyList = covData['varyList']
2479        Vcov = getVCov(names,varyList,covMatrix)
2480        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
2481        if sig < sigVals[M-3]:
2482            sig = sigVals[M-3]
2483   
2484    return Val,sig
2485       
2486def ValEsd(value,esd=0,nTZ=False):
2487    '''Format a floating point number with a given level of precision or
2488    with in crystallographic format with a "esd", as value(esd). If esd is
2489    negative the number is formatted with the level of significant figures
2490    appropriate if abs(esd) were the esd, but the esd is not included.
2491    if the esd is zero, approximately 6 significant figures are printed.
2492    nTZ=True causes "extra" zeros to be removed after the decimal place.
2493    for example:
2494
2495      * "1.235(3)" for value=1.2346 & esd=0.003
2496      * "1.235(3)e4" for value=12346. & esd=30
2497      * "1.235(3)e6" for value=0.12346e7 & esd=3000
2498      * "1.235" for value=1.2346 & esd=-0.003
2499      * "1.240" for value=1.2395 & esd=-0.003
2500      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
2501      * "1.23460" for value=1.2346 & esd=0.0
2502
2503    :param float value: number to be formatted
2504    :param float esd: uncertainty or if esd < 0, specifies level of
2505      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
2506    :param bool nTZ: True to remove trailing zeros (default is False)
2507    :returns: value(esd) or value as a string
2508
2509    '''
2510    # Note: this routine is Python 3 compatible -- I think
2511    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
2512    if math.isnan(value): # invalid value, bail out
2513        return '?'
2514    if math.isnan(esd): # invalid esd, treat as zero
2515        esd = 0
2516        esdoff = 5
2517#    if esd < 1.e-5:
2518#        esd = 0
2519#        esdoff = 5
2520    elif esd != 0:
2521        # transform the esd to a one or two digit integer
2522        l = math.log10(abs(esd)) % 1.
2523        if l < math.log10(cutoff): l+= 1.       
2524        intesd = int(round(10**l)) # esd as integer
2525        # determine the number of digits offset for the esd
2526        esdoff = int(round(math.log10(intesd*1./abs(esd))))
2527    else:
2528        esdoff = 5
2529    valoff = 0
2530    if abs(value) < abs(esdoff): # value is effectively zero
2531        pass
2532    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
2533        # where the digit offset is to the left of the decimal place or where too many
2534        # digits are needed
2535        if abs(value) > 1:
2536            valoff = int(math.log10(abs(value)))
2537        elif abs(value) > 0:
2538            valoff = int(math.log10(abs(value))-0.9999999)
2539        else:
2540            valoff = 0
2541    if esd != 0:
2542        if valoff+esdoff < 0:
2543            valoff = esdoff = 0
2544        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
2545    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
2546        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
2547    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
2548        if abs(value) > 0:           
2549            extra = -math.log10(abs(value))
2550        else:
2551            extra = 0
2552        if extra > 0: extra += 1
2553        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
2554    if esd > 0:
2555        out += ("({:d})").format(intesd)  # add the esd
2556    elif nTZ and '.' in out:
2557        out = out.rstrip('0')  # strip zeros to right of decimal
2558        out = out.rstrip('.')  # and decimal place when not needed
2559    if valoff != 0:
2560        out += ("e{:d}").format(valoff) # add an exponent, when needed
2561    return out
2562   
2563###############################################################################
2564##### Protein validation - "ERRATV2" analysis
2565###############################################################################
2566
2567def validProtein(Phase):
2568   
2569    def sumintact(intact):
2570        return {'CC':intact['CC'],'NN':intact['NN'],'OO':intact['OO'],
2571        'CN':(intact['CN']+intact['NC']),'CO':(intact['CO']+intact['OC']),
2572        'NO':(intact['NO']+intact['ON'])}
2573       
2574    resNames = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
2575        'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
2576# data taken from erratv2.ccp
2577    b1 = np.array([     [0,0,0,0,0,0],
2578          [0,   5040.279078850848200,   3408.805141583649400,   4152.904423767300600,   4236.200004171890200,   5054.781210204625500], 
2579          [0,   3408.805141583648900,   8491.906094010220800,   5958.881777877950300,   1521.387352718486200,   4304.078200827221700], 
2580          [0,   4152.904423767301500,   5958.881777877952100,   7637.167089335050100,   6620.715738223072500,   5287.691183798410700], 
2581          [0,   4236.200004171890200,   1521.387352718486200,   6620.715738223072500,   18368.343774298410000,  4050.797811118806700], 
2582          [0,   5054.781210204625500,   4304.078200827220800,   5287.691183798409800,   4050.797811118806700,   6666.856740479164700]])
2583    avg = np.array([0.,0.192765509919262, 0.195575208778518, 0.275322406824210, 0.059102357035642, 0.233154192767480])
2584    General = Phase['General']
2585    Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
2586    cx,ct,cs,cia = General['AtomPtrs']
2587    Atoms = Phase['Atoms']
2588    cartAtoms = []
2589    xyzmin = 999.*np.ones(3)
2590    xyzmax = -999.*np.ones(3)
2591    #select residue atoms, skipping main chain C & N; make cartesian
2592    for atom in Atoms:
2593        if atom[1] in resNames:
2594            if atom[4].strip() in ['S','Se']:
2595                atom[4] = 'O'
2596            if atom[3].strip() in ['C','N']:    #skip main chain C & N atoms
2597                continue
2598            cartAtoms.append(atom[:cx+3])
2599            cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
2600    XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
2601    xyzmin = np.array([np.min(XYZ.T[i]) for i in [0,1,2]])
2602    xyzmax = np.array([np.max(XYZ.T[i]) for i in [0,1,2]])
2603    nbox = list(np.array(np.ceil((xyzmax-xyzmin)/4.),dtype=int))+[15,]
2604    Boxes = np.zeros(nbox,dtype=int)
2605    iBox = np.array([np.trunc((XYZ.T[i]-xyzmin[i])/4.) for i in [0,1,2]],dtype=int).T
2606    for ib,box in enumerate(iBox):  #put in a try for too many atoms in box (IndexError)?
2607        Boxes[box[0],box[1],box[2],0] += 1
2608        Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
2609    indices = (-1,0,1)
2610    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) 
2611    dsmax = 3.75**2
2612    chains = []
2613    resIntAct = []
2614    chainIntAct = []
2615    res = []
2616    resNames = []
2617    resname = []
2618    newChain = True
2619    intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2620    for ia,atom in enumerate(cartAtoms):
2621        if atom[2] not in chains:   #get chain id & save residue sequence from last chain
2622            chains.append(atom[2])
2623            if len(resIntAct):
2624                resIntAct.append(sumintact(intact))
2625                chainIntAct.append(resIntAct)
2626                resNames += resname
2627                res = []
2628                resname = []
2629                resIntAct = []
2630                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2631                newChain = True
2632        if atom[0] not in res:  #new residue, get residue no.
2633            res.append(atom[0])
2634            resname.append('%s-%s%s'%(atom[2],atom[0],atom[1]))
2635            if not newChain:
2636                resIntAct.append(sumintact(intact))
2637            intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2638            newChain = False
2639        ibox = iBox[ia]         #box location of atom
2640        tgts = []
2641        for unit in Units:      #assemble list of all possible target atoms
2642            jbox = ibox+unit
2643            if np.all(jbox>=0) and np.all((jbox-nbox[:3])<0):               
2644                tgts += list(Boxes[jbox[0],jbox[1],jbox[2]])
2645        tgts = list(set(tgts))
2646        tgts = [tgt for tgt in tgts if atom[1:3] != cartAtoms[tgt][1:3]]    #exclude same residue
2647        tgts = [tgt for tgt in tgts if cartAtoms[tgt][3].strip() not in ['C','N']]  #exclude main chain
2648        tgts = [tgt for tgt in tgts if np.sum((XYZ[ia]-XYZ[tgt])**2) <= dsmax]
2649        for tgt in tgts:
2650            dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
2651            mult = 1.0
2652            if dsqt > 3.25:
2653                mult = 2.*(3.75-dsqt)
2654            intype = atom[4].strip()+cartAtoms[tgt][4].strip()
2655            intact[intype] += mult
2656    resNames += resname
2657    resIntAct.append(sumintact(intact))
2658    chainIntAct.append(resIntAct)
2659    chainProb = []
2660    for ich,chn in enumerate(chains):
2661        IntAct = chainIntAct[ich]
2662        nRes = len(IntAct)
2663        Probs = [0.,0.,0.,0.]
2664        for i in range(4,nRes-4):
2665            mtrx = np.zeros(6)
2666            summ = 0.
2667            for j in range(i-4,i+4):
2668                summ += np.sum(np.array(IntAct[j].values()))
2669                print 'CC',IntAct[j]['CC'],'CN',IntAct[j]['CN'],'CO',IntAct[j]['CO']
2670                print 'NN',IntAct[j]['NN'],'NO',IntAct[j]['NO'],'OO',IntAct[j]['OO']
2671                mtrx[1] += IntAct[j]['CC']
2672                mtrx[2] += IntAct[j]['CN']
2673                mtrx[3] += IntAct[j]['CO']
2674                mtrx[4] += IntAct[j]['NN']
2675                mtrx[5] += IntAct[j]['NO']
2676            mtrx /= summ
2677            mtrx -= avg
2678            prob = np.inner(np.inner(mtrx,b1),mtrx)/36.
2679            print i, mtrx
2680            Probs.append(prob)
2681        Probs += 4*[0.,]
2682        chainProb += Probs
2683    return resNames,chainProb
2684   
2685################################################################################
2686##### Texture fitting stuff
2687################################################################################
2688
2689def FitTexture(General,Gangls,refData,keyList,pgbar):
2690    import pytexture as ptx
2691    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2692   
2693    def printSpHarm(textureData,SHtextureSig):
2694        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
2695        names = ['omega','chi','phi']
2696        namstr = '  names :'
2697        ptstr =  '  values:'
2698        sigstr = '  esds  :'
2699        for name in names:
2700            namstr += '%12s'%('Sample '+name)
2701            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2702            if 'Sample '+name in SHtextureSig:
2703                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2704            else:
2705                sigstr += 12*' '
2706        print namstr
2707        print ptstr
2708        print sigstr
2709        print '\n Texture coefficients:'
2710        SHcoeff = textureData['SH Coeff'][1]
2711        SHkeys = SHcoeff.keys()
2712        nCoeff = len(SHcoeff)
2713        nBlock = nCoeff/10+1
2714        iBeg = 0
2715        iFin = min(iBeg+10,nCoeff)
2716        for block in range(nBlock):
2717            namstr = '  names :'
2718            ptstr =  '  values:'
2719            sigstr = '  esds  :'
2720            for name in SHkeys[iBeg:iFin]:
2721                if 'C' in name:
2722                    namstr += '%12s'%(name)
2723                    ptstr += '%12.3f'%(SHcoeff[name])
2724                    if name in SHtextureSig:
2725                        sigstr += '%12.3f'%(SHtextureSig[name])
2726                    else:
2727                        sigstr += 12*' '
2728            print namstr
2729            print ptstr
2730            print sigstr
2731            iBeg += 10
2732            iFin = min(iBeg+10,nCoeff)
2733           
2734    def Dict2Values(parmdict, varylist):
2735        '''Use before call to leastsq to setup list of values for the parameters
2736        in parmdict, as selected by key in varylist'''
2737        return [parmdict[key] for key in varylist] 
2738       
2739    def Values2Dict(parmdict, varylist, values):
2740        ''' Use after call to leastsq to update the parameter dictionary with
2741        values corresponding to keys in varylist'''
2742        parmdict.update(zip(varylist,values))
2743       
2744    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2745        parmDict.update(zip(varyList,values))
2746        Mat = np.empty(0)
2747        sumObs = 0
2748        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2749        for hist in Gangls.keys():
2750            Refs = refData[hist]
2751            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2752            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2753#            wt = 1./np.max(Refs[:,4],.25)
2754            sumObs += np.sum(wt*Refs[:,5])
2755            Refs[:,6] = 1.
2756            H = Refs[:,:3]
2757            phi,beta = G2lat.CrsAng(H,cell,SGData)
2758            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2759            for item in parmDict:
2760                if 'C' in item:
2761                    L,M,N = eval(item.strip('C'))
2762                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2763                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2764                    Lnorm = G2lat.Lnorm(L)
2765                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2766            mat = wt*(Refs[:,5]-Refs[:,6])
2767            Mat = np.concatenate((Mat,mat))
2768        sumD = np.sum(np.abs(Mat))
2769        R = min(100.,100.*sumD/sumObs)
2770        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2771        print ' Residual: %.3f%%'%(R)
2772        return Mat
2773       
2774    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2775        Mat = np.empty(0)
2776        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2777        for hist in Gangls.keys():
2778            mat = np.zeros((len(varyList),len(refData[hist])))
2779            Refs = refData[hist]
2780            H = Refs[:,:3]
2781            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2782#            wt = 1./np.max(Refs[:,4],.25)
2783            phi,beta = G2lat.CrsAng(H,cell,SGData)
2784            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2785            for j,item in enumerate(varyList):
2786                if 'C' in item:
2787                    L,M,N = eval(item.strip('C'))
2788                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2789                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2790                    Lnorm = G2lat.Lnorm(L)
2791                    mat[j] = -wt*Lnorm*Kcl*Ksl
2792                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2793                        try:
2794                            l = varyList.index(itema)
2795                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2796                        except ValueError:
2797                            pass
2798            if len(Mat):
2799                Mat = np.concatenate((Mat,mat.T))
2800            else:
2801                Mat = mat.T
2802        print 'deriv'
2803        return Mat
2804
2805    print ' Fit texture for '+General['Name']
2806    SGData = General['SGData']
2807    cell = General['Cell'][1:7]
2808    Texture = General['SH Texture']
2809    if not Texture['Order']:
2810        return 'No spherical harmonics coefficients'
2811    varyList = []
2812    parmDict = copy.copy(Texture['SH Coeff'][1])
2813    for item in ['Sample omega','Sample chi','Sample phi']:
2814        parmDict[item] = Texture[item][1]
2815        if Texture[item][0]:
2816            varyList.append(item)
2817    if Texture['SH Coeff'][0]:
2818        varyList += Texture['SH Coeff'][1].keys()
2819    while True:
2820        begin = time.time()
2821        values =  np.array(Dict2Values(parmDict, varyList))
2822        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
2823            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
2824        ncyc = int(result[2]['nfev']/2)
2825        if ncyc:
2826            runtime = time.time()-begin   
2827            chisq = np.sum(result[2]['fvec']**2)
2828            Values2Dict(parmDict, varyList, result[0])
2829            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
2830            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
2831            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2832            try:
2833                sig = np.sqrt(np.diag(result[1])*GOF)
2834                if np.any(np.isnan(sig)):
2835                    print '*** Least squares aborted - some invalid esds possible ***'
2836                break                   #refinement succeeded - finish up!
2837            except ValueError:          #result[1] is None on singular matrix
2838                print '**** Refinement failed - singular matrix ****'
2839                return None
2840        else:
2841            break
2842   
2843    if ncyc:
2844        for parm in parmDict:
2845            if 'C' in parm:
2846                Texture['SH Coeff'][1][parm] = parmDict[parm]
2847            else:
2848                Texture[parm][1] = parmDict[parm] 
2849        sigDict = dict(zip(varyList,sig))
2850        printSpHarm(Texture,sigDict)
2851       
2852    return None
2853   
2854################################################################################
2855##### Fourier & charge flip stuff
2856################################################################################
2857
2858def adjHKLmax(SGData,Hmax):
2859    '''default doc string
2860   
2861    :param type name: description
2862   
2863    :returns: type name: description
2864   
2865    '''
2866    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2867        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
2868        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
2869        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2870    else:
2871        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
2872        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
2873        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2874
2875def OmitMap(data,reflDict,pgbar=None):
2876    '''default doc string
2877   
2878    :param type name: description
2879   
2880    :returns: type name: description
2881   
2882    '''
2883    generalData = data['General']
2884    if not generalData['Map']['MapType']:
2885        print '**** ERROR - Fourier map not defined'
2886        return
2887    mapData = generalData['Map']
2888    dmin = mapData['Resolution']
2889    SGData = generalData['SGData']
2890    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2891    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2892    cell = generalData['Cell'][1:8]       
2893    A = G2lat.cell2A(cell[:6])
2894    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2895    adjHKLmax(SGData,Hmax)
2896    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2897    time0 = time.time()
2898    for iref,ref in enumerate(reflDict['RefList']):
2899        if ref[4] >= dmin:
2900            Fosq,Fcsq,ph = ref[8:11]
2901            Uniq = np.inner(ref[:3],SGMT)
2902            Phi = np.inner(ref[:3],SGT)
2903            for i,hkl in enumerate(Uniq):        #uses uniq
2904                hkl = np.asarray(hkl,dtype='i')
2905                dp = 360.*Phi[i]                #and phi
2906                a = cosd(ph+dp)
2907                b = sind(ph+dp)
2908                phasep = complex(a,b)
2909                phasem = complex(a,-b)
2910                if '2Fo-Fc' in mapData['MapType']:
2911                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
2912                else:
2913                    F = np.sqrt(Fosq)
2914                h,k,l = hkl+Hmax
2915                Fhkl[h,k,l] = F*phasep
2916                h,k,l = -hkl+Hmax
2917                Fhkl[h,k,l] = F*phasem
2918    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2919    M = np.mgrid[0:4,0:4,0:4]
2920    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
2921    iBeg = blkIds*rho0.shape/4
2922    iFin = (blkIds+1)*rho0.shape/4
2923    rho_omit = np.zeros_like(rho0)
2924    nBlk = 0
2925    for iB,iF in zip(iBeg,iFin):
2926        rho1 = np.copy(rho0)
2927        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
2928        Fnew = fft.ifftshift(fft.ifftn(rho1))
2929        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
2930        phase = Fnew/np.absolute(Fnew)
2931        OFhkl = np.absolute(Fhkl)*phase
2932        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
2933        rho_omit[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = np.copy(rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]])
2934        nBlk += 1
2935        pgbar.Update(nBlk)
2936    mapData['rho'] = np.real(rho_omit)/cell[6]
2937    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2938    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2939    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2940    return mapData
2941   
2942def FourierMap(data,reflDict):
2943    '''default doc string
2944   
2945    :param type name: description
2946   
2947    :returns: type name: description
2948   
2949    '''
2950    generalData = data['General']
2951    mapData = generalData['Map']
2952    dmin = mapData['Resolution']
2953    SGData = generalData['SGData']
2954    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2955    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2956    cell = generalData['Cell'][1:8]       
2957    A = G2lat.cell2A(cell[:6])
2958    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2959    adjHKLmax(SGData,Hmax)
2960    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2961#    Fhkl[0,0,0] = generalData['F000X']
2962    time0 = time.time()
2963    for iref,ref in enumerate(reflDict['RefList']):
2964        if ref[4] > dmin:
2965            Fosq,Fcsq,ph = ref[8:11]
2966            Uniq = np.inner(ref[:3],SGMT)
2967            Phi = np.inner(ref[:3],SGT)
2968            for i,hkl in enumerate(Uniq):        #uses uniq
2969                hkl = np.asarray(hkl,dtype='i')
2970                dp = 360.*Phi[i]                #and phi
2971                a = cosd(ph+dp)
2972                b = sind(ph+dp)
2973                phasep = complex(a,b)
2974                phasem = complex(a,-b)
2975                if 'Fobs' in mapData['MapType']:
2976                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
2977                    h,k,l = hkl+Hmax
2978                    Fhkl[h,k,l] = F*phasep
2979                    h,k,l = -hkl+Hmax
2980                    Fhkl[h,k,l] = F*phasem
2981                elif 'Fcalc' in mapData['MapType']:
2982                    F = np.sqrt(Fcsq)
2983                    h,k,l = hkl+Hmax
2984                    Fhkl[h,k,l] = F*phasep
2985                    h,k,l = -hkl+Hmax
2986                    Fhkl[h,k,l] = F*phasem
2987                elif 'delt-F' in mapData['MapType']:
2988                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
2989                    h,k,l = hkl+Hmax
2990                    Fhkl[h,k,l] = dF*phasep
2991                    h,k,l = -hkl+Hmax
2992                    Fhkl[h,k,l] = dF*phasem
2993                elif '2*Fo-Fc' in mapData['MapType']:
2994                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
2995                    h,k,l = hkl+Hmax
2996                    Fhkl[h,k,l] = F*phasep
2997                    h,k,l = -hkl+Hmax
2998                    Fhkl[h,k,l] = F*phasem
2999                elif 'Patterson' in mapData['MapType']:
3000                    h,k,l = hkl+Hmax
3001                    Fhkl[h,k,l] = complex(Fosq,0.)
3002                    h,k,l = -hkl+Hmax
3003                    Fhkl[h,k,l] = complex(Fosq,0.)
3004    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3005    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
3006    mapData['Type'] = reflDict['Type']
3007    mapData['rho'] = np.real(rho)
3008    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3009    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3010   
3011def Fourier4DMap(data,reflDict):
3012    '''default doc string
3013   
3014    :param type name: description
3015   
3016    :returns: type name: description
3017   
3018    '''
3019    generalData = data['General']
3020    map4DData = generalData['4DmapData']
3021    mapData = generalData['Map']
3022    dmin = mapData['Resolution']
3023    SGData = generalData['SGData']
3024    SSGData = generalData['SSGData']
3025    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3026    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3027    cell = generalData['Cell'][1:8]       
3028    A = G2lat.cell2A(cell[:6])
3029    maxM = 4
3030    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
3031    adjHKLmax(SGData,Hmax)
3032    Hmax = np.asarray(Hmax,dtype='i')+1
3033    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3034    time0 = time.time()
3035    for iref,ref in enumerate(reflDict['RefList']):
3036        if ref[5] > dmin:
3037            Fosq,Fcsq,ph = ref[9:12]
3038            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
3039            Uniq = np.inner(ref[:4],SSGMT)
3040            Phi = np.inner(ref[:4],SSGT)
3041            for i,hkl in enumerate(Uniq):        #uses uniq
3042                hkl = np.asarray(hkl,dtype='i')
3043                dp = 360.*Phi[i]                #and phi
3044                a = cosd(ph+dp)
3045                b = sind(ph+dp)
3046                phasep = complex(a,b)
3047                phasem = complex(a,-b)
3048                if 'Fobs' in mapData['MapType']:
3049                    F = np.sqrt(Fosq)
3050                    h,k,l,m = hkl+Hmax
3051                    Fhkl[h,k,l,m] = F*phasep
3052                    h,k,l,m = -hkl+Hmax
3053                    Fhkl[h,k,l,m] = F*phasem
3054                elif 'Fcalc' in mapData['MapType']:
3055                    F = np.sqrt(Fcsq)
3056                    h,k,l,m = hkl+Hmax
3057                    Fhkl[h,k,l,m] = F*phasep
3058                    h,k,l,m = -hkl+Hmax
3059                    Fhkl[h,k,l,m] = F*phasem                   
3060                elif 'delt-F' in mapData['MapType']:
3061                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
3062                    h,k,l,m = hkl+Hmax
3063                    Fhkl[h,k,l,m] = dF*phasep
3064                    h,k,l,m = -hkl+Hmax
3065                    Fhkl[h,k,l,m] = dF*phasem
3066    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
3067    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
3068    map4DData['rho'] = np.real(SSrho)
3069    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3070    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3071    map4DData['Type'] = reflDict['Type']
3072    mapData['Type'] = reflDict['Type']
3073    mapData['rho'] = np.real(rho)
3074    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3075    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3076    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
3077
3078# map printing for testing purposes
3079def printRho(SGLaue,rho,rhoMax):                         
3080    '''default doc string
3081   
3082    :param type name: description
3083   
3084    :returns: type name: description
3085   
3086    '''
3087    dim = len(rho.shape)
3088    if dim == 2:
3089        ix,jy = rho.shape
3090        for j in range(jy):
3091            line = ''
3092            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3093                line += (jy-j)*'  '
3094            for i in range(ix):
3095                r = int(100*rho[i,j]/rhoMax)
3096                line += '%4d'%(r)
3097            print line+'\n'
3098    else:
3099        ix,jy,kz = rho.shape
3100        for k in range(kz):
3101            print 'k = ',k
3102            for j in range(jy):
3103                line = ''
3104                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3105                    line += (jy-j)*'  '
3106                for i in range(ix):
3107                    r = int(100*rho[i,j,k]/rhoMax)
3108                    line += '%4d'%(r)
3109                print line+'\n'
3110## keep this
3111               
3112def findOffset(SGData,A,Fhkl):   
3113    '''default doc string
3114   
3115    :param type name: description
3116   
3117    :returns: type name: description
3118   
3119    '''
3120    if SGData['SpGrp'] == 'P 1':
3121        return [0,0,0]   
3122    hklShape = Fhkl.shape
3123    hklHalf = np.array(hklShape)/2
3124    sortHKL = np.argsort(Fhkl.flatten())
3125    Fdict = {}
3126    for hkl in sortHKL:
3127        HKL = np.unravel_index(hkl,hklShape)
3128        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
3129        if F == 0.:
3130            break
3131        Fdict['%.6f'%(np.absolute(F))] = hkl
3132    Flist = np.flipud(np.sort(Fdict.keys()))
3133    F = str(1.e6)
3134    i = 0
3135    DH = []
3136    Dphi = []
3137    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3138    for F in Flist:
3139        hkl = np.unravel_index(Fdict[F],hklShape)
3140        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
3141            continue
3142        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
3143        Uniq = np.array(Uniq,dtype='i')
3144        Phi = np.array(Phi)
3145        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
3146        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3147        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
3148        ang0 = np.angle(Fh0,deg=True)/360.
3149        for H,phi in zip(Uniq,Phi)[1:]:
3150            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
3151            dH = H-hkl
3152            dang = ang-ang0
3153            DH.append(dH)
3154            Dphi.append((dang+.5) % 1.0)
3155        if i > 20 or len(DH) > 30:
3156            break
3157        i += 1
3158    DH = np.array(DH)
3159    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
3160    Dphi = np.array(Dphi)
3161    steps = np.array(hklShape)
3162    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
3163    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
3164    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
3165    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3166    hist,bins = np.histogram(Mmap,bins=1000)
3167#    for i,item in enumerate(hist[:10]):
3168#        print item,bins[i]
3169    chisq = np.min(Mmap)
3170    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3171    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
3172#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3173    return DX
3174   
3175def ChargeFlip(data,reflDict,pgbar):
3176    '''default doc string
3177   
3178    :param type name: description
3179   
3180    :returns: type name: description
3181   
3182    '''
3183    generalData = data['General']
3184    mapData = generalData['Map']
3185    flipData = generalData['Flip']
3186    FFtable = {}
3187    if 'None' not in flipData['Norm element']:
3188        normElem = flipData['Norm element'].upper()
3189        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3190        for ff in FFs:
3191            if ff['Symbol'] == normElem:
3192                FFtable.update(ff)
3193    dmin = flipData['Resolution']
3194    SGData = generalData['SGData']
3195    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3196    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3197    cell = generalData['Cell'][1:8]       
3198    A = G2lat.cell2A(cell[:6])
3199    Vol = cell[6]
3200    im = 0
3201    if generalData['Type'] in ['modulated','magnetic',]:
3202        im = 1
3203    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3204    adjHKLmax(SGData,Hmax)
3205    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3206    time0 = time.time()
3207    for iref,ref in enumerate(reflDict['RefList']):
3208        dsp = ref[4+im]
3209        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
3210            continue
3211        if dsp > dmin:
3212            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3213            if FFtable:
3214                SQ = 0.25/dsp**2
3215                ff *= G2el.ScatFac(FFtable,SQ)[0]
3216            if ref[8+im] > 0.:         #use only +ve Fobs**2
3217                E = np.sqrt(ref[8+im])/ff
3218            else:
3219                E = 0.
3220            ph = ref[10]
3221            ph = rn.uniform(0.,360.)
3222            Uniq = np.inner(ref[:3],SGMT)
3223            Phi = np.inner(ref[:3],SGT)
3224            for i,hkl in enumerate(Uniq):        #uses uniq
3225                hkl = np.asarray(hkl,dtype='i')
3226                dp = 360.*Phi[i]                #and phi
3227                a = cosd(ph+dp)
3228                b = sind(ph+dp)
3229                phasep = complex(a,b)
3230                phasem = complex(a,-b)
3231                h,k,l = hkl+Hmax
3232                Ehkl[h,k,l] = E*phasep
3233                h,k,l = -hkl+Hmax
3234                Ehkl[h,k,l] = E*phasem
3235#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3236    testHKL = np.array(flipData['testHKL'])+Hmax
3237    CEhkl = copy.copy(Ehkl)
3238    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3239    Emask = ma.getmask(MEhkl)
3240    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3241    Ncyc = 0
3242    old = np.seterr(all='raise')
3243    twophases = []
3244    while True:       
3245        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3246        CEsig = np.std(CErho)
3247        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3248        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3249        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3250        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3251        phase = CFhkl/np.absolute(CFhkl)
3252        twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
3253        CEhkl = np.absolute(Ehkl)*phase
3254        Ncyc += 1
3255        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3256        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3257        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3258        if Rcf < 5.:
3259            break
3260        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3261        if not GoOn or Ncyc > 10000:
3262            break
3263    np.seterr(**old)
3264    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3265    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
3266    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3267    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3268       
3269    mapData['Rcf'] = Rcf
3270    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3271    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3272    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3273    mapData['Type'] = reflDict['Type']
3274    return mapData,twophases
3275   
3276def findSSOffset(SGData,SSGData,A,Fhklm):   
3277    '''default doc string
3278   
3279    :param type name: description
3280   
3281    :returns: type name: description
3282   
3283    '''
3284    if SGData['SpGrp'] == 'P 1':
3285        return [0,0,0,0]   
3286    hklmShape = Fhklm.shape
3287    hklmHalf = np.array(hklmShape)/2
3288    sortHKLM = np.argsort(Fhklm.flatten())
3289    Fdict = {}
3290    for hklm in sortHKLM:
3291        HKLM = np.unravel_index(hklm,hklmShape)
3292        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
3293        if F == 0.:
3294            break
3295        Fdict['%.6f'%(np.absolute(F))] = hklm
3296    Flist = np.flipud(np.sort(Fdict.keys()))
3297    F = str(1.e6)
3298    i = 0
3299    DH = []
3300    Dphi = []
3301    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3302    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3303    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3304    for F in Flist:
3305        hklm = np.unravel_index(Fdict[F],hklmShape)
3306        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
3307            continue
3308        Uniq = np.inner(hklm-hklmHalf,SSGMT)
3309        Phi = np.inner(hklm-hklmHalf,SSGT)
3310        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
3311        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3312        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
3313        ang0 = np.angle(Fh0,deg=True)/360.
3314        for H,phi in zip(Uniq,Phi)[1:]:
3315            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
3316            dH = H-hklm
3317            dang = ang-ang0
3318            DH.append(dH)
3319            Dphi.append((dang+.5) % 1.0)
3320        if i > 20 or len(DH) > 30:
3321            break
3322        i += 1
3323    DH = np.array(DH)
3324    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
3325    Dphi = np.array(Dphi)
3326    steps = np.array(hklmShape)
3327    X,Y,Z,T = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2],0:1:1./steps[3]]
3328    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
3329    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
3330    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3331    hist,bins = np.histogram(Mmap,bins=1000)
3332#    for i,item in enumerate(hist[:10]):
3333#        print item,bins[i]
3334    chisq = np.min(Mmap)
3335    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3336    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
3337#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3338    return DX
3339   
3340def SSChargeFlip(data,reflDict,pgbar):
3341    '''default doc string
3342   
3343    :param type name: description
3344   
3345    :returns: type name: description
3346   
3347    '''
3348    generalData = data['General']
3349    mapData = generalData['Map']
3350    map4DData = {}
3351    flipData = generalData['Flip']
3352    FFtable = {}
3353    if 'None' not in flipData['Norm element']:
3354        normElem = flipData['Norm element'].upper()
3355        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3356        for ff in FFs:
3357            if ff['Symbol'] == normElem:
3358                FFtable.update(ff)
3359    dmin = flipData['Resolution']
3360    SGData = generalData['SGData']
3361    SSGData = generalData['SSGData']
3362    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3363    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3364    cell = generalData['Cell'][1:8]       
3365    A = G2lat.cell2A(cell[:6])
3366    Vol = cell[6]
3367    maxM = 4
3368    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3369    adjHKLmax(SGData,Hmax)
3370    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3371    time0 = time.time()
3372    for iref,ref in enumerate(reflDict['RefList']):
3373        dsp = ref[5]
3374        if dsp > dmin:
3375            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3376            if FFtable:
3377                SQ = 0.25/dsp**2
3378                ff *= G2el.ScatFac(FFtable,SQ)[0]
3379            if ref[9] > 0.:         #use only +ve Fobs**2
3380                E = np.sqrt(ref[9])/ff
3381            else:
3382                E = 0.
3383            ph = ref[11]
3384            ph = rn.uniform(0.,360.)
3385            Uniq = np.inner(ref[:4],SSGMT)
3386            Phi = np.inner(ref[:4],SSGT)
3387            for i,hklm in enumerate(Uniq):        #uses uniq
3388                hklm = np.asarray(hklm,dtype='i')
3389                dp = 360.*Phi[i]                #and phi
3390                a = cosd(ph+dp)
3391                b = sind(ph+dp)
3392                phasep = complex(a,b)
3393                phasem = complex(a,-b)
3394                h,k,l,m = hklm+Hmax
3395                Ehkl[h,k,l,m] = E*phasep
3396                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
3397                Ehkl[h,k,l,m] = E*phasem
3398#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3399    CEhkl = copy.copy(Ehkl)
3400    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3401    Emask = ma.getmask(MEhkl)
3402    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3403    Ncyc = 0
3404    old = np.seterr(all='raise')
3405    while True:       
3406        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3407        CEsig = np.std(CErho)
3408        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3409        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3410        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3411        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3412        phase = CFhkl/np.absolute(CFhkl)
3413        CEhkl = np.absolute(Ehkl)*phase
3414        Ncyc += 1
3415        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3416        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3417        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3418        if Rcf < 5.:
3419            break
3420        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3421        if not GoOn or Ncyc > 10000:
3422            break
3423    np.seterr(**old)
3424    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3425    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
3426    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
3427    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3428    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3429
3430    mapData['Rcf'] = Rcf
3431    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3432    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3433    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3434    mapData['Type'] = reflDict['Type']
3435
3436    map4DData['Rcf'] = Rcf
3437    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2),roll[3],axis=3))
3438    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3439    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3440    map4DData['Type'] = reflDict['Type']
3441    return mapData,map4DData
3442   
3443def getRho(xyz,mapData):
3444    ''' get scattering density at a point by 8-point interpolation
3445    param xyz: coordinate to be probed
3446    param: mapData: dict of map data
3447   
3448    :returns: density at xyz
3449    '''
3450    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3451    if not len(mapData):
3452        return 0.0
3453    rho = copy.copy(mapData['rho'])     #don't mess up original
3454    if not len(rho):
3455        return 0.0
3456    mapShape = np.array(rho.shape)
3457    mapStep = 1./mapShape
3458    X = np.array(xyz)%1.    #get into unit cell
3459    I = np.array(X*mapShape,dtype='int')
3460    D = X-I*mapStep         #position inside map cell
3461    D12 = D[0]*D[1]
3462    D13 = D[0]*D[2]
3463    D23 = D[1]*D[2]
3464    D123 = np.prod(D)
3465    Rho = rollMap(rho,-I)    #shifts map so point is in corner
3466    R = Rho[0,0,0]*(1.-np.sum(D))+Rho[1,0,0]*D[0]+Rho[0,1,0]*D[1]+Rho[0,0,1]*D[2]+  \
3467        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
3468        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
3469        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
3470    return R
3471       
3472def SearchMap(generalData,drawingData,Neg=False):
3473    '''Does a search of a density map for peaks meeting the criterion of peak
3474    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3475    mapData is data['General']['mapData']; the map is also in mapData.
3476
3477    :param generalData: the phase data structure; includes the map
3478    :param drawingData: the drawing data structure
3479    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
3480
3481    :returns: (peaks,mags,dzeros) where
3482
3483        * peaks : ndarray
3484          x,y,z positions of the peaks found in the map
3485        * mags : ndarray
3486          the magnitudes of the peaks
3487        * dzeros : ndarray
3488          the distance of the peaks from  the unit cell origin
3489        * dcent : ndarray
3490          the distance of the peaks from  the unit cell center
3491
3492    '''       
3493    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3494   
3495    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
3496   
3497#    def noDuplicate(xyz,peaks,Amat):
3498#        XYZ = np.inner(Amat,xyz)
3499#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3500#            print ' Peak',xyz,' <0.5A from another peak'
3501#            return False
3502#        return True
3503#                           
3504    def fixSpecialPos(xyz,SGData,Amat):
3505        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3506        X = []
3507        xyzs = [equiv[0] for equiv in equivs]
3508        for x in xyzs:
3509            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3510                X.append(x)
3511        if len(X) > 1:
3512            return np.average(X,axis=0)
3513        else:
3514            return xyz
3515       
3516    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3517        Mag,x0,y0,z0,sig = parms
3518        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
3519#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
3520        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
3521       
3522    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3523        Mag,x0,y0,z0,sig = parms
3524        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3525        return M
3526       
3527    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3528        Mag,x0,y0,z0,sig = parms
3529        dMdv = np.zeros(([5,]+list(rX.shape)))
3530        delt = .01
3531        for i in range(5):
3532            parms[i] -= delt
3533            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3534            parms[i] += 2.*delt
3535            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3536            parms[i] -= delt
3537            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3538        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3539        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3540        dMdv = np.reshape(dMdv,(5,rX.size))
3541        Hess = np.inner(dMdv,dMdv)
3542       
3543        return Vec,Hess
3544       
3545    SGData = generalData['SGData']
3546    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3547    peaks = []
3548    mags = []
3549    dzeros = []
3550    dcent = []
3551    try:
3552        mapData = generalData['Map']
3553        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3554        if Neg:
3555            rho = -copy.copy(mapData['rho'])     #flip +/-
3556        else:
3557            rho = copy.copy(mapData['rho'])     #don't mess up original
3558        mapHalf = np.array(rho.shape)/2
3559        res = mapData['Resolution']
3560        incre = np.array(rho.shape,dtype=np.float)
3561        step = max(1.0,1./res)+1
3562        steps = np.array((3*[step,]),dtype='int32')
3563    except KeyError:
3564        print '**** ERROR - Fourier map not defined'
3565        return peaks,mags
3566    rhoMask = ma.array(rho,mask=(rho<contLevel))
3567    indices = (-1,0,1)
3568    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3569    for roll in rolls:
3570        if np.any(roll):
3571            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3572    indx = np.transpose(rhoMask.nonzero())
3573    peaks = indx/incre
3574    mags = rhoMask[rhoMask.nonzero()]
3575    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3576        rho = rollMap(rho,ind)
3577        rMM = mapHalf-steps
3578        rMP = mapHalf+steps+1
3579        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3580        peakInt = np.sum(rhoPeak)*res**3
3581        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3582        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
3583        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3584            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3585        x1 = result[0]
3586        if not np.any(x1 < 0):
3587            peak = (np.array(x1[1:4])-ind)/incre
3588        peak = fixSpecialPos(peak,SGData,Amat)
3589        rho = rollMap(rho,-ind)
3590    cent = np.ones(3)*.5     
3591    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3592    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3593    if Neg:     #want negative magnitudes for negative peaks
3594        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3595    else:
3596        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3597   
3598def sortArray(data,pos,reverse=False):
3599    '''data is a list of items
3600    sort by pos in list; reverse if True
3601    '''
3602    T = []
3603    for i,M in enumerate(data):
3604        try:
3605            T.append((M[pos],i))
3606        except IndexError:
3607            return data
3608    D = dict(zip(T,data))
3609    T.sort()
3610    if reverse:
3611        T.reverse()
3612    X = []
3613    for key in T:
3614        X.append(D[key])
3615    return X
3616
3617def PeaksEquiv(data,Ind):
3618    '''Find the equivalent map peaks for those selected. Works on the
3619    contents of data['Map Peaks'].
3620
3621    :param data: the phase data structure
3622    :param list Ind: list of selected peak indices
3623    :returns: augmented list of peaks including those related by symmetry to the
3624      ones in Ind
3625
3626    '''       
3627    def Duplicate(xyz,peaks,Amat):
3628        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3629            return True
3630        return False
3631                           
3632    generalData = data['General']
3633    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3634    SGData = generalData['SGData']
3635    mapPeaks = data['Map Peaks']
3636    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3637    Indx = {}
3638    for ind in Ind:
3639        xyz = np.array(mapPeaks[ind][1:4])
3640        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3641        for jnd,xyz in enumerate(XYZ):       
3642            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3643    Ind = []
3644    for ind in Indx:
3645        if Indx[ind]:
3646            Ind.append(ind)
3647    return Ind
3648               
3649def PeaksUnique(data,Ind):
3650    '''Finds the symmetry unique set of peaks from those selected. Works on the
3651    contents of data['Map Peaks'].
3652
3653    :param data: the phase data structure
3654    :param list Ind: list of selected peak indices
3655    :returns: the list of symmetry unique peaks from among those given in Ind
3656
3657    '''       
3658#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3659
3660    def noDuplicate(xyz,peaks,Amat):
3661        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3662            return False
3663        return True
3664                           
3665    generalData = data['General']
3666    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3667    SGData = generalData['SGData']
3668    mapPeaks = data['Map Peaks']
3669    Indx = {}
3670    XYZ = {}
3671    for ind in Ind:
3672        XYZ[ind] = np.array(mapPeaks[ind][1:4])
3673        Indx[ind] = True
3674    for ind in Ind:
3675        if Indx[ind]:
3676            xyz = XYZ[ind]
3677            for jnd in Ind:
3678                if ind != jnd and Indx[jnd]:                       
3679                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3680                    xyzs = np.array([equiv[0] for equiv in Equiv])
3681                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3682    Ind = []
3683    for ind in Indx:
3684        if Indx[ind]:
3685            Ind.append(ind)
3686    return Ind
3687   
3688################################################################################
3689##### single peak fitting profile fxn stuff
3690################################################################################
3691
3692def getCWsig(ins,pos):
3693    '''get CW peak profile sigma^2
3694   
3695    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3696      as values only
3697    :param float pos: 2-theta of peak
3698    :returns: float getCWsig: peak sigma^2
3699   
3700    '''
3701    tp = tand(pos/2.0)
3702    return ins['U']*tp**2+ins['V']*tp+ins['W']
3703   
3704def getCWsigDeriv(pos):
3705    '''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3706   
3707    :param float pos: 2-theta of peak
3708   
3709    :returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3710   
3711    '''
3712    tp = tand(pos/2.0)
3713    return tp**2,tp,1.0
3714   
3715def getCWgam(ins,pos):
3716    '''get CW peak profile gamma
3717   
3718    :param dict ins: instrument parameters with at least 'X' & 'Y'
3719      as values only
3720    :param float pos: 2-theta of peak
3721    :returns: float getCWgam: peak gamma
3722   
3723    '''
3724    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
3725   
3726def getCWgamDeriv(pos):
3727    '''get derivatives of CW peak profile gamma wrt X & Y
3728   
3729    :param float pos: 2-theta of peak
3730   
3731    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3732   
3733    '''
3734    return 1./cosd(pos/2.0),tand(pos/2.0)
3735   
3736def getTOFsig(ins,dsp):
3737    '''get TOF peak profile sigma^2
3738   
3739    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3740      as values only
3741    :param float dsp: d-spacing of peak
3742   
3743    :returns: float getTOFsig: peak sigma^2
3744   
3745    '''
3746    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
3747   
3748def getTOFsigDeriv(dsp):
3749    '''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3750   
3751    :param float dsp: d-spacing of peak
3752   
3753    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3754   
3755    '''
3756    return 1.0,dsp**2,dsp**4,1./dsp**2
3757   
3758def getTOFgamma(ins,dsp):
3759    '''get TOF peak profile gamma
3760   
3761    :param dict ins: instrument parameters with at least 'X' & 'Y'
3762      as values only
3763    :param float dsp: d-spacing of peak
3764   
3765    :returns: float getTOFgamma: peak gamma
3766   
3767    '''
3768    return ins['X']*dsp+ins['Y']*dsp**2
3769   
3770def getTOFgammaDeriv(dsp):
3771    '''get derivatives of TOF peak profile gamma wrt X & Y
3772   
3773    :param float dsp: d-spacing of peak
3774   
3775    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3776   
3777    '''
3778    return dsp,dsp**2
3779   
3780def getTOFbeta(ins,dsp):
3781    '''get TOF peak profile beta
3782   
3783    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3784      as values only
3785    :param float dsp: d-spacing of peak
3786   
3787    :returns: float getTOFbeta: peak beat
3788   
3789    '''
3790    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
3791   
3792def getTOFbetaDeriv(dsp):
3793    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3794   
3795    :param float dsp: d-spacing of peak
3796   
3797    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3798   
3799    '''
3800    return 1.0,1./dsp**4,1./dsp**2
3801   
3802def getTOFalpha(ins,dsp):
3803    '''get TOF peak profile alpha
3804   
3805    :param dict ins: instrument parameters with at least 'alpha'
3806      as values only
3807    :param float dsp: d-spacing of peak
3808   
3809    :returns: flaot getTOFalpha: peak alpha
3810   
3811    '''
3812    return ins['alpha']/dsp
3813   
3814def getTOFalphaDeriv(dsp):
3815    '''get derivatives of TOF peak profile beta wrt alpha
3816   
3817    :param float dsp: d-spacing of peak
3818   
3819    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
3820   
3821    '''
3822    return 1./dsp
3823   
3824def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
3825    '''set starting peak parameters for single peak fits from plot selection or auto selection
3826   
3827    :param dict Parms: instrument parameters dictionary
3828    :param dict Parms2: table lookup for TOF profile coefficients
3829    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
3830    :param float mag: peak top magnitude from pick
3831    :param bool ifQ: True if pos in Q
3832    :param bool useFit: True if use fitted CW Parms values (not defaults)
3833   
3834    :returns: list XY: peak list entry:
3835        for CW: [pos,0,mag,1,sig,0,gam,0]
3836        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3837        NB: mag refinement set by default, all others off
3838   
3839    '''
3840    ind = 0
3841    if useFit:
3842        ind = 1
3843    ins = {}
3844    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
3845        for x in ['U','V','W','X','Y']:
3846            ins[x] = Parms[x][ind]
3847        if ifQ:                              #qplot - convert back to 2-theta
3848            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
3849        sig = getCWsig(ins,pos)
3850        gam = getCWgam(ins,pos)           
3851        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
3852    else:
3853        if ifQ:
3854            dsp = 2.*np.pi/pos
3855            pos = Parms['difC']*dsp
3856        else:
3857            dsp = pos/Parms['difC'][1]
3858        if 'Pdabc' in Parms2:
3859            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
3860                ins[x] = Parms[x][ind]
3861            Pdabc = Parms2['Pdabc'].T
3862            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
3863            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
3864        else:
3865            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
3866                ins[x] = Parms[x][ind]
3867            alp = getTOFalpha(ins,dsp)
3868            bet = getTOFbeta(ins,dsp)
3869        sig = getTOFsig(ins,dsp)
3870        gam = getTOFgamma(ins,dsp)
3871        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3872    return XY
3873   
3874################################################################################
3875##### MC/SA stuff
3876################################################################################
3877
3878#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
3879# Original Author: Travis Oliphant 2002
3880# Bug-fixes in 2006 by Tim Leslie
3881
3882
3883import numpy
3884from numpy import asarray, exp, squeeze, sign, \
3885     all, shape, array, where
3886from numpy import random
3887
3888#__all__ = ['anneal']
3889
3890_double_min = numpy.finfo(float).min
3891_double_max = numpy.finfo(float).max
3892class base_schedule(object):
3893    def __init__(self):
3894        self.dwell = 20
3895        self.lower = 0.
3896        self.upper = 1.
3897        self.Ninit = 50
3898        self.accepted = 0
3899        self.tests = 0
3900        self.feval = 0
3901        self.k = 0
3902        self.T = None
3903
3904    def init(self, **options):
3905        self.__dict__.update(options)
3906        self.lower = asarray(self.lower)
3907        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
3908        self.upper = asarray(self.upper)
3909        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
3910        self.k = 0
3911        self.accepted = 0
3912        self.feval = 0
3913        self.tests = 0
3914
3915    def getstart_temp(self, best_state):
3916        """ Find a matching starting temperature and starting parameters vector
3917        i.e. find x0 such that func(x0) = T0.
3918
3919        :param best_state: _state
3920            A _state object to store the function value and x0 found.
3921
3922        :returns: x0 : array
3923            The starting parameters vector.
3924        """
3925
3926        assert(not self.dims is None)
3927        lrange = self.lower
3928        urange = self.upper
3929        fmax = _double_min
3930        fmin = _double_max
3931        for _ in range(self.Ninit):
3932            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
3933            fval = self.func(x0, *self.args)
3934            self.feval += 1
3935            if fval > fmax:
3936                fmax = fval
3937            if fval < fmin:
3938                fmin = fval
3939                best_state.cost = fval
3940                best_state.x = array(x0)
3941
3942        self.T0 = (fmax-fmin)*1.5
3943        return best_state.x
3944       
3945    def set_range(self,x0,frac):
3946        delrange = frac*(self.upper-self.lower)
3947        self.upper = x0+delrange
3948        self.lower = x0-delrange
3949
3950    def accept_test(self, dE):
3951        T = self.T
3952        self.tests += 1
3953        if dE < 0:
3954            self.accepted += 1
3955            return 1
3956        p = exp(-dE*1.0/T)
3957        if (p > random.uniform(0.0, 1.0)):
3958            self.accepted += 1
3959            return 1
3960        return 0
3961
3962    def update_guess(self, x0):
3963        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3964
3965    def update_temp(self, x0):
3966        pass
3967
3968class fast_sa(base_schedule):
3969    def init(self, **options):
3970        self.__dict__.update(options)
3971
3972    def update_guess(self, x0):
3973        x0 = asarray(x0)
3974        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
3975        T = self.T
3976        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
3977        xnew = xc*(self.upper - self.lower)+self.lower
3978        return xnew
3979
3980    def update_temp(self):
3981        self.T = self.T0*exp(-self.c * self.k**(self.quench))
3982        self.k += 1
3983        return
3984
3985class log_sa(base_schedule):        #OK
3986
3987    def init(self,**options):
3988        self.__dict__.update(options)
3989       
3990    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
3991        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
3992        T = self.T
3993        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
3994        xnew = xc*(self.upper - self.lower)+self.lower
3995        return xnew
3996       
3997    def update_temp(self):
3998        self.k += 1
3999        self.T = self.T0*self.slope**self.k
4000       
4001class _state(object):
4002    def __init__(self):
4003        self.x = None
4004        self.cost = None
4005
4006def makeTsched(data):
4007    if data['Algorithm'] == 'fast':
4008        sched = fast_sa()
4009        sched.quench = data['fast parms'][0]
4010        sched.c = data['fast parms'][1]
4011    elif data['Algorithm'] == 'log':
4012        sched = log_sa()
4013        sched.slope = data['log slope']
4014    sched.T0 = data['Annealing'][0]
4015    if not sched.T0:
4016        sched.T0 = 50.
4017    Tf = data['Annealing'][1]
4018    if not Tf:
4019        Tf = 0.001
4020    Tsched = [sched.T0,]
4021    while Tsched[-1] > Tf:
4022        sched.update_temp()
4023        Tsched.append(sched.T)
4024    return Tsched[1:]
4025   
4026def anneal(func, x0, args=(), schedule='fast', 
4027           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4028           feps=1e-6, quench=1.0, c=1.0,
4029           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4030           ranRange=0.10,autoRan=False,dlg=None):
4031    """Minimize a function using simulated annealing.
4032
4033    Schedule is a schedule class implementing the annealing schedule.
4034    Available ones are 'fast', 'cauchy', 'boltzmann'
4035
4036    :param callable func: f(x, \*args)
4037        Function to be optimized.
4038    :param ndarray x0:
4039        Initial guess.
4040    :param tuple args:
4041        Extra parameters to `func`.
4042    :param base_schedule schedule:
4043        Annealing schedule to use (a class).
4044    :param float T0:
4045        Initial Temperature (estimated as 1.2 times the largest
4046        cost-function deviation over random points in the range).
4047    :param float Tf:
4048        Final goal temperature.
4049    :param int maxeval:
4050        Maximum function evaluations.
4051    :param int maxaccept:
4052        Maximum changes to accept.
4053    :param int maxiter:
4054        Maximum cooling iterations.
4055    :param float feps:
4056        Stopping relative error tolerance for the function value in
4057        last four coolings.
4058    :param float quench,c:
4059        Parameters to alter fast_sa schedule.
4060    :param float/ndarray lower,upper:
4061        Lower and upper bounds on `x`.
4062    :param int dwell:
4063        The number of times to search the space at each temperature.
4064    :param float slope:
4065        Parameter for log schedule
4066    :param bool ranStart=False:
4067        True for set 10% of ranges about x
4068
4069    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4070
4071     * xmin (ndarray): Point giving smallest value found.
4072     * Jmin (float): Minimum value of function found.
4073     * T (float): Final temperature.
4074     * feval (int): Number of function evaluations.
4075     * iters (int): Number of cooling iterations.
4076     * accept (int): Number of tests accepted.
4077     * retval (int): Flag indicating stopping condition:
4078
4079              *  0: Points no longer changing
4080              *  1: Cooled to final temperature
4081              *  2: Maximum function evaluations
4082              *  3: Maximum cooling iterations reached
4083              *  4: Maximum accepted query locations reached
4084              *  5: Final point not the minimum amongst encountered points
4085
4086    *Notes*:
4087    Simulated annealing is a random algorithm which uses no derivative
4088    information from the function being optimized. In practice it has
4089    been more useful in discrete optimization than continuous
4090    optimization, as there are usually better algorithms for continuous
4091    optimization problems.
4092
4093    Some experimentation by trying the difference temperature
4094    schedules and altering their parameters is likely required to
4095    obtain good performance.
4096
4097    The randomness in the algorithm comes from random sampling in numpy.
4098    To obtain the same results you can call numpy.random.seed with the
4099    same seed immediately before calling scipy.optimize.anneal.
4100
4101    We give a brief description of how the three temperature schedules
4102    generate new points and vary their temperature. Temperatures are
4103    only updated with iterations in the outer loop. The inner loop is
4104    over xrange(dwell), and new points are generated for
4105    every iteration in the inner loop. (Though whether the proposed
4106    new points are accepted is probabilistic.)
4107
4108    For readability, let d denote the dimension of the inputs to func.
4109    Also, let x_old denote the previous state, and k denote the
4110    iteration number of the outer loop. All other variables not
4111    defined below are input variables to scipy.optimize.anneal itself.
4112
4113    In the 'fast' schedule the updates are ::
4114
4115        u ~ Uniform(0, 1, size=d)
4116        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4117        xc = y * (upper - lower)
4118        x_new = x_old + xc
4119
4120        T_new = T0 * exp(-c * k**quench)
4121
4122
4123    """
4124    x0 = asarray(x0)
4125    lower = asarray(lower)
4126    upper = asarray(upper)
4127
4128    schedule = eval(schedule+'_sa()')
4129    #   initialize the schedule
4130    schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4131        c=c, quench=quench, dwell=dwell, slope=slope)
4132
4133    current_state, last_state, best_state = _state(), _state(), _state()
4134    if ranStart:
4135        schedule.set_range(x0,ranRange)
4136    if T0 is None:
4137        x0 = schedule.getstart_temp(best_state)
4138    else:
4139        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4140        best_state.x = None
4141        best_state.cost = numpy.Inf
4142
4143    last_state.x = asarray(x0).copy()
4144    fval = func(x0,*args)
4145    schedule.feval += 1
4146    last_state.cost = fval
4147    if last_state.cost < best_state.cost:
4148        best_state.cost = fval
4149        best_state.x = asarray(x0).copy()
4150    schedule.T = schedule.T0
4151    fqueue = [100, 300, 500, 700]
4152    iters = 1
4153    keepGoing = True
4154    bestn = 0
4155    while keepGoing:
4156        retval = 0
4157        for n in xrange(dwell):
4158            current_state.x = schedule.update_guess(last_state.x)
4159            current_state.cost = func(current_state.x,*args)
4160            schedule.feval += 1
4161
4162            dE = current_state.cost - last_state.cost
4163            if schedule.accept_test(dE):
4164                last_state.x = current_state.x.copy()
4165                last_state.cost = current_state.cost
4166                if last_state.cost < best_state.cost:
4167                    best_state.x = last_state.x.copy()
4168                    best_state.cost = last_state.cost
4169                    bestn = n
4170                    if best_state.cost < 1.0 and autoRan:
4171                        schedule.set_range(x0,best_state.cost/2.)                       
4172        if dlg:
4173            GoOn = dlg.Update(min(100.,best_state.cost*100),
4174                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4175                    'Best trial:',bestn,  \
4176                    'MC/SA Residual:',best_state.cost*100,'%', \
4177                    ))[0]
4178            if not GoOn:
4179                best_state.x = last_state.x.copy()
4180                best_state.cost = last_state.cost
4181                retval = 5
4182        schedule.update_temp()
4183        iters += 1
4184        # Stopping conditions
4185        # 0) last saved values of f from each cooling step
4186        #     are all very similar (effectively cooled)
4187        # 1) Tf is set and we are below it
4188        # 2) maxeval is set and we are past it
4189        # 3) maxiter is set and we are past it
4190        # 4) maxaccept is set and we are past it
4191        # 5) user canceled run via progress bar
4192
4193        fqueue.append(squeeze(last_state.cost))
4194        fqueue.pop(0)
4195        af = asarray(fqueue)*1.0
4196        if retval == 5:
4197            print ' User terminated run; incomplete MC/SA'
4198            keepGoing = False
4199            break
4200        if all(abs((af-af[0])/af[0]) < feps):
4201            retval = 0
4202            if abs(af[-1]-best_state.cost) > feps*10:
4203                retval = 5
4204                print " Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1)
4205            break
4206        if (Tf is not None) and (schedule.T < Tf):
4207#            print ' Minimum T reached in %d steps'%(iters-1)
4208            retval = 1
4209            break
4210        if (maxeval is not None) and (schedule.feval > maxeval):
4211            retval = 2
4212            break
4213        if (iters > maxiter):
4214            print  " Warning: Maximum number of iterations exceeded."
4215            retval = 3
4216            break
4217        if (maxaccept is not None) and (schedule.accepted > maxaccept):
4218            retval = 4
4219            break
4220
4221    return best_state.x, best_state.cost, schedule.T, \
4222           schedule.feval, iters, schedule.accepted, retval
4223
4224def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4225    outlist = []
4226    timelist = []
4227    nsflist = []
4228    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
4229    for n in range(iCyc):
4230        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False)         #mcsa result,time,rcov
4231        outlist.append(result[0])
4232        timelist.append(result[1])
4233        nsflist.append(result[2])
4234        print ' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
4235    out_q.put(outlist)
4236    out_t.put(timelist)
4237    out_n.put(nsflist)
4238
4239def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4240    import multiprocessing as mp
4241   
4242    out_q = mp.Queue()
4243    out_t = mp.Queue()
4244    out_n = mp.Queue()
4245    procs = []
4246    totsftime = 0.
4247    totnsf = 0
4248    iCyc = np.zeros(nprocs)
4249    for i in range(nCyc):
4250        iCyc[i%nprocs] += 1
4251    for i in range(nprocs):
4252        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4253        procs.append(p)
4254        p.start()
4255    resultlist = []
4256    for i in range(nprocs):
4257        resultlist += out_q.get()
4258        totsftime += np.sum(out_t.get())
4259        totnsf += np.sum(out_n.get())
4260    for p in procs:
4261        p.join()
4262    return resultlist,totsftime,totnsf
4263
4264def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4265    '''default doc string
4266   
4267    :param type name: description
4268   
4269    :returns: type name: description
4270    '''
4271   
4272    class RandomDisplacementBounds(object):
4273        """random displacement with bounds"""
4274        def __init__(self, xmin, xmax, stepsize=0.5):
4275            self.xmin = xmin
4276            self.xmax = xmax
4277            self.stepsize = stepsize
4278   
4279        def __call__(self, x):
4280            """take a random step but ensure the new position is within the bounds"""
4281            while True:
4282                # this could be done in a much more clever way, but it will work for example purposes
4283                steps = self.xmax-self.xmin
4284                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
4285                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4286                    break
4287            return xnew
4288   
4289    global tsum,nsum
4290    tsum = 0.
4291    nsum = 0
4292   
4293    def getMDparms(item,pfx,parmDict,varyList):
4294        parmDict[pfx+'MDaxis'] = item['axis']
4295        parmDict[pfx+'MDval'] = item['Coef'][0]
4296        if item['Coef'][1]:
4297            varyList += [pfx+'MDval',]
4298            limits = item['Coef'][2]
4299            lower.append(limits[0])
4300            upper.append(limits[1])
4301                       
4302    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4303        parmDict[pfx+'Atype'] = item['atType']
4304        aTypes |= set([item['atType'],]) 
4305        pstr = ['Ax','Ay','Az']
4306        XYZ = [0,0,0]
4307        for i in range(3):
4308            name = pfx+pstr[i]
4309            parmDict[name] = item['Pos'][0][i]
4310            XYZ[i] = parmDict[name]
4311            if item['Pos'][1][i]:
4312                varyList += [name,]
4313                limits = item['Pos'][2][i]
4314                lower.append(limits[0])
4315                upper.append(limits[1])
4316        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4317           
4318    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4319        parmDict[mfx+'MolCent'] = item['MolCent']
4320        parmDict[mfx+'RBId'] = item['RBId']
4321        pstr = ['Px','Py','Pz']
4322        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4323        for i in range(3):
4324            name = mfx+pstr[i]
4325            parmDict[name] = item['Pos'][0][i]
4326            if item['Pos'][1][i]:
4327                varyList += [name,]
4328                limits = item['Pos'][2][i]
4329                lower.append(limits[0])
4330                upper.append(limits[1])
4331        AV = item['Ori'][0]
4332        A = AV[0]
4333        V = AV[1:]
4334        for i in range(4):
4335            name = mfx+ostr[i]
4336            if i:
4337                parmDict[name] = V[i-1]
4338            else:
4339                parmDict[name] = A
4340            if item['Ovar'] == 'AV':
4341                varyList += [name,]
4342                limits = item['Ori'][2][i]
4343                lower.append(limits[0])
4344                upper.append(limits[1])
4345            elif item['Ovar'] == 'A' and not i:
4346                varyList += [name,]
4347                limits = item['Ori'][2][i]
4348                lower.append(limits[0])
4349                upper.append(limits[1])
4350        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
4351            for i in range(len(item['Tor'][0])):
4352                name = mfx+'Tor'+str(i)
4353                parmDict[name] = item['Tor'][0][i]
4354                if item['Tor'][1][i]:
4355                    varyList += [name,]
4356                    limits = item['Tor'][2][i]
4357                    lower.append(limits[0])
4358                    upper.append(limits[1])
4359        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
4360        aTypes |= set(atypes)
4361        atNo += len(atypes)
4362        return atNo
4363               
4364    def GetAtomM(Xdata,SGData):
4365        Mdata = []
4366        for xyz in Xdata:
4367            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
4368        return np.array(Mdata)
4369       
4370    def GetAtomT(RBdata,parmDict):
4371        'Needs a doc string'
4372        atNo = parmDict['atNo']
4373        nfixAt = parmDict['nfixAt']
4374        Tdata = atNo*[' ',]
4375        for iatm in range(nfixAt):
4376            parm = ':'+str(iatm)+':Atype'
4377            if parm in parmDict:
4378                Tdata[iatm] = aTypes.index(parmDict[parm])
4379        iatm = nfixAt
4380        for iObj in range(parmDict['nObj']):
4381            pfx = str(iObj)+':'
4382            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4383                if parmDict[pfx+'Type'] == 'Vector':
4384                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4385                    nAtm = len(RBRes['rbVect'][0])
4386                else:       #Residue
4387                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4388                    nAtm = len(RBRes['rbXYZ'])
4389                for i in range(nAtm):
4390                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
4391                    iatm += 1
4392            elif parmDict[pfx+'Type'] == 'Atom':
4393                atNo = parmDict[pfx+'atNo']
4394                parm = pfx+'Atype'              #remove extra ':'
4395                if parm in parmDict:
4396                    Tdata[atNo] = aTypes.index(parmDict[parm])
4397                iatm += 1
4398            else:
4399                continue        #skips March Dollase
4400        return Tdata
4401       
4402    def GetAtomX(RBdata,parmDict):
4403        'Needs a doc string'
4404        Bmat = parmDict['Bmat']
4405        atNo = parmDict['atNo']
4406        nfixAt = parmDict['nfixAt']
4407        Xdata = np.zeros((3,atNo))
4408        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
4409        for iatm in range(nfixAt):
4410            for key in keys:
4411                parm = ':'+str(iatm)+key
4412                if parm in parmDict:
4413                    keys[key][iatm] = parmDict[parm]
4414        iatm = nfixAt
4415        for iObj in range(parmDict['nObj']):
4416            pfx = str(iObj)+':'
4417            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4418                if parmDict[pfx+'Type'] == 'Vector':
4419                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4420                    vecs = RBRes['rbVect']
4421                    mags = RBRes['VectMag']
4422                    Cart = np.zeros_like(vecs[0])
4423                    for vec,mag in zip(vecs,mags):
4424                        Cart += vec*mag
4425                elif parmDict[pfx+'Type'] == 'Residue':
4426                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4427                    Cart = np.array(RBRes['rbXYZ'])
4428                    for itor,seq in enumerate(RBRes['rbSeq']):
4429                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
4430                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
4431                if parmDict[pfx+'MolCent'][1]:
4432                    Cart -= parmDict[pfx+'MolCent'][0]
4433                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4434                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
4435                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
4436                iatm += len(Cart)
4437            elif parmDict[pfx+'Type'] == 'Atom':
4438                atNo = parmDict[pfx+'atNo']
4439                for key in keys:
4440                    parm = pfx+key[1:]              #remove extra ':'
4441                    if parm in parmDict:
4442                        keys[key][atNo] = parmDict[parm]
4443                iatm += 1
4444            else:
4445                continue        #skips March Dollase
4446        return Xdata.T
4447       
4448    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
4449        allX = np.inner(Xdata,SGM)+SGT
4450        allT = np.repeat(Tdata,allX.shape[1])
4451        allM = np.repeat(Mdata,allX.shape[1])
4452        allX = np.reshape(allX,(-1,3))
4453        return allT,allM,allX
4454
4455    def getAllX(Xdata,SGM,SGT):
4456        allX = np.inner(Xdata,SGM)+SGT
4457        allX = np.reshape(allX,(-1,3))
4458        return allX
4459       
4460    def normQuaternions(RBdata,parmDict,varyList,values):
4461        for iObj in range(parmDict['nObj']):
4462            pfx = str(iObj)+':'
4463            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4464                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4465                A,V = Q2AVdeg(Qori)
4466                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
4467                    if i:
4468                        parmDict[pfx+name] = V[i-1]
4469                    else:
4470                        parmDict[pfx+name] = A
4471       
4472    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
4473        ''' Compute structure factors for all h,k,l for phase
4474        input:
4475            refList: [ref] where each ref = h,k,l,m,d,...
4476            rcov:   array[nref,nref] covariance terms between Fo^2 values
4477            ifInv:  bool True if centrosymmetric
4478            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
4479            RBdata: [dict] rigid body dictionary
4480            varyList: [list] names of varied parameters in MC/SA (not used here)           
4481            ParmDict: [dict] problem parameters
4482        puts result F^2 in each ref[5] in refList
4483        returns:
4484            delt-F*rcov*delt-F/sum(Fo^2)
4485        '''   
4486           
4487        global tsum,nsum
4488        t0 = time.time()
4489        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
4490        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
4491        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
4492        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
4493        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
4494        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
4495        refList[5] = Aterm
4496        if not ifInv:
4497            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
4498        if len(cosTable):        #apply MD correction
4499            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
4500        sumFcsq = np.sum(refList[5])
4501        scale = parmDict['sumFosq']/sumFcsq
4502        refList[5] *= scale
4503        refList[6] = refList[4]-refList[5]
4504        M = np.inner(refList[6],np.inner(rcov,refList[6]))
4505        tsum += (time.time()-t0)
4506        nsum += 1
4507        return np.sqrt(M/np.sum(refList[4]**2))
4508   
4509    def MCSAcallback(x, f, fmin,accept):
4510        return not pgbar.Update(min(100.,fmin*100),
4511            newmsg='%s%8.4f%s'%('MC/SA Residual:',fmin*100,'%'))[0]
4512
4513    sq2pi = np.sqrt(2*np.pi)
4514    sq4pi = np.sqrt(4*np.pi)
4515    generalData = data['General']
4516    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
4517    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
4518    SGData = generalData['SGData']
4519    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
4520    SGMT = np.array([SGData['SGOps'][i][0