source: trunk/GSASIImath.py @ 3287

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

fix atom transformation issues in modulation drawings.

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