source: trunk/GSASIImath.py @ 3219

Last change on this file since 3219 was 3219, checked in by vondreele, 5 years ago

fix import mag commensurate structures - all drawn correctly now
modulated mags not correct yet

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