source: trunk/GSASIImath.py @ 3738

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

fix mag form factor table so Ce & Ce+2 are both seen. Comment section moved to end of dict.
fix export powder data as csv & text - python 3 vs 2 problem
fix p 62'2' operator problem - reverse order of generators
fixes for some incommensurate mag problems - still more to be found

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