source: trunk/GSASIImath.py @ 2860

Last change on this file since 2860 was 2860, checked in by vondreele, 6 years ago

a bit more of protein validation

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