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

Last change on this file since 2957 was 2957, checked in by toby, 5 years ago

use grid scroll bars for Peak List and Index Peak List displays; increase convergence on SVD lam increase

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