source: trunk/GSASIImath.py @ 3877

Last change on this file since 3877 was 3877, checked in by vondreele, 3 years ago

micsf again

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