source: trunk/GSASIImath.py @ 3044

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

remove local basinhopping

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 183.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2017-09-03 19:35:45 +0000 (Sun, 03 Sep 2017) $
5# $Author: vondreele $
6# $Revision: 3044 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 3044 2017-09-03 19:35:45Z vondreele $
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: 3044 $")
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,old):
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 from errat.f
2577    b1_old = np.array([ 
2578        [1154.343,  600.213, 1051.018, 1132.885,  960.738],
2579        [600.213, 1286.818, 1282.042,  957.156,  612.789],
2580        [1051.018, 1282.042, 3519.471,  991.974, 1226.491],
2581        [1132.885,  957.156,  991.974, 1798.672,  820.355],
2582        [960.738,  612.789, 1226.491,  820.355, 2428.966]
2583        ])
2584    avg_old = np.array([ 0.225, 0.281, 0.071, 0.237, 0.044])    #Table 1 3.5A Obsd. Fr. p 1513
2585# data taken from erratv2.ccp
2586    b1 = np.array([
2587          [5040.279078850848200,        3408.805141583649400,   4152.904423767300600,   4236.200004171890200,   5054.781210204625500], 
2588          [3408.805141583648900,        8491.906094010220800,   5958.881777877950300,   1521.387352718486200,   4304.078200827221700], 
2589          [4152.904423767301500,        5958.881777877952100,   7637.167089335050100,   6620.715738223072500,   5287.691183798410700], 
2590          [4236.200004171890200,        1521.387352718486200,   6620.715738223072500,   18368.343774298410000,  4050.797811118806700], 
2591          [5054.781210204625500,        4304.078200827220800,   5287.691183798409800,   4050.797811118806700,   6666.856740479164700]])
2592    avg = np.array([0.192765509919262, 0.195575208778518, 0.275322406824210, 0.059102357035642, 0.233154192767480])
2593    General = Phase['General']
2594    Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
2595    cx,ct,cs,cia = General['AtomPtrs']
2596    Atoms = Phase['Atoms']
2597    cartAtoms = []
2598    xyzmin = 999.*np.ones(3)
2599    xyzmax = -999.*np.ones(3)
2600    #select residue atoms,S,Se --> O make cartesian
2601    for atom in Atoms:
2602        if atom[1] in resNames:
2603            if atom[4].strip() in ['S','Se']:
2604                if not old:
2605                    continue        #S,Se skipped for erratv2?
2606                atom[3] = 'Os'
2607                atom[4] = 'O'
2608            cartAtoms.append(atom[:cx+3])
2609            cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
2610    XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
2611    xyzmin = np.array([np.min(XYZ.T[i]) for i in [0,1,2]])
2612    xyzmax = np.array([np.max(XYZ.T[i]) for i in [0,1,2]])
2613    nbox = list(np.array(np.ceil((xyzmax-xyzmin)/4.),dtype=int))+[15,]
2614    Boxes = np.zeros(nbox,dtype=int)
2615    iBox = np.array([np.trunc((XYZ.T[i]-xyzmin[i])/4.) for i in [0,1,2]],dtype=int).T
2616    for ib,box in enumerate(iBox):  #put in a try for too many atoms in box (IndexError)?
2617        Boxes[box[0],box[1],box[2],0] += 1
2618        Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
2619    #Box content checks with errat.f ibox1 array
2620    indices = (-1,0,1)
2621    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) 
2622    dsmax = 3.75**2
2623    if old:
2624        dsmax = 3.5**2
2625    chains = []
2626    resIntAct = []
2627    chainIntAct = []
2628    res = []
2629    resNames = []
2630    resname = []
2631    newChain = True
2632    intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2633    for ia,atom in enumerate(cartAtoms):
2634        if atom[2] not in chains:   #get chain id & save residue sequence from last chain
2635            chains.append(atom[2])
2636            if len(resIntAct):
2637                resIntAct.append(sumintact(intact))
2638                chainIntAct.append(resIntAct)
2639                resNames += resname
2640                res = []
2641                resname = []
2642                resIntAct = []
2643                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2644                newChain = True
2645        if atom[0] not in res:  #new residue, get residue no.
2646            res.append(atom[0])
2647            resname.append('%s-%s%s'%(atom[2],atom[0],atom[1]))
2648            if not newChain:
2649                resIntAct.append(sumintact(intact))
2650            intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2651            newChain = False
2652        ibox = iBox[ia]         #box location of atom
2653        tgts = []
2654        for unit in Units:      #assemble list of all possible target atoms
2655            jbox = ibox+unit
2656            if np.all(jbox>=0) and np.all((jbox-nbox[:3])<0):               
2657                tgts += list(Boxes[jbox[0],jbox[1],jbox[2]])
2658        tgts = list(set(tgts))
2659        tgts = [tgt for tgt in tgts if np.sum((XYZ[ia]-XYZ[tgt])**2) < dsmax]
2660        tgts = [tgt for tgt in tgts if atom[:3] != cartAtoms[tgt][:3]]    #exclude same residue
2661        ires = int(atom[0])
2662        if old:
2663            if atom[3].strip() == 'C':
2664                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2665            elif atom[3].strip() == 'N':
2666                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() in ['C','CA'] and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2667            elif atom[3].strip() == 'CA':
2668                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2669        else:
2670            if atom[3].strip() == 'C':
2671                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) == ires+1)]
2672            elif atom[3].strip() == 'N':
2673                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'C' and int(cartAtoms[tgt][0]) == ires-1)]
2674        for tgt in tgts:
2675            dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
2676            mult = 1.0
2677            if dsqt > 3.25 and not old:
2678                mult = 2.*(3.75-dsqt)
2679            intype = atom[4].strip()+cartAtoms[tgt][4].strip()
2680            intact[intype] += mult
2681    resNames += resname
2682    resIntAct.append(sumintact(intact))
2683    chainIntAct.append(resIntAct)
2684    chainProb = []
2685    for ich,chn in enumerate(chains):
2686        IntAct = chainIntAct[ich]
2687        nRes = len(IntAct)
2688        Probs = [0.,0.,0.,0.]
2689        for i in range(4,nRes-4):
2690            mtrx = np.zeros(5)
2691            summ = 0.
2692            for j in range(i-4,i+5):
2693                summ += np.sum(np.array(IntAct[j].values()))
2694                if old:
2695                    mtrx[0] += IntAct[j]['CC']
2696                    mtrx[1] += IntAct[j]['CO']
2697                    mtrx[2] += IntAct[j]['NN']
2698                    mtrx[3] += IntAct[j]['NO']
2699                    mtrx[4] += IntAct[j]['OO']
2700                else:
2701                    mtrx[0] += IntAct[j]['CC']
2702                    mtrx[1] += IntAct[j]['CN']
2703                    mtrx[2] += IntAct[j]['CO']
2704                    mtrx[3] += IntAct[j]['NN']
2705                    mtrx[4] += IntAct[j]['NO']
2706            mtrx /= summ
2707            if old:
2708                mtrx -= avg_old
2709                prob = np.inner(np.inner(mtrx,b1_old),mtrx)
2710            else:
2711                mtrx -= avg
2712                prob = np.inner(np.inner(mtrx,b1),mtrx)
2713            Probs.append(prob)
2714            print i,prob
2715        Probs += 4*[0.,]
2716        chainProb += Probs
2717    return resNames,chainProb
2718   
2719################################################################################
2720##### Texture fitting stuff
2721################################################################################
2722
2723def FitTexture(General,Gangls,refData,keyList,pgbar):
2724    import pytexture as ptx
2725    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2726   
2727    def printSpHarm(textureData,SHtextureSig):
2728        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
2729        names = ['omega','chi','phi']
2730        namstr = '  names :'
2731        ptstr =  '  values:'
2732        sigstr = '  esds  :'
2733        for name in names:
2734            namstr += '%12s'%('Sample '+name)
2735            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2736            if 'Sample '+name in SHtextureSig:
2737                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2738            else:
2739                sigstr += 12*' '
2740        print namstr
2741        print ptstr
2742        print sigstr
2743        print '\n Texture coefficients:'
2744        SHcoeff = textureData['SH Coeff'][1]
2745        SHkeys = SHcoeff.keys()
2746        nCoeff = len(SHcoeff)
2747        nBlock = nCoeff/10+1
2748        iBeg = 0
2749        iFin = min(iBeg+10,nCoeff)
2750        for block in range(nBlock):
2751            namstr = '  names :'
2752            ptstr =  '  values:'
2753            sigstr = '  esds  :'
2754            for name in SHkeys[iBeg:iFin]:
2755                if 'C' in name:
2756                    namstr += '%12s'%(name)
2757                    ptstr += '%12.3f'%(SHcoeff[name])
2758                    if name in SHtextureSig:
2759                        sigstr += '%12.3f'%(SHtextureSig[name])
2760                    else:
2761                        sigstr += 12*' '
2762            print namstr
2763            print ptstr
2764            print sigstr
2765            iBeg += 10
2766            iFin = min(iBeg+10,nCoeff)
2767           
2768    def Dict2Values(parmdict, varylist):
2769        '''Use before call to leastsq to setup list of values for the parameters
2770        in parmdict, as selected by key in varylist'''
2771        return [parmdict[key] for key in varylist] 
2772       
2773    def Values2Dict(parmdict, varylist, values):
2774        ''' Use after call to leastsq to update the parameter dictionary with
2775        values corresponding to keys in varylist'''
2776        parmdict.update(zip(varylist,values))
2777       
2778    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2779        parmDict.update(zip(varyList,values))
2780        Mat = np.empty(0)
2781        sumObs = 0
2782        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2783        for hist in Gangls.keys():
2784            Refs = refData[hist]
2785            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2786            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2787#            wt = 1./np.max(Refs[:,4],.25)
2788            sumObs += np.sum(wt*Refs[:,5])
2789            Refs[:,6] = 1.
2790            H = Refs[:,:3]
2791            phi,beta = G2lat.CrsAng(H,cell,SGData)
2792            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2793            for item in parmDict:
2794                if 'C' in item:
2795                    L,M,N = eval(item.strip('C'))
2796                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2797                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2798                    Lnorm = G2lat.Lnorm(L)
2799                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2800            mat = wt*(Refs[:,5]-Refs[:,6])
2801            Mat = np.concatenate((Mat,mat))
2802        sumD = np.sum(np.abs(Mat))
2803        R = min(100.,100.*sumD/sumObs)
2804        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2805        print ' Residual: %.3f%%'%(R)
2806        return Mat
2807       
2808    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2809        Mat = np.empty(0)
2810        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2811        for hist in Gangls.keys():
2812            mat = np.zeros((len(varyList),len(refData[hist])))
2813            Refs = refData[hist]
2814            H = Refs[:,:3]
2815            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2816#            wt = 1./np.max(Refs[:,4],.25)
2817            phi,beta = G2lat.CrsAng(H,cell,SGData)
2818            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2819            for j,item in enumerate(varyList):
2820                if 'C' in item:
2821                    L,M,N = eval(item.strip('C'))
2822                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2823                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2824                    Lnorm = G2lat.Lnorm(L)
2825                    mat[j] = -wt*Lnorm*Kcl*Ksl
2826                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2827                        try:
2828                            l = varyList.index(itema)
2829                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2830                        except ValueError:
2831                            pass
2832            if len(Mat):
2833                Mat = np.concatenate((Mat,mat.T))
2834            else:
2835                Mat = mat.T
2836        print 'deriv'
2837        return Mat
2838
2839    print ' Fit texture for '+General['Name']
2840    SGData = General['SGData']
2841    cell = General['Cell'][1:7]
2842    Texture = General['SH Texture']
2843    if not Texture['Order']:
2844        return 'No spherical harmonics coefficients'
2845    varyList = []
2846    parmDict = copy.copy(Texture['SH Coeff'][1])
2847    for item in ['Sample omega','Sample chi','Sample phi']:
2848        parmDict[item] = Texture[item][1]
2849        if Texture[item][0]:
2850            varyList.append(item)
2851    if Texture['SH Coeff'][0]:
2852        varyList += Texture['SH Coeff'][1].keys()
2853    while True:
2854        begin = time.time()
2855        values =  np.array(Dict2Values(parmDict, varyList))
2856        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
2857            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
2858        ncyc = int(result[2]['nfev']/2)
2859        if ncyc:
2860            runtime = time.time()-begin   
2861            chisq = np.sum(result[2]['fvec']**2)
2862            Values2Dict(parmDict, varyList, result[0])
2863            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
2864            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
2865            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2866            try:
2867                sig = np.sqrt(np.diag(result[1])*GOF)
2868                if np.any(np.isnan(sig)):
2869                    print '*** Least squares aborted - some invalid esds possible ***'
2870                break                   #refinement succeeded - finish up!
2871            except ValueError:          #result[1] is None on singular matrix
2872                print '**** Refinement failed - singular matrix ****'
2873                return None
2874        else:
2875            break
2876   
2877    if ncyc:
2878        for parm in parmDict:
2879            if 'C' in parm:
2880                Texture['SH Coeff'][1][parm] = parmDict[parm]
2881            else:
2882                Texture[parm][1] = parmDict[parm] 
2883        sigDict = dict(zip(varyList,sig))
2884        printSpHarm(Texture,sigDict)
2885       
2886    return None
2887   
2888################################################################################
2889##### Fourier & charge flip stuff
2890################################################################################
2891
2892def adjHKLmax(SGData,Hmax):
2893    '''default doc string
2894   
2895    :param type name: description
2896   
2897    :returns: type name: description
2898   
2899    '''
2900    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2901        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
2902        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
2903        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2904    else:
2905        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
2906        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
2907        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2908
2909def OmitMap(data,reflDict,pgbar=None):
2910    '''default doc string
2911   
2912    :param type name: description
2913   
2914    :returns: type name: description
2915   
2916    '''
2917    generalData = data['General']
2918    if not generalData['Map']['MapType']:
2919        print '**** ERROR - Fourier map not defined'
2920        return
2921    mapData = generalData['Map']
2922    dmin = mapData['Resolution']
2923    SGData = generalData['SGData']
2924    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2925    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2926    cell = generalData['Cell'][1:8]       
2927    A = G2lat.cell2A(cell[:6])
2928    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2929    adjHKLmax(SGData,Hmax)
2930    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2931    time0 = time.time()
2932    for iref,ref in enumerate(reflDict['RefList']):
2933        if ref[4] >= dmin:
2934            Fosq,Fcsq,ph = ref[8:11]
2935            Uniq = np.inner(ref[:3],SGMT)
2936            Phi = np.inner(ref[:3],SGT)
2937            for i,hkl in enumerate(Uniq):        #uses uniq
2938                hkl = np.asarray(hkl,dtype='i')
2939                dp = 360.*Phi[i]                #and phi
2940                a = cosd(ph+dp)
2941                b = sind(ph+dp)
2942                phasep = complex(a,b)
2943                phasem = complex(a,-b)
2944                if '2Fo-Fc' in mapData['MapType']:
2945                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
2946                else:
2947                    F = np.sqrt(Fosq)
2948                h,k,l = hkl+Hmax
2949                Fhkl[h,k,l] = F*phasep
2950                h,k,l = -hkl+Hmax
2951                Fhkl[h,k,l] = F*phasem
2952    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2953    M = np.mgrid[0:4,0:4,0:4]
2954    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
2955    iBeg = blkIds*rho0.shape/4
2956    iFin = (blkIds+1)*rho0.shape/4
2957    rho_omit = np.zeros_like(rho0)
2958    nBlk = 0
2959    for iB,iF in zip(iBeg,iFin):
2960        rho1 = np.copy(rho0)
2961        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
2962        Fnew = fft.ifftshift(fft.ifftn(rho1))
2963        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
2964        phase = Fnew/np.absolute(Fnew)
2965        OFhkl = np.absolute(Fhkl)*phase
2966        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
2967        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]])
2968        nBlk += 1
2969        pgbar.Update(nBlk)
2970    mapData['rho'] = np.real(rho_omit)/cell[6]
2971    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2972    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2973    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2974    return mapData
2975   
2976def FourierMap(data,reflDict):
2977    '''default doc string
2978   
2979    :param type name: description
2980   
2981    :returns: type name: description
2982   
2983    '''
2984    generalData = data['General']
2985    mapData = generalData['Map']
2986    dmin = mapData['Resolution']
2987    SGData = generalData['SGData']
2988    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2989    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2990    cell = generalData['Cell'][1:8]       
2991    A = G2lat.cell2A(cell[:6])
2992    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2993    adjHKLmax(SGData,Hmax)
2994    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2995#    Fhkl[0,0,0] = generalData['F000X']
2996    time0 = time.time()
2997    for iref,ref in enumerate(reflDict['RefList']):
2998        if ref[4] > dmin:
2999            Fosq,Fcsq,ph = ref[8:11]
3000            Uniq = np.inner(ref[:3],SGMT)
3001            Phi = np.inner(ref[:3],SGT)
3002            for i,hkl in enumerate(Uniq):        #uses uniq
3003                hkl = np.asarray(hkl,dtype='i')
3004                dp = 360.*Phi[i]                #and phi
3005                a = cosd(ph+dp)
3006                b = sind(ph+dp)
3007                phasep = complex(a,b)
3008                phasem = complex(a,-b)
3009                if 'Fobs' in mapData['MapType']:
3010                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
3011                    h,k,l = hkl+Hmax
3012                    Fhkl[h,k,l] = F*phasep
3013                    h,k,l = -hkl+Hmax
3014                    Fhkl[h,k,l] = F*phasem
3015                elif 'Fcalc' in mapData['MapType']:
3016                    F = np.sqrt(Fcsq)
3017                    h,k,l = hkl+Hmax
3018                    Fhkl[h,k,l] = F*phasep
3019                    h,k,l = -hkl+Hmax
3020                    Fhkl[h,k,l] = F*phasem
3021                elif 'delt-F' in mapData['MapType']:
3022                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3023                    h,k,l = hkl+Hmax
3024                    Fhkl[h,k,l] = dF*phasep
3025                    h,k,l = -hkl+Hmax
3026                    Fhkl[h,k,l] = dF*phasem
3027                elif '2*Fo-Fc' in mapData['MapType']:
3028                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3029                    h,k,l = hkl+Hmax
3030                    Fhkl[h,k,l] = F*phasep
3031                    h,k,l = -hkl+Hmax
3032                    Fhkl[h,k,l] = F*phasem
3033                elif 'Patterson' in mapData['MapType']:
3034                    h,k,l = hkl+Hmax
3035                    Fhkl[h,k,l] = complex(Fosq,0.)
3036                    h,k,l = -hkl+Hmax
3037                    Fhkl[h,k,l] = complex(Fosq,0.)
3038    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3039    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
3040    mapData['Type'] = reflDict['Type']
3041    mapData['rho'] = np.real(rho)
3042    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3043    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3044   
3045def Fourier4DMap(data,reflDict):
3046    '''default doc string
3047   
3048    :param type name: description
3049   
3050    :returns: type name: description
3051   
3052    '''
3053    generalData = data['General']
3054    map4DData = generalData['4DmapData']
3055    mapData = generalData['Map']
3056    dmin = mapData['Resolution']
3057    SGData = generalData['SGData']
3058    SSGData = generalData['SSGData']
3059    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3060    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3061    cell = generalData['Cell'][1:8]       
3062    A = G2lat.cell2A(cell[:6])
3063    maxM = 4
3064    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
3065    adjHKLmax(SGData,Hmax)
3066    Hmax = np.asarray(Hmax,dtype='i')+1
3067    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3068    time0 = time.time()
3069    for iref,ref in enumerate(reflDict['RefList']):
3070        if ref[5] > dmin:
3071            Fosq,Fcsq,ph = ref[9:12]
3072            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
3073            Uniq = np.inner(ref[:4],SSGMT)
3074            Phi = np.inner(ref[:4],SSGT)
3075            for i,hkl in enumerate(Uniq):        #uses uniq
3076                hkl = np.asarray(hkl,dtype='i')
3077                dp = 360.*Phi[i]                #and phi
3078                a = cosd(ph+dp)
3079                b = sind(ph+dp)
3080                phasep = complex(a,b)
3081                phasem = complex(a,-b)
3082                if 'Fobs' in mapData['MapType']:
3083                    F = np.sqrt(Fosq)
3084                    h,k,l,m = hkl+Hmax
3085                    Fhkl[h,k,l,m] = F*phasep
3086                    h,k,l,m = -hkl+Hmax
3087                    Fhkl[h,k,l,m] = F*phasem
3088                elif 'Fcalc' in mapData['MapType']:
3089                    F = np.sqrt(Fcsq)
3090                    h,k,l,m = hkl+Hmax
3091                    Fhkl[h,k,l,m] = F*phasep
3092                    h,k,l,m = -hkl+Hmax
3093                    Fhkl[h,k,l,m] = F*phasem                   
3094                elif 'delt-F' in mapData['MapType']:
3095                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
3096                    h,k,l,m = hkl+Hmax
3097                    Fhkl[h,k,l,m] = dF*phasep
3098                    h,k,l,m = -hkl+Hmax
3099                    Fhkl[h,k,l,m] = dF*phasem
3100    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
3101    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
3102    map4DData['rho'] = np.real(SSrho)
3103    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3104    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3105    map4DData['Type'] = reflDict['Type']
3106    mapData['Type'] = reflDict['Type']
3107    mapData['rho'] = np.real(rho)
3108    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3109    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3110    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
3111
3112# map printing for testing purposes
3113def printRho(SGLaue,rho,rhoMax):                         
3114    '''default doc string
3115   
3116    :param type name: description
3117   
3118    :returns: type name: description
3119   
3120    '''
3121    dim = len(rho.shape)
3122    if dim == 2:
3123        ix,jy = rho.shape
3124        for j in range(jy):
3125            line = ''
3126            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3127                line += (jy-j)*'  '
3128            for i in range(ix):
3129                r = int(100*rho[i,j]/rhoMax)
3130                line += '%4d'%(r)
3131            print line+'\n'
3132    else:
3133        ix,jy,kz = rho.shape
3134        for k in range(kz):
3135            print 'k = ',k
3136            for j in range(jy):
3137                line = ''
3138                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3139                    line += (jy-j)*'  '
3140                for i in range(ix):
3141                    r = int(100*rho[i,j,k]/rhoMax)
3142                    line += '%4d'%(r)
3143                print line+'\n'
3144## keep this
3145               
3146def findOffset(SGData,A,Fhkl):   
3147    '''default doc string
3148   
3149    :param type name: description
3150   
3151    :returns: type name: description
3152   
3153    '''
3154    if SGData['SpGrp'] == 'P 1':
3155        return [0,0,0]   
3156    hklShape = Fhkl.shape
3157    hklHalf = np.array(hklShape)/2
3158    sortHKL = np.argsort(Fhkl.flatten())
3159    Fdict = {}
3160    for hkl in sortHKL:
3161        HKL = np.unravel_index(hkl,hklShape)
3162        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
3163        if F == 0.:
3164            break
3165        Fdict['%.6f'%(np.absolute(F))] = hkl
3166    Flist = np.flipud(np.sort(Fdict.keys()))
3167    F = str(1.e6)
3168    i = 0
3169    DH = []
3170    Dphi = []
3171    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3172    for F in Flist:
3173        hkl = np.unravel_index(Fdict[F],hklShape)
3174        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
3175            continue
3176        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
3177        Uniq = np.array(Uniq,dtype='i')
3178        Phi = np.array(Phi)
3179        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
3180        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3181        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
3182        ang0 = np.angle(Fh0,deg=True)/360.
3183        for H,phi in zip(Uniq,Phi)[1:]:
3184            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
3185            dH = H-hkl
3186            dang = ang-ang0
3187            DH.append(dH)
3188            Dphi.append((dang+.5) % 1.0)
3189        if i > 20 or len(DH) > 30:
3190            break
3191        i += 1
3192    DH = np.array(DH)
3193    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
3194    Dphi = np.array(Dphi)
3195    steps = np.array(hklShape)
3196    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
3197    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
3198    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
3199    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3200    hist,bins = np.histogram(Mmap,bins=1000)
3201#    for i,item in enumerate(hist[:10]):
3202#        print item,bins[i]
3203    chisq = np.min(Mmap)
3204    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3205    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
3206#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3207    return DX
3208   
3209def ChargeFlip(data,reflDict,pgbar):
3210    '''default doc string
3211   
3212    :param type name: description
3213   
3214    :returns: type name: description
3215   
3216    '''
3217    generalData = data['General']
3218    mapData = generalData['Map']
3219    flipData = generalData['Flip']
3220    FFtable = {}
3221    if 'None' not in flipData['Norm element']:
3222        normElem = flipData['Norm element'].upper()
3223        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3224        for ff in FFs:
3225            if ff['Symbol'] == normElem:
3226                FFtable.update(ff)
3227    dmin = flipData['Resolution']
3228    SGData = generalData['SGData']
3229    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3230    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3231    cell = generalData['Cell'][1:8]       
3232    A = G2lat.cell2A(cell[:6])
3233    Vol = cell[6]
3234    im = 0
3235    if generalData['Type'] in ['modulated','magnetic',]:
3236        im = 1
3237    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3238    adjHKLmax(SGData,Hmax)
3239    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3240    time0 = time.time()
3241    for iref,ref in enumerate(reflDict['RefList']):
3242        dsp = ref[4+im]
3243        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
3244            continue
3245        if dsp > dmin:
3246            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3247            if FFtable:
3248                SQ = 0.25/dsp**2
3249                ff *= G2el.ScatFac(FFtable,SQ)[0]
3250            if ref[8+im] > 0.:         #use only +ve Fobs**2
3251                E = np.sqrt(ref[8+im])/ff
3252            else:
3253                E = 0.
3254            ph = ref[10]
3255            ph = rn.uniform(0.,360.)
3256            Uniq = np.inner(ref[:3],SGMT)
3257            Phi = np.inner(ref[:3],SGT)
3258            for i,hkl in enumerate(Uniq):        #uses uniq
3259                hkl = np.asarray(hkl,dtype='i')
3260                dp = 360.*Phi[i]                #and phi
3261                a = cosd(ph+dp)
3262                b = sind(ph+dp)
3263                phasep = complex(a,b)
3264                phasem = complex(a,-b)
3265                h,k,l = hkl+Hmax
3266                Ehkl[h,k,l] = E*phasep
3267                h,k,l = -hkl+Hmax
3268                Ehkl[h,k,l] = E*phasem
3269#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3270    testHKL = np.array(flipData['testHKL'])+Hmax
3271    CEhkl = copy.copy(Ehkl)
3272    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3273    Emask = ma.getmask(MEhkl)
3274    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3275    Ncyc = 0
3276    old = np.seterr(all='raise')
3277    twophases = []
3278    while True:       
3279        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3280        CEsig = np.std(CErho)
3281        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3282        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3283        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3284        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3285        phase = CFhkl/np.absolute(CFhkl)
3286        twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
3287        CEhkl = np.absolute(Ehkl)*phase
3288        Ncyc += 1
3289        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3290        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3291        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3292        if Rcf < 5.:
3293            break
3294        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3295        if not GoOn or Ncyc > 10000:
3296            break
3297    np.seterr(**old)
3298    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3299    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
3300    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3301    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3302       
3303    mapData['Rcf'] = Rcf
3304    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3305    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3306    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3307    mapData['Type'] = reflDict['Type']
3308    return mapData,twophases
3309   
3310def findSSOffset(SGData,SSGData,A,Fhklm):   
3311    '''default doc string
3312   
3313    :param type name: description
3314   
3315    :returns: type name: description
3316   
3317    '''
3318    if SGData['SpGrp'] == 'P 1':
3319        return [0,0,0,0]   
3320    hklmShape = Fhklm.shape
3321    hklmHalf = np.array(hklmShape)/2
3322    sortHKLM = np.argsort(Fhklm.flatten())
3323    Fdict = {}
3324    for hklm in sortHKLM:
3325        HKLM = np.unravel_index(hklm,hklmShape)
3326        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
3327        if F == 0.:
3328            break
3329        Fdict['%.6f'%(np.absolute(F))] = hklm
3330    Flist = np.flipud(np.sort(Fdict.keys()))
3331    F = str(1.e6)
3332    i = 0
3333    DH = []
3334    Dphi = []
3335    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3336    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3337    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3338    for F in Flist:
3339        hklm = np.unravel_index(Fdict[F],hklmShape)
3340        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
3341            continue
3342        Uniq = np.inner(hklm-hklmHalf,SSGMT)
3343        Phi = np.inner(hklm-hklmHalf,SSGT)
3344        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
3345        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3346        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
3347        ang0 = np.angle(Fh0,deg=True)/360.
3348        for H,phi in zip(Uniq,Phi)[1:]:
3349            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
3350            dH = H-hklm
3351            dang = ang-ang0
3352            DH.append(dH)
3353            Dphi.append((dang+.5) % 1.0)
3354        if i > 20 or len(DH) > 30:
3355            break
3356        i += 1
3357    DH = np.array(DH)
3358    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
3359    Dphi = np.array(Dphi)
3360    steps = np.array(hklmShape)
3361    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]]
3362    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
3363    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
3364    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3365    hist,bins = np.histogram(Mmap,bins=1000)
3366#    for i,item in enumerate(hist[:10]):
3367#        print item,bins[i]
3368    chisq = np.min(Mmap)
3369    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3370    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
3371#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3372    return DX
3373   
3374def SSChargeFlip(data,reflDict,pgbar):
3375    '''default doc string
3376   
3377    :param type name: description
3378   
3379    :returns: type name: description
3380   
3381    '''
3382    generalData = data['General']
3383    mapData = generalData['Map']
3384    map4DData = {}
3385    flipData = generalData['Flip']
3386    FFtable = {}
3387    if 'None' not in flipData['Norm element']:
3388        normElem = flipData['Norm element'].upper()
3389        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3390        for ff in FFs:
3391            if ff['Symbol'] == normElem:
3392                FFtable.update(ff)
3393    dmin = flipData['Resolution']
3394    SGData = generalData['SGData']
3395    SSGData = generalData['SSGData']
3396    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3397    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3398    cell = generalData['Cell'][1:8]       
3399    A = G2lat.cell2A(cell[:6])
3400    Vol = cell[6]
3401    maxM = 4
3402    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3403    adjHKLmax(SGData,Hmax)
3404    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3405    time0 = time.time()
3406    for iref,ref in enumerate(reflDict['RefList']):
3407        dsp = ref[5]
3408        if dsp > dmin:
3409            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3410            if FFtable:
3411                SQ = 0.25/dsp**2
3412                ff *= G2el.ScatFac(FFtable,SQ)[0]
3413            if ref[9] > 0.:         #use only +ve Fobs**2
3414                E = np.sqrt(ref[9])/ff
3415            else:
3416                E = 0.
3417            ph = ref[11]
3418            ph = rn.uniform(0.,360.)
3419            Uniq = np.inner(ref[:4],SSGMT)
3420            Phi = np.inner(ref[:4],SSGT)
3421            for i,hklm in enumerate(Uniq):        #uses uniq
3422                hklm = np.asarray(hklm,dtype='i')
3423                dp = 360.*Phi[i]                #and phi
3424                a = cosd(ph+dp)
3425                b = sind(ph+dp)
3426                phasep = complex(a,b)
3427                phasem = complex(a,-b)
3428                h,k,l,m = hklm+Hmax
3429                Ehkl[h,k,l,m] = E*phasep
3430                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
3431                Ehkl[h,k,l,m] = E*phasem
3432#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3433    CEhkl = copy.copy(Ehkl)
3434    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3435    Emask = ma.getmask(MEhkl)
3436    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3437    Ncyc = 0
3438    old = np.seterr(all='raise')
3439    while True:       
3440        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3441        CEsig = np.std(CErho)
3442        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3443        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3444        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3445        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3446        phase = CFhkl/np.absolute(CFhkl)
3447        CEhkl = np.absolute(Ehkl)*phase
3448        Ncyc += 1
3449        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3450        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3451        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3452        if Rcf < 5.:
3453            break
3454        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3455        if not GoOn or Ncyc > 10000:
3456            break
3457    np.seterr(**old)
3458    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3459    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
3460    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
3461    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3462    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3463
3464    mapData['Rcf'] = Rcf
3465    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3466    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3467    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3468    mapData['Type'] = reflDict['Type']
3469
3470    map4DData['Rcf'] = Rcf
3471    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))
3472    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3473    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3474    map4DData['Type'] = reflDict['Type']
3475    return mapData,map4DData
3476   
3477def getRho(xyz,mapData):
3478    ''' get scattering density at a point by 8-point interpolation
3479    param xyz: coordinate to be probed
3480    param: mapData: dict of map data
3481   
3482    :returns: density at xyz
3483    '''
3484    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3485    if not len(mapData):
3486        return 0.0
3487    rho = copy.copy(mapData['rho'])     #don't mess up original
3488    if not len(rho):
3489        return 0.0
3490    mapShape = np.array(rho.shape)
3491    mapStep = 1./mapShape
3492    X = np.array(xyz)%1.    #get into unit cell
3493    I = np.array(X*mapShape,dtype='int')
3494    D = X-I*mapStep         #position inside map cell
3495    D12 = D[0]*D[1]
3496    D13 = D[0]*D[2]
3497    D23 = D[1]*D[2]
3498    D123 = np.prod(D)
3499    Rho = rollMap(rho,-I)    #shifts map so point is in corner
3500    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]+  \
3501        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
3502        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
3503        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
3504    return R
3505       
3506def SearchMap(generalData,drawingData,Neg=False):
3507    '''Does a search of a density map for peaks meeting the criterion of peak
3508    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3509    mapData is data['General']['mapData']; the map is also in mapData.
3510
3511    :param generalData: the phase data structure; includes the map
3512    :param drawingData: the drawing data structure
3513    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
3514
3515    :returns: (peaks,mags,dzeros) where
3516
3517        * peaks : ndarray
3518          x,y,z positions of the peaks found in the map
3519        * mags : ndarray
3520          the magnitudes of the peaks
3521        * dzeros : ndarray
3522          the distance of the peaks from  the unit cell origin
3523        * dcent : ndarray
3524          the distance of the peaks from  the unit cell center
3525
3526    '''       
3527    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3528   
3529    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
3530   
3531#    def noDuplicate(xyz,peaks,Amat):
3532#        XYZ = np.inner(Amat,xyz)
3533#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3534#            print ' Peak',xyz,' <0.5A from another peak'
3535#            return False
3536#        return True
3537#                           
3538    def fixSpecialPos(xyz,SGData,Amat):
3539        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3540        X = []
3541        xyzs = [equiv[0] for equiv in equivs]
3542        for x in xyzs:
3543            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3544                X.append(x)
3545        if len(X) > 1:
3546            return np.average(X,axis=0)
3547        else:
3548            return xyz
3549       
3550    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3551        Mag,x0,y0,z0,sig = parms
3552        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
3553#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
3554        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
3555       
3556    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3557        Mag,x0,y0,z0,sig = parms
3558        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3559        return M
3560       
3561    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3562        Mag,x0,y0,z0,sig = parms
3563        dMdv = np.zeros(([5,]+list(rX.shape)))
3564        delt = .01
3565        for i in range(5):
3566            parms[i] -= delt
3567            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3568            parms[i] += 2.*delt
3569            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3570            parms[i] -= delt
3571            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3572        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3573        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3574        dMdv = np.reshape(dMdv,(5,rX.size))
3575        Hess = np.inner(dMdv,dMdv)
3576       
3577        return Vec,Hess
3578       
3579    SGData = generalData['SGData']
3580    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3581    peaks = []
3582    mags = []
3583    dzeros = []
3584    dcent = []
3585    try:
3586        mapData = generalData['Map']
3587        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3588        if Neg:
3589            rho = -copy.copy(mapData['rho'])     #flip +/-
3590        else:
3591            rho = copy.copy(mapData['rho'])     #don't mess up original
3592        mapHalf = np.array(rho.shape)/2
3593        res = mapData['Resolution']
3594        incre = np.array(rho.shape,dtype=np.float)
3595        step = max(1.0,1./res)+1
3596        steps = np.array((3*[step,]),dtype='int32')
3597    except KeyError:
3598        print '**** ERROR - Fourier map not defined'
3599        return peaks,mags
3600    rhoMask = ma.array(rho,mask=(rho<contLevel))
3601    indices = (-1,0,1)
3602    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3603    for roll in rolls:
3604        if np.any(roll):
3605            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3606    indx = np.transpose(rhoMask.nonzero())
3607    peaks = indx/incre
3608    mags = rhoMask[rhoMask.nonzero()]
3609    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3610        rho = rollMap(rho,ind)
3611        rMM = mapHalf-steps
3612        rMP = mapHalf+steps+1
3613        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3614        peakInt = np.sum(rhoPeak)*res**3
3615        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3616        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
3617        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3618            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3619        x1 = result[0]
3620        if not np.any(x1 < 0):
3621            peak = (np.array(x1[1:4])-ind)/incre
3622        peak = fixSpecialPos(peak,SGData,Amat)
3623        rho = rollMap(rho,-ind)
3624    cent = np.ones(3)*.5     
3625    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3626    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3627    if Neg:     #want negative magnitudes for negative peaks
3628        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3629    else:
3630        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3631   
3632def sortArray(data,pos,reverse=False):
3633    '''data is a list of items
3634    sort by pos in list; reverse if True
3635    '''
3636    T = []
3637    for i,M in enumerate(data):
3638        try:
3639            T.append((M[pos],i))
3640        except IndexError:
3641            return data
3642    D = dict(zip(T,data))
3643    T.sort()
3644    if reverse:
3645        T.reverse()
3646    X = []
3647    for key in T:
3648        X.append(D[key])
3649    return X
3650
3651def PeaksEquiv(data,Ind):
3652    '''Find the equivalent map peaks for those selected. Works on the
3653    contents of data['Map Peaks'].
3654
3655    :param data: the phase data structure
3656    :param list Ind: list of selected peak indices
3657    :returns: augmented list of peaks including those related by symmetry to the
3658      ones in Ind
3659
3660    '''       
3661    def Duplicate(xyz,peaks,Amat):
3662        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3663            return True
3664        return False
3665                           
3666    generalData = data['General']
3667    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3668    SGData = generalData['SGData']
3669    mapPeaks = data['Map Peaks']
3670    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3671    Indx = {}
3672    for ind in Ind:
3673        xyz = np.array(mapPeaks[ind][1:4])
3674        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3675        for jnd,xyz in enumerate(XYZ):       
3676            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3677    Ind = []
3678    for ind in Indx:
3679        if Indx[ind]:
3680            Ind.append(ind)
3681    return Ind
3682               
3683def PeaksUnique(data,Ind):
3684    '''Finds the symmetry unique set of peaks from those selected. Works on the
3685    contents of data['Map Peaks'].
3686
3687    :param data: the phase data structure
3688    :param list Ind: list of selected peak indices
3689    :returns: the list of symmetry unique peaks from among those given in Ind
3690
3691    '''       
3692#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3693
3694    def noDuplicate(xyz,peaks,Amat):
3695        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3696            return False
3697        return True
3698                           
3699    generalData = data['General']
3700    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3701    SGData = generalData['SGData']
3702    mapPeaks = data['Map Peaks']
3703    Indx = {}
3704    XYZ = {}
3705    for ind in Ind:
3706        XYZ[ind] = np.array(mapPeaks[ind][1:4])
3707        Indx[ind] = True
3708    for ind in Ind:
3709        if Indx[ind]:
3710            xyz = XYZ[ind]
3711            for jnd in Ind:
3712                if ind != jnd and Indx[jnd]:                       
3713                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3714                    xyzs = np.array([equiv[0] for equiv in Equiv])
3715                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3716    Ind = []
3717    for ind in Indx:
3718        if Indx[ind]:
3719            Ind.append(ind)
3720    return Ind
3721   
3722################################################################################
3723##### single peak fitting profile fxn stuff
3724################################################################################
3725
3726def getCWsig(ins,pos):
3727    '''get CW peak profile sigma^2
3728   
3729    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3730      as values only
3731    :param float pos: 2-theta of peak
3732    :returns: float getCWsig: peak sigma^2
3733   
3734    '''
3735    tp = tand(pos/2.0)
3736    return ins['U']*tp**2+ins['V']*tp+ins['W']
3737   
3738def getCWsigDeriv(pos):
3739    '''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3740   
3741    :param float pos: 2-theta of peak
3742   
3743    :returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3744   
3745    '''
3746    tp = tand(pos/2.0)
3747    return tp**2,tp,1.0
3748   
3749def getCWgam(ins,pos):
3750    '''get CW peak profile gamma
3751   
3752    :param dict ins: instrument parameters with at least 'X' & 'Y'
3753      as values only
3754    :param float pos: 2-theta of peak
3755    :returns: float getCWgam: peak gamma
3756   
3757    '''
3758    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
3759   
3760def getCWgamDeriv(pos):
3761    '''get derivatives of CW peak profile gamma wrt X & Y
3762   
3763    :param float pos: 2-theta of peak
3764   
3765    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3766   
3767    '''
3768    return 1./cosd(pos/2.0),tand(pos/2.0)
3769   
3770def getTOFsig(ins,dsp):
3771    '''get TOF peak profile sigma^2
3772   
3773    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3774      as values only
3775    :param float dsp: d-spacing of peak
3776   
3777    :returns: float getTOFsig: peak sigma^2
3778   
3779    '''
3780    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
3781   
3782def getTOFsigDeriv(dsp):
3783    '''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3784   
3785    :param float dsp: d-spacing of peak
3786   
3787    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3788   
3789    '''
3790    return 1.0,dsp**2,dsp**4,1./dsp**2
3791   
3792def getTOFgamma(ins,dsp):
3793    '''get TOF peak profile gamma
3794   
3795    :param dict ins: instrument parameters with at least 'X' & 'Y'
3796      as values only
3797    :param float dsp: d-spacing of peak
3798   
3799    :returns: float getTOFgamma: peak gamma
3800   
3801    '''
3802    return ins['X']*dsp+ins['Y']*dsp**2
3803   
3804def getTOFgammaDeriv(dsp):
3805    '''get derivatives of TOF peak profile gamma wrt X & Y
3806   
3807    :param float dsp: d-spacing of peak
3808   
3809    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3810   
3811    '''
3812    return dsp,dsp**2
3813   
3814def getTOFbeta(ins,dsp):
3815    '''get TOF peak profile beta
3816   
3817    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3818      as values only
3819    :param float dsp: d-spacing of peak
3820   
3821    :returns: float getTOFbeta: peak beat
3822   
3823    '''
3824    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
3825   
3826def getTOFbetaDeriv(dsp):
3827    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3828   
3829    :param float dsp: d-spacing of peak
3830   
3831    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3832   
3833    '''
3834    return 1.0,1./dsp**4,1./dsp**2
3835   
3836def getTOFalpha(ins,dsp):
3837    '''get TOF peak profile alpha
3838   
3839    :param dict ins: instrument parameters with at least 'alpha'
3840      as values only
3841    :param float dsp: d-spacing of peak
3842   
3843    :returns: flaot getTOFalpha: peak alpha
3844   
3845    '''
3846    return ins['alpha']/dsp
3847   
3848def getTOFalphaDeriv(dsp):
3849    '''get derivatives of TOF peak profile beta wrt alpha
3850   
3851    :param float dsp: d-spacing of peak
3852   
3853    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
3854   
3855    '''
3856    return 1./dsp
3857   
3858def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
3859    '''set starting peak parameters for single peak fits from plot selection or auto selection
3860   
3861    :param dict Parms: instrument parameters dictionary
3862    :param dict Parms2: table lookup for TOF profile coefficients
3863    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
3864    :param float mag: peak top magnitude from pick
3865    :param bool ifQ: True if pos in Q
3866    :param bool useFit: True if use fitted CW Parms values (not defaults)
3867   
3868    :returns: list XY: peak list entry:
3869        for CW: [pos,0,mag,1,sig,0,gam,0]
3870        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3871        NB: mag refinement set by default, all others off
3872   
3873    '''
3874    ind = 0
3875    if useFit:
3876        ind = 1
3877    ins = {}
3878    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
3879        for x in ['U','V','W','X','Y']:
3880            ins[x] = Parms[x][ind]
3881        if ifQ:                              #qplot - convert back to 2-theta
3882            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
3883        sig = getCWsig(ins,pos)
3884        gam = getCWgam(ins,pos)           
3885        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
3886    else:
3887        if ifQ:
3888            dsp = 2.*np.pi/pos
3889            pos = Parms['difC']*dsp
3890        else:
3891            dsp = pos/Parms['difC'][1]
3892        if 'Pdabc' in Parms2:
3893            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
3894                ins[x] = Parms[x][ind]
3895            Pdabc = Parms2['Pdabc'].T
3896            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
3897            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
3898        else:
3899            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
3900                ins[x] = Parms[x][ind]
3901            alp = getTOFalpha(ins,dsp)
3902            bet = getTOFbeta(ins,dsp)
3903        sig = getTOFsig(ins,dsp)
3904        gam = getTOFgamma(ins,dsp)
3905        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3906    return XY
3907   
3908################################################################################
3909##### MC/SA stuff
3910################################################################################
3911
3912#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
3913# Original Author: Travis Oliphant 2002
3914# Bug-fixes in 2006 by Tim Leslie
3915
3916
3917import numpy
3918from numpy import asarray, exp, squeeze, sign, \
3919     all, shape, array, where
3920from numpy import random
3921
3922#__all__ = ['anneal']
3923
3924_double_min = numpy.finfo(float).min
3925_double_max = numpy.finfo(float).max
3926class base_schedule(object):
3927    def __init__(self):
3928        self.dwell = 20
3929        self.lower = 0.
3930        self.upper = 1.
3931        self.Ninit = 50
3932        self.accepted = 0
3933        self.tests = 0
3934        self.feval = 0
3935        self.k = 0
3936        self.T = None
3937
3938    def init(self, **options):
3939        self.__dict__.update(options)
3940        self.lower = asarray(self.lower)
3941        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
3942        self.upper = asarray(self.upper)
3943        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
3944        self.k = 0
3945        self.accepted = 0
3946        self.feval = 0
3947        self.tests = 0
3948
3949    def getstart_temp(self, best_state):
3950        """ Find a matching starting temperature and starting parameters vector
3951        i.e. find x0 such that func(x0) = T0.
3952
3953        :param best_state: _state
3954            A _state object to store the function value and x0 found.
3955
3956        :returns: x0 : array
3957            The starting parameters vector.
3958        """
3959
3960        assert(not self.dims is None)
3961        lrange = self.lower
3962        urange = self.upper
3963        fmax = _double_min
3964        fmin = _double_max
3965        for _ in range(self.Ninit):
3966            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
3967            fval = self.func(x0, *self.args)
3968            self.feval += 1
3969            if fval > fmax:
3970                fmax = fval
3971            if fval < fmin:
3972                fmin = fval
3973                best_state.cost = fval
3974                best_state.x = array(x0)
3975
3976        self.T0 = (fmax-fmin)*1.5
3977        return best_state.x
3978       
3979    def set_range(self,x0,frac):
3980        delrange = frac*(self.upper-self.lower)
3981        self.upper = x0+delrange
3982        self.lower = x0-delrange
3983
3984    def accept_test(self, dE):
3985        T = self.T
3986        self.tests += 1
3987        if dE < 0:
3988            self.accepted += 1
3989            return 1
3990        p = exp(-dE*1.0/T)
3991        if (p > random.uniform(0.0, 1.0)):
3992            self.accepted += 1
3993            return 1
3994        return 0
3995
3996    def update_guess(self, x0):
3997        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3998
3999    def update_temp(self, x0):
4000        pass
4001
4002class fast_sa(base_schedule):
4003    def init(self, **options):
4004        self.__dict__.update(options)
4005
4006    def update_guess(self, x0):
4007        x0 = asarray(x0)
4008        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4009        T = self.T
4010        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4011        xnew = xc*(self.upper - self.lower)+self.lower
4012        return xnew
4013
4014    def update_temp(self):
4015        self.T = self.T0*exp(-self.c * self.k**(self.quench))
4016        self.k += 1
4017        return
4018
4019class log_sa(base_schedule):        #OK
4020
4021    def init(self,**options):
4022        self.__dict__.update(options)
4023       
4024    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
4025        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4026        T = self.T
4027        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4028        xnew = xc*(self.upper - self.lower)+self.lower
4029        return xnew
4030       
4031    def update_temp(self):
4032        self.k += 1
4033        self.T = self.T0*self.slope**self.k
4034       
4035class _state(object):
4036    def __init__(self):
4037        self.x = None
4038        self.cost = None
4039
4040def makeTsched(data):
4041    if data['Algorithm'] == 'fast':
4042        sched = fast_sa()
4043        sched.quench = data['fast parms'][0]
4044        sched.c = data['fast parms'][1]
4045    elif data['Algorithm'] == 'log':
4046        sched = log_sa()
4047        sched.slope = data['log slope']
4048    sched.T0 = data['Annealing'][0]
4049    if not sched.T0:
4050        sched.T0 = 50.
4051    Tf = data['Annealing'][1]
4052    if not Tf:
4053        Tf = 0.001
4054    Tsched = [sched.T0,]
4055    while Tsched[-1] > Tf:
4056        sched.update_temp()
4057        Tsched.append(sched.T)
4058    return Tsched[1:]
4059   
4060def anneal(func, x0, args=(), schedule='fast', 
4061           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4062           feps=1e-6, quench=1.0, c=1.0,
4063           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4064           ranRange=0.10,autoRan=False,dlg=None):
4065    ''' Scipy license:
4066        Copyright (c) 2001, 2002 Enthought, Inc.
4067    All rights reserved.
4068   
4069    Copyright (c) 2003-2016 SciPy Developers.
4070    All rights reserved.
4071   
4072    Redistribution and use in source and binary forms, with or without
4073    modification, are permitted provided that the following conditions are met:
4074   
4075      a. Redistributions of source code must retain the above copyright notice,
4076         this list of conditions and the following disclaimer.
4077      b. Redistributions in binary form must reproduce the above copyright
4078         notice, this list of conditions and the following disclaimer in the
4079         documentation and/or other materials provided with the distribution.
4080      c. Neither the name of Enthought nor the names of the SciPy Developers
4081         may be used to endorse or promote products derived from this software
4082         without specific prior written permission.
4083    '''
4084
4085    """Minimize a function using simulated annealing.
4086
4087    Schedule is a schedule class implementing the annealing schedule.
4088    Available ones are 'fast', 'cauchy', 'boltzmann'
4089
4090    :param callable func: f(x, \*args)
4091        Function to be optimized.
4092    :param ndarray x0:
4093        Initial guess.
4094    :param tuple args:
4095        Extra parameters to `func`.
4096    :param base_schedule schedule:
4097        Annealing schedule to use (a class).
4098    :param float T0:
4099        Initial Temperature (estimated as 1.2 times the largest
4100        cost-function deviation over random points in the range).
4101    :param float Tf:
4102        Final goal temperature.
4103    :param int maxeval:
4104        Maximum function evaluations.
4105    :param int maxaccept:
4106        Maximum changes to accept.
4107    :param int maxiter:
4108        Maximum cooling iterations.
4109    :param float feps:
4110        Stopping relative error tolerance for the function value in
4111        last four coolings.
4112    :param float quench,c:
4113        Parameters to alter fast_sa schedule.
4114    :param float/ndarray lower,upper:
4115        Lower and upper bounds on `x`.
4116    :param int dwell:
4117        The number of times to search the space at each temperature.
4118    :param float slope:
4119        Parameter for log schedule
4120    :param bool ranStart=False:
4121        True for set 10% of ranges about x
4122
4123    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4124
4125     * xmin (ndarray): Point giving smallest value found.
4126     * Jmin (float): Minimum value of function found.
4127     * T (float): Final temperature.
4128     * feval (int): Number of function evaluations.
4129     * iters (int): Number of cooling iterations.
4130     * accept (int): Number of tests accepted.
4131     * retval (int): Flag indicating stopping condition:
4132
4133              *  0: Points no longer changing
4134              *  1: Cooled to final temperature
4135              *  2: Maximum function evaluations
4136              *  3: Maximum cooling iterations reached
4137              *  4: Maximum accepted query locations reached
4138              *  5: Final point not the minimum amongst encountered points
4139
4140    *Notes*:
4141    Simulated annealing is a random algorithm which uses no derivative
4142    information from the function being optimized. In practice it has
4143    been more useful in discrete optimization than continuous
4144    optimization, as there are usually better algorithms for continuous
4145    optimization problems.
4146
4147    Some experimentation by trying the difference temperature
4148    schedules and altering their parameters is likely required to
4149    obtain good performance.
4150
4151    The randomness in the algorithm comes from random sampling in numpy.
4152    To obtain the same results you can call numpy.random.seed with the
4153    same seed immediately before calling scipy.optimize.anneal.
4154
4155    We give a brief description of how the three temperature schedules
4156    generate new points and vary their temperature. Temperatures are
4157    only updated with iterations in the outer loop. The inner loop is
4158    over xrange(dwell), and new points are generated for
4159    every iteration in the inner loop. (Though whether the proposed
4160    new points are accepted is probabilistic.)
4161
4162    For readability, let d denote the dimension of the inputs to func.
4163    Also, let x_old denote the previous state, and k denote the
4164    iteration number of the outer loop. All other variables not
4165    defined below are input variables to scipy.optimize.anneal itself.
4166
4167    In the 'fast' schedule the updates are ::
4168
4169        u ~ Uniform(0, 1, size=d)
4170        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4171        xc = y * (upper - lower)
4172        x_new = x_old + xc
4173
4174        T_new = T0 * exp(-c * k**quench)
4175
4176
4177    """
4178    x0 = asarray(x0)
4179    lower = asarray(lower)
4180    upper = asarray(upper)
4181
4182    schedule = eval(schedule+'_sa()')
4183    #   initialize the schedule
4184    schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4185        c=c, quench=quench, dwell=dwell, slope=slope)
4186
4187    current_state, last_state, best_state = _state(), _state(), _state()
4188    if ranStart:
4189        schedule.set_range(x0,ranRange)
4190    if T0 is None:
4191        x0 = schedule.getstart_temp(best_state)
4192    else:
4193        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4194        best_state.x = None
4195        best_state.cost = numpy.Inf
4196
4197    last_state.x = asarray(x0).copy()
4198    fval = func(x0,*args)
4199    schedule.feval += 1
4200    last_state.cost = fval
4201    if last_state.cost < best_state.cost:
4202        best_state.cost = fval
4203        best_state.x = asarray(x0).copy()
4204    schedule.T = schedule.T0
4205    fqueue = [100, 300, 500, 700]
4206    iters = 1
4207    keepGoing = True
4208    bestn = 0
4209    while keepGoing:
4210        retval = 0
4211        for n in xrange(dwell):
4212            current_state.x = schedule.update_guess(last_state.x)
4213            current_state.cost = func(current_state.x,*args)
4214            schedule.feval += 1
4215
4216            dE = current_state.cost - last_state.cost
4217            if schedule.accept_test(dE):
4218                last_state.x = current_state.x.copy()
4219                last_state.cost = current_state.cost
4220                if last_state.cost < best_state.cost:
4221                    best_state.x = last_state.x.copy()
4222                    best_state.cost = last_state.cost
4223                    bestn = n
4224                    if best_state.cost < 1.0 and autoRan:
4225                        schedule.set_range(x0,best_state.cost/2.)                       
4226        if dlg:
4227            GoOn = dlg.Update(min(100.,best_state.cost*100),
4228                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4229                    'Best trial:',bestn,  \
4230                    'MC/SA Residual:',best_state.cost*100,'%', \
4231                    ))[0]
4232            if not GoOn:
4233                best_state.x = last_state.x.copy()
4234                best_state.cost = last_state.cost
4235                retval = 5
4236        schedule.update_temp()
4237        iters += 1
4238        # Stopping conditions
4239        # 0) last saved values of f from each cooling step
4240        #     are all very similar (effectively cooled)
4241        # 1) Tf is set and we are below it
4242        # 2) maxeval is set and we are past it
4243        # 3) maxiter is set and we are past it
4244        # 4) maxaccept is set and we are past it
4245        # 5) user canceled run via progress bar
4246
4247        fqueue.append(squeeze(last_state.cost))
4248        fqueue.pop(0)
4249        af = asarray(fqueue)*1.0
4250        if retval == 5:
4251            print ' User terminated run; incomplete MC/SA'
4252            keepGoing = False
4253            break
4254        if all(abs((af-af[0])/af[0]) < feps):
4255            retval = 0
4256            if abs(af[-1]-best_state.cost) > feps*10:
4257                retval = 5
4258                print " Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1)
4259            break
4260        if (Tf is not None) and (schedule.T < Tf):
4261#            print ' Minimum T reached in %d steps'%(iters-1)
4262            retval = 1
4263            break
4264        if (maxeval is not None) and (schedule.feval > maxeval):
4265            retval = 2
4266            break
4267        if (iters > maxiter):
4268            print  " Warning: Maximum number of iterations exceeded."
4269            retval = 3
4270            break
4271        if (maxaccept is not None) and (schedule.accepted > maxaccept):
4272            retval = 4
4273            break
4274
4275    return best_state.x, best_state.cost, schedule.T, \
4276           schedule.feval, iters, schedule.accepted, retval
4277
4278def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4279    outlist = []
4280    timelist = []
4281    nsflist = []
4282    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
4283    for n in range(iCyc):
4284        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False)         #mcsa result,time,rcov
4285        outlist.append(result[0])
4286        timelist.append(result[1])
4287        nsflist.append(result[2])
4288        print ' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
4289    out_q.put(outlist)
4290    out_t.put(timelist)
4291    out_n.put(nsflist)
4292
4293def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4294    import multiprocessing as mp
4295   
4296    out_q = mp.Queue()
4297    out_t = mp.Queue()
4298    out_n = mp.Queue()
4299    procs = []
4300    totsftime = 0.
4301    totnsf = 0
4302    iCyc = np.zeros(nprocs)
4303    for i in range(nCyc):
4304        iCyc[i%nprocs] += 1
4305    for i in range(nprocs):
4306        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4307        procs.append(p)
4308        p.start()
4309    resultlist = []
4310    for i in range(nprocs):
4311        resultlist += out_q.get()
4312        totsftime += np.sum(out_t.get())
4313        totnsf += np.sum(out_n.get())
4314    for p in procs:
4315        p.join()
4316    return resultlist,totsftime,totnsf
4317
4318def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4319    '''default doc string
4320   
4321    :param type name: description
4322   
4323    :returns: type name: description
4324    '''
4325   
4326    class RandomDisplacementBounds(object):
4327        """random displacement with bounds"""
4328        def __init__(self, xmin, xmax, stepsize=0.5):
4329            self.xmin = xmin
4330            self.xmax = xmax
4331            self.stepsize = stepsize
4332   
4333        def __call__(self, x):
4334            """take a random step but ensure the new position is within the bounds"""
4335            while True:
4336                # this could be done in a much more clever way, but it will work for example purposes
4337                steps = self.xmax-self.xmin
4338                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
4339                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4340                    break
4341            return xnew
4342   
4343    global tsum,nsum
4344    tsum = 0.
4345    nsum = 0
4346   
4347    def getMDparms(item,pfx,parmDict,varyList):
4348        parmDict[pfx+'MDaxis'] = item['axis']
4349        parmDict[pfx+'MDval'] = item['Coef'][0]
4350        if item['Coef'][1]:
4351            varyList += [pfx+'MDval',]
4352            limits = item['Coef'][2]
4353            lower.append(limits[0])
4354            upper.append(limits[1])
4355                       
4356    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4357        parmDict[pfx+'Atype'] = item['atType']
4358        aTypes |= set([item['atType'],]) 
4359        pstr = ['Ax','Ay','Az']
4360        XYZ = [0,0,0]
4361        for i in range(3):
4362            name = pfx+pstr[i]
4363            parmDict[name] = item['Pos'][0][i]
4364            XYZ[i] = parmDict[name]
4365            if item['Pos'][1][i]:
4366                varyList += [name,]
4367                limits = item['Pos'][2][i]
4368                lower.append(limits[0])
4369                upper.append(limits[1])
4370        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4371           
4372    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4373        parmDict[mfx+'MolCent'] = item['MolCent']
4374        parmDict[mfx+'RBId'] = item['RBId']
4375        pstr = ['Px','Py','Pz']
4376        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4377        for i in range(3):
4378            name = mfx+pstr[i]
4379            parmDict[name] = item['Pos'][0][i]
4380            if item['Pos'][1][i]:
4381                varyList += [name,]
4382                limits = item['Pos'][2][i]
4383                lower.append(limits[0])
4384                upper.append(limits[1])
4385        AV = item['Ori'][0]
4386        A = AV[0]
4387        V = AV[1:]
4388        for i in range(4):
4389            name = mfx+ostr[i]
4390            if i:
4391                parmDict[name] = V[i-1]
4392            else:
4393                parmDict[name] = A
4394            if item['Ovar'] == 'AV':
4395                varyList += [name,]
4396                limits = item['Ori'][2][i]
4397                lower.append(limits[0])
4398                upper.append(limits[1])
4399            elif item['Ovar'] == 'A' and not i:
4400                varyList += [name,]
4401                limits = item['Ori'][2][i]
4402                lower.append(limits[0])
4403                upper.append(limits[1])
4404        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
4405            for i in range(len(item['Tor'][0])):
4406                name = mfx+'Tor'+str(i)
4407                parmDict[name] = item['Tor'][0][i]
4408                if item['Tor'][1][i]:
4409                    varyList += [name,]
4410                    limits = item['Tor'][2][i]
4411                    lower.append(limits[0])
4412                    upper.append(limits[1])
4413        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
4414        aTypes |= set(atypes)
4415        atNo += len(atypes)
4416        return atNo
4417               
4418    def GetAtomM(Xdata,SGData):
4419        Mdata = []
4420        for xyz in Xdata:
4421            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
4422        return np.array(Mdata)
4423       
4424    def GetAtomT(RBdata,parmDict):
4425        'Needs a doc string'
4426        atNo = parmDict['atNo']
4427        nfixAt = parmDict['nfixAt']
4428        Tdata = atNo*[' ',]
4429        for iatm in range(nfixAt):
4430            parm = ':'+str(iatm)+':Atype'
4431            if parm in parmDict:
4432                Tdata[iatm] = aTypes.index(parmDict[parm])
4433        iatm = nfixAt
4434        for iObj in range(parmDict['nObj']):
4435            pfx = str(iObj)+':'
4436            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4437                if parmDict[pfx+'Type'] == 'Vector':
4438                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4439                    nAtm = len(RBRes['rbVect'][0])
4440                else:       #Residue
4441                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4442                    nAtm = len(RBRes['rbXYZ'])
4443                for i in range(nAtm):
4444                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
4445                    iatm += 1
4446            elif parmDict[pfx+'Type'] == 'Atom':
4447                atNo = parmDict[pfx+'atNo']
4448                parm = pfx+'Atype'              #remove extra ':'
4449                if parm in parmDict:
4450                    Tdata[atNo] = aTypes.index(parmDict[parm])
4451                iatm += 1
4452            else:
4453                continue        #skips March Dollase
4454        return Tdata
4455       
4456    def GetAtomX(RBdata,parmDict):
4457        'Needs a doc string'
4458        Bmat = parmDict['Bmat']
4459        atNo = parmDict['atNo']
4460        nfixAt = parmDict['nfixAt']
4461        Xdata = np.zeros((3,atNo))
4462        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
4463        for iatm in range(nfixAt):
4464            for key in keys:
4465                parm = ':'+str(iatm)+key
4466                if parm in parmDict:
4467                    keys[key][iatm] = parmDict[parm]
4468        iatm = nfixAt
4469        for iObj in range(parmDict['nObj']):
4470            pfx = str(iObj)+':'
4471            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4472                if parmDict[pfx+'Type'] == 'Vector':
4473                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4474                    vecs = RBRes['rbVect']
4475                    mags = RBRes['VectMag']
4476                    Cart = np.zeros_like(vecs[0])
4477                    for vec,mag in zip(vecs,mags):
4478                        Cart += vec*mag
4479                elif parmDict[pfx+'Type'] == 'Residue':
4480                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4481                    Cart = np.array(RBRes['rbXYZ'])
4482                    for itor,seq in enumerate(RBRes['rbSeq']):
4483                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
4484                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
4485                if parmDict[pfx+'MolCent'][1]:
4486                    Cart -= parmDict[pfx+'MolCent'][0]
4487                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4488                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
4489                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
4490                iatm += len(Cart)
4491            elif parmDict[pfx+'Type'] == 'Atom':
4492                atNo = parmDict[pfx+'atNo']
4493                for key in keys:
4494                    parm = pfx+key[1:]              #remove extra ':'
4495                    if parm in parmDict:
4496                        keys[key][atNo] = parmDict[parm]
4497                iatm += 1
4498            else:
4499                continue        #skips March Dollase
4500        return Xdata.T
4501       
4502    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
4503        allX = np.inner(Xdata,SGM)+SGT
4504        allT = np.repeat(Tdata,allX.shape[1])
4505        allM = np.repeat(Mdata,allX.shape[1])
4506        allX = np.reshape(allX,(-1,3))
4507        return allT,allM,allX
4508
4509    def getAllX(Xdata,SGM,SGT):
4510        allX = np.inner(Xdata,SGM)+SGT
4511        allX = np.reshape(allX,(-1,3))
4512        return allX
4513       
4514    def normQuaternions(RBdata,parmDict,varyList,values):
4515        for iObj in range(parmDict['nObj']):
4516            pfx = str(iObj)+':'
4517            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4518                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4519                A,V = Q2AVdeg(Qori)
4520                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
4521                    if i:
4522                        parmDict[pfx+name] = V[i-1]
4523                    else:
4524                        parmDict[pfx+name] = A
4525       
4526    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
4527        ''' Compute structure factors for all h,k,l for phase