source: trunk/GSASIImath.py @ 3568

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

finish interface with k-SUBGROUPSMAG output - selection of space group & transformation to make magnetic phase

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 189.5 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2018-08-23 20:40:51 +0000 (Thu, 23 Aug 2018) $
5# $Author: vondreele $
6# $Revision: 3568 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 3568 2018-08-23 20:40:51Z 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: 3568 $")
27import GSASIIElem as G2el
28import GSASIIlattice as G2lat
29import GSASIIspc as G2spc
30import GSASIIpwd as G2pwd
31import numpy.fft as fft
32import scipy.optimize as so
33try:
34    import pypowder as pwd
35except ImportError:
36    print ('pypowder is not available - profile calcs. not allowed')
37
38sind = lambda x: np.sin(x*np.pi/180.)
39cosd = lambda x: np.cos(x*np.pi/180.)
40tand = lambda x: np.tan(x*np.pi/180.)
41asind = lambda x: 180.*np.arcsin(x)/np.pi
42acosd = lambda x: 180.*np.arccos(x)/np.pi
43atand = lambda x: 180.*np.arctan(x)/np.pi
44atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
45twopi = 2.0*np.pi
46twopisq = 2.0*np.pi**2
47nxs = np.newaxis
48   
49################################################################################
50##### Hessian least-squares Levenberg-Marquardt routine
51################################################################################
52
53def pinv(a, rcond=1e-15 ):
54    """
55    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
56    Modified from numpy.linalg.pinv; assumes a is Hessian & returns no. zeros found
57    Calculate the generalized inverse of a matrix using its
58    singular-value decomposition (SVD) and including all
59    *large* singular values.
60
61    :param array a: (M, M) array_like - here assumed to be LS Hessian
62      Matrix to be pseudo-inverted.
63    :param float rcond: Cutoff for small singular values.
64      Singular values smaller (in modulus) than
65      `rcond` * largest_singular_value (again, in modulus)
66      are set to zero.
67
68    :returns: B : (M, M) ndarray
69      The pseudo-inverse of `a`
70
71    Raises: LinAlgError
72      If the SVD computation does not converge.
73
74    Notes:
75      The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
76      defined as: "the matrix that 'solves' [the least-squares problem]
77      :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
78      :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
79
80    It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
81    value decomposition of A, then
82    :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
83    orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
84    of A's so-called singular values, (followed, typically, by
85    zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
86    consisting of the reciprocals of A's singular values
87    (again, followed by zeros). [1]
88
89    References:
90    .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL, Academic Press, Inc., 1980, pp. 139-142.
91
92    """
93    u, s, vt = nl.svd(a, 0)
94    cutoff = rcond*np.maximum.reduce(s)
95    s = np.where(s>cutoff,1./s,0.)
96    nzero = s.shape[0]-np.count_nonzero(s)
97#    res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))
98    res = np.dot(vt.T,s[:,nxs]*u.T)
99    return res,nzero
100
101def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
102   
103    """
104    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
105    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`   
106    where :math:`x = arg min(\sum_{y=0}^{N_{obs}} (func(y)^2,axis=0))`
107
108    :param function func: callable method or function
109        should take at least one (possibly length N vector) argument and
110        returns M floating point numbers.
111    :param np.ndarray x0: The starting estimate for the minimization of length N
112    :param function Hess: callable method or function
113        A required function or method to compute the weighted vector and Hessian for func.
114        It must be a symmetric NxN array
115    :param tuple args: Any extra arguments to func are placed in this tuple.
116    :param float ftol: Relative error desired in the sum of squares.
117    :param float xtol: Relative tolerance of zeros in the SVD solution in nl.pinv.
118    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
119        until other limits are met (ftol, xtol)
120    :param int lamda: initial Marquardt lambda=10**lamda
121    :param bool Print: True for printing results (residuals & times) by cycle
122
123    :returns: (x,cov_x,infodict) where
124
125      * x : ndarray
126        The solution (or the result of the last iteration for an unsuccessful
127        call).
128      * cov_x : ndarray
129        Uses the fjac and ipvt optional outputs to construct an
130        estimate of the jacobian around the solution.  ``None`` if a
131        singular matrix encountered (indicates very flat curvature in
132        some direction).  This matrix must be multiplied by the
133        residual standard deviation to get the covariance of the
134        parameter estimates -- see curve_fit.
135      * infodict : dict
136        a dictionary of optional outputs with the keys:
137
138         * 'fvec' : the function evaluated at the output
139         * 'num cyc':
140         * 'nfev':
141         * 'lamMax':
142         * 'psing':
143         * 'SVD0':
144           
145    """
146               
147    ifConverged = False
148    deltaChi2 = -10.
149    x0 = np.array(x0, ndmin=1)      #might be redundant?
150    n = len(x0)
151    if type(args) != type(()):
152        args = (args,)
153       
154    icycle = 0
155    One = np.ones((n,n))
156    lam = 10.**lamda
157    lamMax = lam
158    nfev = 0
159    if Print:
160        print (' Hessian Levenburg-Marquardt SVD refinement on %d variables:'%(n))
161    Lam = np.zeros((n,n))
162    while icycle < maxcyc:
163        time0 = time.time()
164        M = func(x0,*args)
165        Nobs = len(M)
166        nfev += 1
167        chisq0 = np.sum(M**2)
168        Yvec,Amat = Hess(x0,*args)
169        Adiag = np.sqrt(np.diag(Amat))
170        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
171        if np.any(psing):                #hard singularity in matrix
172            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
173        Anorm = np.outer(Adiag,Adiag)
174        Yvec /= Adiag
175        Amat /= Anorm
176        if Print:
177            print ('initial chi^2 %.5g on %d obs.'%(chisq0,Nobs))
178        chitol = ftol
179        while True:
180            Lam = np.eye(Amat.shape[0])*lam
181            Amatlam = Amat*(One+Lam)
182            try:
183                Ainv,Nzeros = pinv(Amatlam,xtol)    #do Moore-Penrose inversion (via SVD)
184            except nl.LinAlgError:
185                print ('ouch #1 bad SVD inversion; change parameterization')
186                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
187                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
188            Xvec = np.inner(Ainv,Yvec)      #solve
189            Xvec /= Adiag
190            M2 = func(x0+Xvec,*args)
191            nfev += 1
192            chisq1 = np.sum(M2**2)
193            if chisq1 > chisq0*(1.+chitol):     #TODO put Alan Coehlo's criteria for lambda here?
194                lam *= 10.
195                if Print:
196                    print ('new chi^2 %.5g on %d obs., %d SVD zeros ; matrix modification needed; lambda now %.1e'  \
197                           %(chisq1,Nobs,Nzeros,lam))
198            else:
199                x0 += Xvec
200                lam /= 10.
201                break
202            if lam > 10.:
203                print ('ouch #3 chisq1 %g.4 stuck > chisq0 %g.4'%(chisq1,chisq0))
204                break
205            chitol *= 2
206        lamMax = max(lamMax,lam)
207        deltaChi2 = (chisq0-chisq1)/chisq0
208        if Print:
209            print (' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g'%(
210                icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
211        if deltaChi2 < ftol:
212            ifConverged = True
213            if Print: print ("converged")
214            break
215        icycle += 1
216    M = func(x0,*args)
217    nfev += 1
218    Yvec,Amat = Hess(x0,*args)
219    Adiag = np.sqrt(np.diag(Amat))
220    Anorm = np.outer(Adiag,Adiag)
221    Lam = np.eye(Amat.shape[0])*lam
222    Amatlam = Amat/Anorm       
223    try:
224        Bmat,Nzero = pinv(Amatlam,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
225        if Print: 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):               #ToDo: something odd here, but it works
1703                    if icent < 0:                       
1704                        mom += np.sum(posFourier(.25-tauT,np.array(scof),np.array(ccof)),axis=1)
1705                    else:
1706                        mom += np.sum(posFourier(tauT*icent,np.array(scof),np.array(ccof)),axis=1)
1707            if len(Sadp):
1708                scof = []
1709                ccof = []
1710                waveType = Sadp[0]
1711                for i,sadp in enumerate(Sadp[1:]):
1712                    scof.append(sadp[0][:6])
1713                    ccof.append(sadp[0][6:])
1714                uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1715            if atom[cia] == 'A':                   
1716                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
1717                drawatom[dcx:dcx+3] = X
1718                drawatom[dci-6:dci] = U
1719            else:
1720                X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
1721                drawatom[dcx:dcx+3] = X
1722            if generalData['Type'] == 'magnetic':
1723                M = G2spc.ApplyStringOpsMom(opr,SGData,atmom+mom)
1724                drawatom[dcx+3:dcx+6] = M
1725    return drawAtoms,Fade
1726   
1727# gauleg.py Gauss Legendre numerical quadrature, x and w computation
1728# integrate from a to b using n evaluations of the function f(x) 
1729# usage: from gauleg import gaulegf         
1730#        x,w = gaulegf( a, b, n)                               
1731#        area = 0.0                                           
1732#        for i in range(1,n+1):          #  yes, 1..n                   
1733#          area += w[i]*f(x[i])                                   
1734
1735def gaulegf(a, b, n):
1736    x = range(n+1) # x[0] unused
1737    w = range(n+1) # w[0] unused
1738    eps = 3.0E-14
1739    m = (n+1)/2
1740    xm = 0.5*(b+a)
1741    xl = 0.5*(b-a)
1742    for i in range(1,m+1):
1743        z = math.cos(3.141592654*(i-0.25)/(n+0.5))
1744        while True:
1745            p1 = 1.0
1746            p2 = 0.0
1747            for j in range(1,n+1):
1748                p3 = p2
1749                p2 = p1
1750                p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
1751       
1752            pp = n*(z*p1-p2)/(z*z-1.0)
1753            z1 = z
1754            z = z1 - p1/pp
1755            if abs(z-z1) <= eps:
1756                break
1757
1758        x[i] = xm - xl*z
1759        x[n+1-i] = xm + xl*z
1760        w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
1761        w[n+1-i] = w[i]
1762    return np.array(x), np.array(w)
1763# end gaulegf
1764   
1765   
1766def BessJn(nmax,x):
1767    ''' compute Bessel function J(n,x) from scipy routine & recurrance relation
1768    returns sequence of J(n,x) for n in range [-nmax...0...nmax]
1769   
1770    :param  integer nmax: maximul order for Jn(x)
1771    :param  float x: argument for Jn(x)
1772   
1773    :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
1774   
1775    '''
1776    import scipy.special as sp
1777    bessJn = np.zeros(2*nmax+1)
1778    bessJn[nmax] = sp.j0(x)
1779    bessJn[nmax+1] = sp.j1(x)
1780    bessJn[nmax-1] = -bessJn[nmax+1]
1781    for i in range(2,nmax+1):
1782        bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
1783        bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i
1784    return bessJn
1785   
1786def BessIn(nmax,x):
1787    ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
1788    returns sequence of I(n,x) for n in range [-nmax...0...nmax]
1789   
1790    :param  integer nmax: maximul order for In(x)
1791    :param  float x: argument for In(x)
1792   
1793    :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
1794   
1795    '''
1796    import scipy.special as sp
1797    bessIn = np.zeros(2*nmax+1)
1798    bessIn[nmax] = sp.i0(x)
1799    bessIn[nmax+1] = sp.i1(x)
1800    bessIn[nmax-1] = bessIn[nmax+1]
1801    for i in range(2,nmax+1):
1802        bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x
1803        bessIn[nmax-i] = bessIn[i+nmax]
1804    return bessIn
1805       
1806   
1807################################################################################
1808##### distance, angle, planes, torsion stuff
1809################################################################################
1810
1811def CalcDist(distance_dict, distance_atoms, parmDict):
1812    if not len(parmDict):
1813        return 0.
1814    pId = distance_dict['pId']
1815    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1816    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1817    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1818    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1819    inv = 1
1820    symNo = distance_dict['symNo']
1821    if symNo < 0:
1822        inv = -1
1823        symNo *= -1
1824    cen = symNo//100
1825    op = symNo%100-1
1826    M,T = distance_dict['SGData']['SGOps'][op]
1827    D = T*inv+distance_dict['SGData']['SGCen'][cen]
1828    D += distance_dict['cellNo']
1829    Txyz = np.inner(M*inv,Txyz)+D
1830    dist = np.sqrt(np.sum(np.inner(Amat,(Txyz-Oxyz))**2))
1831#    GSASIIpath.IPyBreak()
1832    return dist   
1833   
1834def CalcDistDeriv(distance_dict, distance_atoms, parmDict):
1835    if not len(parmDict):
1836        return None
1837    pId = distance_dict['pId']
1838    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1839    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1840    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1841    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1842    symNo = distance_dict['symNo']
1843    Tunit = distance_dict['cellNo']
1844    SGData = distance_dict['SGData']   
1845    deriv = getDistDerv(Oxyz,Txyz,Amat,Tunit,symNo,SGData)
1846    return deriv
1847   
1848def CalcAngle(angle_dict, angle_atoms, parmDict):
1849    if not len(parmDict):
1850        return 0.
1851    pId = angle_dict['pId']
1852    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1853    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1854    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1855    Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1856    Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1857    ABxyz = [Axyz,Bxyz]
1858    symNo = angle_dict['symNo']
1859    vec = np.zeros((2,3))
1860    for i in range(2):
1861        inv = 1
1862        if symNo[i] < 0:
1863            inv = -1
1864        cen = inv*symNo[i]//100
1865        op = inv*symNo[i]%100-1
1866        M,T = angle_dict['SGData']['SGOps'][op]
1867        D = T*inv+angle_dict['SGData']['SGCen'][cen]
1868        D += angle_dict['cellNo'][i]
1869        ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
1870        vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
1871        dist = np.sqrt(np.sum(vec[i]**2))
1872        if not dist:
1873            return 0.
1874        vec[i] /= dist
1875    angle = acosd(np.sum(vec[0]*vec[1]))
1876#    GSASIIpath.IPyBreak()
1877    return angle
1878
1879def CalcAngleDeriv(angle_dict, angle_atoms, parmDict):
1880    if not len(parmDict):
1881        return None
1882    pId = angle_dict['pId']
1883    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1884    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1885    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1886    Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1887    Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1888    symNo = angle_dict['symNo']
1889    Tunit = angle_dict['cellNo']
1890    SGData = angle_dict['SGData']   
1891    deriv = getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData)
1892    return deriv
1893
1894def getSyXYZ(XYZ,ops,SGData):
1895    '''default doc
1896   
1897    :param type name: description
1898   
1899    :returns: type name: description
1900   
1901    '''
1902    XYZout = np.zeros_like(XYZ)
1903    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
1904        if op == '1':
1905            XYZout[i] = xyz
1906        else:
1907            oprs = op.split('+')
1908            unit = [0,0,0]
1909            if len(oprs)>1:
1910                unit = np.array(list(eval(oprs[1])))
1911            syop =int(oprs[0])
1912            inv = syop//abs(syop)
1913            syop *= inv
1914            cent = syop//100
1915            syop %= 100
1916            syop -= 1
1917            M,T = SGData['SGOps'][syop]
1918            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
1919    return XYZout
1920   
1921def getRestDist(XYZ,Amat):
1922    '''default doc string
1923   
1924    :param type name: description
1925   
1926    :returns: type name: description
1927   
1928    '''
1929    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1930   
1931def getRestDeriv(Func,XYZ,Amat,ops,SGData):
1932    '''default doc string
1933   
1934    :param type name: description
1935   
1936    :returns: type name: description
1937   
1938    '''
1939    deriv = np.zeros((len(XYZ),3))
1940    dx = 0.00001
1941    for j,xyz in enumerate(XYZ):
1942        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1943            XYZ[j] -= x
1944            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1945            XYZ[j] += 2*x
1946            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1947            XYZ[j] -= x
1948            deriv[j][i] = (d1-d2)/(2*dx)
1949    return deriv.flatten()
1950
1951def getRestAngle(XYZ,Amat):
1952    '''default doc string
1953   
1954    :param type name: description
1955   
1956    :returns: type name: description
1957   
1958    '''
1959   
1960    def calcVec(Ox,Tx,Amat):
1961        return np.inner(Amat,(Tx-Ox))
1962
1963    VecA = calcVec(XYZ[1],XYZ[0],Amat)
1964    VecA /= np.sqrt(np.sum(VecA**2))
1965    VecB = calcVec(XYZ[1],XYZ[2],Amat)
1966    VecB /= np.sqrt(np.sum(VecB**2))
1967    edge = VecB-VecA
1968    edge = np.sum(edge**2)
1969    angle = (2.-edge)/2.
1970    angle = max(angle,-1.)
1971    return acosd(angle)
1972   
1973def getRestPlane(XYZ,Amat):
1974    '''default doc string
1975   
1976    :param type name: description
1977   
1978    :returns: type name: description
1979   
1980    '''
1981    sumXYZ = np.zeros(3)
1982    for xyz in XYZ:
1983        sumXYZ += xyz
1984    sumXYZ /= len(XYZ)
1985    XYZ = np.array(XYZ)-sumXYZ
1986    XYZ = np.inner(Amat,XYZ).T
1987    Zmat = np.zeros((3,3))
1988    for i,xyz in enumerate(XYZ):
1989        Zmat += np.outer(xyz.T,xyz)
1990    Evec,Emat = nl.eig(Zmat)
1991    Evec = np.sqrt(Evec)/(len(XYZ)-3)
1992    Order = np.argsort(Evec)
1993    return Evec[Order[0]]
1994   
1995def getRestChiral(XYZ,Amat):   
1996    '''default doc string
1997   
1998    :param type name: description
1999   
2000    :returns: type name: description
2001   
2002    '''
2003    VecA = np.empty((3,3))   
2004    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
2005    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
2006    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
2007    return nl.det(VecA)
2008   
2009def getRestTorsion(XYZ,Amat):
2010    '''default doc string
2011   
2012    :param type name: description
2013   
2014    :returns: type name: description
2015   
2016    '''
2017    VecA = np.empty((3,3))
2018    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
2019    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
2020    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
2021    D = nl.det(VecA)
2022    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
2023    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
2024    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
2025    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
2026    Ang = 1.0
2027    if abs(P12) < 1.0 and abs(P23) < 1.0:
2028        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
2029    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
2030    return TOR
2031   
2032def calcTorsionEnergy(TOR,Coeff=[]):
2033    '''default doc string
2034   
2035    :param type name: description
2036   
2037    :returns: type name: description
2038   
2039    '''
2040    sum = 0.
2041    if len(Coeff):
2042        cof = np.reshape(Coeff,(3,3)).T
2043        delt = TOR-cof[1]
2044        delt = np.where(delt<-180.,delt+360.,delt)
2045        delt = np.where(delt>180.,delt-360.,delt)
2046        term = -cof[2]*delt**2
2047        val = cof[0]*np.exp(term/1000.0)
2048        pMax = cof[0][np.argmin(val)]
2049        Eval = np.sum(val)
2050        sum = Eval-pMax
2051    return sum,Eval
2052
2053def getTorsionDeriv(XYZ,Amat,Coeff):
2054    '''default doc string
2055   
2056    :param type name: description
2057   
2058    :returns: type name: description
2059   
2060    '''
2061    deriv = np.zeros((len(XYZ),3))
2062    dx = 0.00001
2063    for j,xyz in enumerate(XYZ):
2064        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2065            XYZ[j] -= x
2066            tor = getRestTorsion(XYZ,Amat)
2067            p1,d1 = calcTorsionEnergy(tor,Coeff)
2068            XYZ[j] += 2*x
2069            tor = getRestTorsion(XYZ,Amat)
2070            p2,d2 = calcTorsionEnergy(tor,Coeff)           
2071            XYZ[j] -= x
2072            deriv[j][i] = (p2-p1)/(2*dx)
2073    return deriv.flatten()
2074
2075def getRestRama(XYZ,Amat):
2076    '''Computes a pair of torsion angles in a 5 atom string
2077   
2078    :param nparray XYZ: crystallographic coordinates of 5 atoms
2079    :param nparray Amat: crystal to cartesian transformation matrix
2080   
2081    :returns: list (phi,psi) two torsion angles in degrees
2082   
2083    '''
2084    phi = getRestTorsion(XYZ[:5],Amat)
2085    psi = getRestTorsion(XYZ[1:],Amat)
2086    return phi,psi
2087   
2088def calcRamaEnergy(phi,psi,Coeff=[]):
2089    '''Computes pseudo potential energy from a pair of torsion angles and a
2090    numerical description of the potential energy surface. Used to create
2091    penalty function in LS refinement:     
2092    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
2093
2094    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
2095   
2096    :param float phi: first torsion angle (:math:`\\phi`)
2097    :param float psi: second torsion angle (:math:`\\psi`)
2098    :param list Coeff: pseudo potential coefficients
2099   
2100    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
2101      sum is used for penalty function.
2102   
2103    '''
2104    sum = 0.
2105    Eval = 0.
2106    if len(Coeff):
2107        cof = Coeff.T
2108        dPhi = phi-cof[1]
2109        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
2110        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
2111        dPsi = psi-cof[2]
2112        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
2113        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
2114        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
2115        val = cof[0]*np.exp(val/1000.)
2116        pMax = cof[0][np.argmin(val)]
2117        Eval = np.sum(val)
2118        sum = Eval-pMax
2119    return sum,Eval
2120
2121def getRamaDeriv(XYZ,Amat,Coeff):
2122    '''Computes numerical derivatives of torsion angle pair pseudo potential
2123    with respect of crystallographic atom coordinates of the 5 atom sequence
2124   
2125    :param nparray XYZ: crystallographic coordinates of 5 atoms
2126    :param nparray Amat: crystal to cartesian transformation matrix
2127    :param list Coeff: pseudo potential coefficients
2128   
2129    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
2130     crystallographic xyz coordinates.
2131   
2132    '''
2133    deriv = np.zeros((len(XYZ),3))
2134    dx = 0.00001
2135    for j,xyz in enumerate(XYZ):
2136        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2137            XYZ[j] -= x
2138            phi,psi = getRestRama(XYZ,Amat)
2139            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
2140            XYZ[j] += 2*x
2141            phi,psi = getRestRama(XYZ,Amat)
2142            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
2143            XYZ[j] -= x
2144            deriv[j][i] = (p2-p1)/(2*dx)
2145    return deriv.flatten()
2146
2147def getRestPolefig(ODFln,SamSym,Grid):
2148    '''default doc string
2149   
2150    :param type name: description
2151   
2152    :returns: type name: description
2153   
2154    '''
2155    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
2156    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
2157    R = np.where(R <= 1.,2.*atand(R),0.0)
2158    Z = np.zeros_like(R)
2159    Z = G2lat.polfcal(ODFln,SamSym,R,P)
2160    Z = np.reshape(Z,(Grid,Grid))
2161    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
2162
2163def getRestPolefigDerv(HKL,Grid,SHCoeff):
2164    '''default doc string
2165   
2166    :param type name: description
2167   
2168    :returns: type name: description
2169   
2170    '''
2171    pass
2172       
2173def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
2174    '''default doc string
2175   
2176    :param type name: description
2177   
2178    :returns: type name: description
2179   
2180    '''
2181    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
2182        TxT = inv*(np.inner(M,Tx)+T)+C+U
2183        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
2184       
2185    inv = Top/abs(Top)
2186    cent = abs(Top)//100
2187    op = abs(Top)%100-1
2188    M,T = SGData['SGOps'][op]
2189    C = SGData['SGCen'][cent]
2190    dx = .00001
2191    deriv = np.zeros(6)
2192    for i in [0,1,2]:
2193        Oxyz[i] -= dx
2194        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2195        Oxyz[i] += 2*dx
2196        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2197        Oxyz[i] -= dx
2198        Txyz[i] -= dx
2199        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2200        Txyz[i] += 2*dx
2201        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2202        Txyz[i] -= dx
2203    return deriv
2204   
2205def getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData):
2206   
2207    def calcAngle(Oxyz,ABxyz,Amat,Tunit,symNo,SGData):
2208        vec = np.zeros((2,3))
2209        for i in range(2):
2210            inv = 1
2211            if symNo[i] < 0:
2212                inv = -1
2213            cen = inv*symNo[i]//100
2214            op = inv*symNo[i]%100-1
2215            M,T = SGData['SGOps'][op]
2216            D = T*inv+SGData['SGCen'][cen]
2217            D += Tunit[i]
2218            ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
2219            vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
2220            dist = np.sqrt(np.sum(vec[i]**2))
2221            if not dist:
2222                return 0.
2223            vec[i] /= dist
2224        angle = acosd(np.sum(vec[0]*vec[1]))
2225    #    GSASIIpath.IPyBreak()
2226        return angle
2227       
2228    dx = .00001
2229    deriv = np.zeros(9)
2230    for i in [0,1,2]:
2231        Oxyz[i] -= dx
2232        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2233        Oxyz[i] += 2*dx
2234        deriv[i] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2235        Oxyz[i] -= dx
2236        Axyz[i] -= dx
2237        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2238        Axyz[i] += 2*dx
2239        deriv[i+3] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2240        Axyz[i] -= dx
2241        Bxyz[i] -= dx
2242        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2243        Bxyz[i] += 2*dx
2244        deriv[i+6] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2245        Bxyz[i] -= dx
2246    return deriv
2247   
2248def getAngSig(VA,VB,Amat,SGData,covData={}):
2249    '''default doc string
2250   
2251    :param type name: description
2252   
2253    :returns: type name: description
2254   
2255    '''
2256    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
2257        TxT = inv*(np.inner(M,Tx)+T)+C+U
2258        return np.inner(Amat,(TxT-Ox))
2259       
2260    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
2261        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
2262        VecA /= np.sqrt(np.sum(VecA**2))
2263        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
2264        VecB /= np.sqrt(np.sum(VecB**2))
2265        edge = VecB-VecA
2266        edge = np.sum(edge**2)
2267        angle = (2.-edge)/2.
2268        angle = max(angle,-1.)
2269        return acosd(angle)
2270       
2271    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
2272    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
2273    invA = invB = 1
2274    invA = TopA//abs(TopA)
2275    invB = TopB//abs(TopB)
2276    centA = abs(TopA)//100
2277    centB = abs(TopB)//100
2278    opA = abs(TopA)%100-1
2279    opB = abs(TopB)%100-1
2280    MA,TA = SGData['SGOps'][opA]
2281    MB,TB = SGData['SGOps'][opB]
2282    CA = SGData['SGCen'][centA]
2283    CB = SGData['SGCen'][centB]
2284    if 'covMatrix' in covData:
2285        covMatrix = covData['covMatrix']
2286        varyList = covData['varyList']
2287        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
2288        dx = .00001
2289        dadx = np.zeros(9)
2290        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2291        for i in [0,1,2]:
2292            OxA[i] -= dx
2293            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2294            OxA[i] += 2*dx
2295            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2296            OxA[i] -= dx
2297           
2298            TxA[i] -= dx
2299            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2300            TxA[i] += 2*dx
2301            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2302            TxA[i] -= dx
2303           
2304            TxB[i] -= dx
2305            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2306            TxB[i] += 2*dx
2307            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2308            TxB[i] -= dx
2309           
2310        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2311        if sigAng < 0.01:
2312            sigAng = 0.0
2313        return Ang,sigAng
2314    else:
2315        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
2316
2317def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
2318    '''default doc string
2319   
2320    :param type name: description
2321   
2322    :returns: type name: description
2323   
2324    '''
2325    def calcDist(Atoms,SyOps,Amat):
2326        XYZ = []
2327        for i,atom in enumerate(Atoms):
2328            Inv,M,T,C,U = SyOps[i]
2329            XYZ.append(np.array(atom[1:4]))
2330            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2331            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2332        V1 = XYZ[1]-XYZ[0]
2333        return np.sqrt(np.sum(V1**2))
2334       
2335    SyOps = []
2336    names = []
2337    for i,atom in enumerate(Oatoms):
2338        names += atom[-1]
2339        Op,unit = Atoms[i][-1]
2340        inv = Op//abs(Op)
2341        m,t = SGData['SGOps'][abs(Op)%100-1]
2342        c = SGData['SGCen'][abs(Op)//100]
2343        SyOps.append([inv,m,t,c,unit])
2344    Dist = calcDist(Oatoms,SyOps,Amat)
2345   
2346    sig = -0.001
2347    if 'covMatrix' in covData:
2348        dx = .00001
2349        dadx = np.zeros(6)
2350        for i in range(6):
2351            ia = i//3
2352            ix = i%3
2353            Oatoms[ia][ix+1] += dx
2354            a0 = calcDist(Oatoms,SyOps,Amat)
2355            Oatoms[ia][ix+1] -= 2*dx
2356            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2357        covMatrix = covData['covMatrix']
2358        varyList = covData['varyList']
2359        DistVcov = getVCov(names,varyList,covMatrix)
2360        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
2361        if sig < 0.001:
2362            sig = -0.001
2363   
2364    return Dist,sig
2365
2366def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
2367    '''default doc string
2368   
2369    :param type name: description
2370   
2371    :returns: type name: description
2372   
2373    '''
2374
2375    def calcAngle(Atoms,SyOps,Amat):
2376        XYZ = []
2377        for i,atom in enumerate(Atoms):
2378            Inv,M,T,C,U = SyOps[i]
2379            XYZ.append(np.array(atom[1:4]))
2380            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2381            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2382        V1 = XYZ[1]-XYZ[0]
2383        V1 /= np.sqrt(np.sum(V1**2))
2384        V2 = XYZ[1]-XYZ[2]
2385        V2 /= np.sqrt(np.sum(V2**2))
2386        V3 = V2-V1
2387        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2388        return acosd(cang)
2389
2390    SyOps = []
2391    names = []
2392    for i,atom in enumerate(Oatoms):
2393        names += atom[-1]
2394        Op,unit = Atoms[i][-1]
2395        inv = Op//abs(Op)
2396        m,t = SGData['SGOps'][abs(Op)%100-1]
2397        c = SGData['SGCen'][abs(Op)//100]
2398        SyOps.append([inv,m,t,c,unit])
2399    Angle = calcAngle(Oatoms,SyOps,Amat)
2400   
2401    sig = -0.01
2402    if 'covMatrix' in covData:
2403        dx = .00001
2404        dadx = np.zeros(9)
2405        for i in range(9):
2406            ia = i//3
2407            ix = i%3
2408            Oatoms[ia][ix+1] += dx
2409            a0 = calcAngle(Oatoms,SyOps,Amat)
2410            Oatoms[ia][ix+1] -= 2*dx
2411            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2412        covMatrix = covData['covMatrix']
2413        varyList = covData['varyList']
2414        AngVcov = getVCov(names,varyList,covMatrix)
2415        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2416        if sig < 0.01:
2417            sig = -0.01
2418   
2419    return Angle,sig
2420
2421def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
2422    '''default doc string
2423   
2424    :param type name: description
2425   
2426    :returns: type name: description
2427   
2428    '''
2429
2430    def calcTorsion(Atoms,SyOps,Amat):
2431       
2432        XYZ = []
2433        for i,atom in enumerate(Atoms):
2434            Inv,M,T,C,U = SyOps[i]
2435            XYZ.append(np.array(atom[1:4]))
2436            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2437            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2438        V1 = XYZ[1]-XYZ[0]
2439        V2 = XYZ[2]-XYZ[1]
2440        V3 = XYZ[3]-XYZ[2]
2441        V1 /= np.sqrt(np.sum(V1**2))
2442        V2 /= np.sqrt(np.sum(V2**2))
2443        V3 /= np.sqrt(np.sum(V3**2))
2444        M = np.array([V1,V2,V3])
2445        D = nl.det(M)
2446        P12 = np.dot(V1,V2)
2447        P13 = np.dot(V1,V3)
2448        P23 = np.dot(V2,V3)
2449        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2450        return Tors
2451           
2452    SyOps = []
2453    names = []
2454    for i,atom in enumerate(Oatoms):
2455        names += atom[-1]
2456        Op,unit = Atoms[i][-1]
2457        inv = Op//abs(Op)
2458        m,t = SGData['SGOps'][abs(Op)%100-1]
2459        c = SGData['SGCen'][abs(Op)//100]
2460        SyOps.append([inv,m,t,c,unit])
2461    Tors = calcTorsion(Oatoms,SyOps,Amat)
2462   
2463    sig = -0.01
2464    if 'covMatrix' in covData:
2465        dx = .00001
2466        dadx = np.zeros(12)
2467        for i in range(12):
2468            ia = i//3
2469            ix = i%3
2470            Oatoms[ia][ix+1] -= dx
2471            a0 = calcTorsion(Oatoms,SyOps,Amat)
2472            Oatoms[ia][ix+1] += 2*dx
2473            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2474            Oatoms[ia][ix+1] -= dx           
2475        covMatrix = covData['covMatrix']
2476        varyList = covData['varyList']
2477        TorVcov = getVCov(names,varyList,covMatrix)
2478        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
2479        if sig < 0.01:
2480            sig = -0.01
2481   
2482    return Tors,sig
2483       
2484def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
2485    '''default doc string
2486   
2487    :param type name: description
2488   
2489    :returns: type name: description
2490   
2491    '''
2492
2493    def calcDist(Atoms,SyOps,Amat):
2494        XYZ = []
2495        for i,atom in enumerate(Atoms):
2496            Inv,M,T,C,U = SyOps[i]
2497            XYZ.append(np.array(atom[1:4]))
2498            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2499            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2500        V1 = XYZ[1]-XYZ[0]
2501        return np.sqrt(np.sum(V1**2))
2502       
2503    def calcAngle(Atoms,SyOps,Amat):
2504        XYZ = []
2505        for i,atom in enumerate(Atoms):
2506            Inv,M,T,C,U = SyOps[i]
2507            XYZ.append(np.array(atom[1:4]))
2508            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2509            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2510        V1 = XYZ[1]-XYZ[0]
2511        V1 /= np.sqrt(np.sum(V1**2))
2512        V2 = XYZ[1]-XYZ[2]
2513        V2 /= np.sqrt(np.sum(V2**2))
2514        V3 = V2-V1
2515        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2516        return acosd(cang)
2517
2518    def calcTorsion(Atoms,SyOps,Amat):
2519       
2520        XYZ = []
2521        for i,atom in enumerate(Atoms):
2522            Inv,M,T,C,U = SyOps[i]
2523            XYZ.append(np.array(atom[1:4]))
2524            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2525            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2526        V1 = XYZ[1]-XYZ[0]
2527        V2 = XYZ[2]-XYZ[1]
2528        V3 = XYZ[3]-XYZ[2]
2529        V1 /= np.sqrt(np.sum(V1**2))
2530        V2 /= np.sqrt(np.sum(V2**2))
2531        V3 /= np.sqrt(np.sum(V3**2))
2532        M = np.array([V1,V2,V3])
2533        D = nl.det(M)
2534        P12 = np.dot(V1,V2)
2535        P13 = np.dot(V1,V3)
2536        P23 = np.dot(V2,V3)
2537        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2538        return Tors
2539           
2540    SyOps = []
2541    names = []
2542    for i,atom in enumerate(Oatoms):
2543        names += atom[-1]
2544        Op,unit = Atoms[i][-1]
2545        inv = Op//abs(Op)
2546        m,t = SGData['SGOps'][abs(Op)%100-1]
2547        c = SGData['SGCen'][abs(Op)//100]
2548        SyOps.append([inv,m,t,c,unit])
2549    M = len(Oatoms)
2550    if M == 2:
2551        Val = calcDist(Oatoms,SyOps,Amat)
2552    elif M == 3:
2553        Val = calcAngle(Oatoms,SyOps,Amat)
2554    else:
2555        Val = calcTorsion(Oatoms,SyOps,Amat)
2556   
2557    sigVals = [-0.001,-0.01,-0.01]
2558    sig = sigVals[M-3]
2559    if 'covMatrix' in covData:
2560        dx = .00001
2561        N = M*3
2562        dadx = np.zeros(N)
2563        for i in range(N):
2564            ia = i//3
2565            ix = i%3
2566            Oatoms[ia][ix+1] += dx
2567            if M == 2:
2568                a0 = calcDist(Oatoms,SyOps,Amat)
2569            elif M == 3:
2570                a0 = calcAngle(Oatoms,SyOps,Amat)
2571            else:
2572                a0 = calcTorsion(Oatoms,SyOps,Amat)
2573            Oatoms[ia][ix+1] -= 2*dx
2574            if M == 2:
2575                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2576            elif M == 3:
2577                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2578            else:
2579                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2580        covMatrix = covData['covMatrix']
2581        varyList = covData['varyList']
2582        Vcov = getVCov(names,varyList,covMatrix)
2583        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
2584        if sig < sigVals[M-3]:
2585            sig = sigVals[M-3]
2586   
2587    return Val,sig
2588       
2589def ValEsd(value,esd=0,nTZ=False):
2590    '''Format a floating point number with a given level of precision or
2591    with in crystallographic format with a "esd", as value(esd). If esd is
2592    negative the number is formatted with the level of significant figures
2593    appropriate if abs(esd) were the esd, but the esd is not included.
2594    if the esd is zero, approximately 6 significant figures are printed.
2595    nTZ=True causes "extra" zeros to be removed after the decimal place.
2596    for example:
2597
2598      * "1.235(3)" for value=1.2346 & esd=0.003
2599      * "1.235(3)e4" for value=12346. & esd=30
2600      * "1.235(3)e6" for value=0.12346e7 & esd=3000
2601      * "1.235" for value=1.2346 & esd=-0.003
2602      * "1.240" for value=1.2395 & esd=-0.003
2603      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
2604      * "1.23460" for value=1.2346 & esd=0.0
2605
2606    :param float value: number to be formatted
2607    :param float esd: uncertainty or if esd < 0, specifies level of
2608      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
2609    :param bool nTZ: True to remove trailing zeros (default is False)
2610    :returns: value(esd) or value as a string
2611
2612    '''
2613    # Note: this routine is Python 3 compatible -- I think
2614    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
2615    if math.isnan(value): # invalid value, bail out
2616        return '?'
2617    if math.isnan(esd): # invalid esd, treat as zero
2618        esd = 0
2619        esdoff = 5
2620#    if esd < 1.e-5:
2621#        esd = 0
2622#        esdoff = 5
2623    elif esd != 0:
2624        # transform the esd to a one or two digit integer
2625        l = math.log10(abs(esd)) % 1.
2626        if l < math.log10(cutoff): l+= 1.       
2627        intesd = int(round(10**l)) # esd as integer
2628        # determine the number of digits offset for the esd
2629        esdoff = int(round(math.log10(intesd*1./abs(esd))))
2630    else:
2631        esdoff = 5
2632    valoff = 0
2633    if abs(value) < abs(esdoff): # value is effectively zero
2634        pass
2635    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
2636        # where the digit offset is to the left of the decimal place or where too many
2637        # digits are needed
2638        if abs(value) > 1:
2639            valoff = int(math.log10(abs(value)))
2640        elif abs(value) > 0:
2641            valoff = int(math.log10(abs(value))-0.9999999)
2642        else:
2643            valoff = 0
2644    if esd != 0:
2645        if valoff+esdoff < 0:
2646            valoff = esdoff = 0
2647        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
2648    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
2649        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
2650    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
2651        if abs(value) > 0:           
2652            extra = -math.log10(abs(value))
2653        else:
2654            extra = 0
2655        if extra > 0: extra += 1
2656        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
2657    if esd > 0:
2658        out += ("({:d})").format(intesd)  # add the esd
2659    elif nTZ and '.' in out:
2660        out = out.rstrip('0')  # strip zeros to right of decimal
2661        out = out.rstrip('.')  # and decimal place when not needed
2662    if valoff != 0:
2663        out += ("e{:d}").format(valoff) # add an exponent, when needed
2664    return out
2665   
2666###############################################################################
2667##### Protein validation - "ERRATV2" analysis
2668###############################################################################
2669
2670def validProtein(Phase,old):
2671   
2672    def sumintact(intact):
2673        return {'CC':intact['CC'],'NN':intact['NN'],'OO':intact['OO'],
2674        'CN':(intact['CN']+intact['NC']),'CO':(intact['CO']+intact['OC']),
2675        'NO':(intact['NO']+intact['ON'])}
2676       
2677    resNames = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
2678        'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
2679# data from errat.f
2680    b1_old = np.array([ 
2681        [1154.343,  600.213, 1051.018, 1132.885,  960.738],
2682        [600.213, 1286.818, 1282.042,  957.156,  612.789],
2683        [1051.018, 1282.042, 3519.471,  991.974, 1226.491],
2684        [1132.885,  957.156,  991.974, 1798.672,  820.355],
2685        [960.738,  612.789, 1226.491,  820.355, 2428.966]
2686        ])
2687    avg_old = np.array([ 0.225, 0.281, 0.071, 0.237, 0.044])    #Table 1 3.5A Obsd. Fr. p 1513
2688# data taken from erratv2.ccp
2689    b1 = np.array([
2690          [5040.279078850848200,        3408.805141583649400,   4152.904423767300600,   4236.200004171890200,   5054.781210204625500], 
2691          [3408.805141583648900,        8491.906094010220800,   5958.881777877950300,   1521.387352718486200,   4304.078200827221700], 
2692          [4152.904423767301500,        5958.881777877952100,   7637.167089335050100,   6620.715738223072500,   5287.691183798410700], 
2693          [4236.200004171890200,        1521.387352718486200,   6620.715738223072500,   18368.343774298410000,  4050.797811118806700], 
2694          [5054.781210204625500,        4304.078200827220800,   5287.691183798409800,   4050.797811118806700,   6666.856740479164700]])
2695    avg = np.array([0.192765509919262, 0.195575208778518, 0.275322406824210, 0.059102357035642, 0.233154192767480])
2696    General = Phase['General']
2697    Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
2698    cx,ct,cs,cia = General['AtomPtrs']
2699    Atoms = Phase['Atoms']
2700    cartAtoms = []
2701    xyzmin = 999.*np.ones(3)
2702    xyzmax = -999.*np.ones(3)
2703    #select residue atoms,S,Se --> O make cartesian
2704    for atom in Atoms:
2705        if atom[1] in resNames:
2706            cartAtoms.append(atom[:cx+3])
2707            if atom[4].strip() in ['S','Se']:
2708                if not old:
2709                    continue        #S,Se skipped for erratv2?
2710                cartAtoms[-1][3] = 'Os'
2711                cartAtoms[-1][4] = 'O'
2712            cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
2713            cartAtoms[-1].append(atom[cia+8])
2714    XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
2715    xyzmin = np.array([np.min(XYZ.T[i]) for i in [0,1,2]])
2716    xyzmax = np.array([np.max(XYZ.T[i]) for i in [0,1,2]])
2717    nbox = list(np.array(np.ceil((xyzmax-xyzmin)/4.),dtype=int))+[15,]
2718    Boxes = np.zeros(nbox,dtype=int)
2719    iBox = np.array([np.trunc((XYZ.T[i]-xyzmin[i])/4.) for i in [0,1,2]],dtype=int).T
2720    for ib,box in enumerate(iBox):  #put in a try for too many atoms in box (IndexError)?
2721        try:
2722            Boxes[box[0],box[1],box[2],0] += 1
2723            Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
2724        except IndexError:
2725            print('too many atoms in box' )
2726            continue
2727    #Box content checks with errat.f $ erratv2.cpp ibox1 arrays
2728    indices = (-1,0,1)
2729    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) 
2730    dsmax = 3.75**2
2731    if old:
2732        dsmax = 3.5**2
2733    chains = []
2734    resIntAct = []
2735    chainIntAct = []
2736    res = []
2737    resNames = []
2738    resIDs = {}
2739    resname = []
2740    resID = {}
2741    newChain = True
2742    intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2743    for ia,atom in enumerate(cartAtoms):
2744        jntact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2745        if atom[2] not in chains:   #get chain id & save residue sequence from last chain
2746            chains.append(atom[2])
2747            if len(resIntAct):
2748                resIntAct.append(sumintact(intact))
2749                chainIntAct.append(resIntAct)
2750                resNames += resname
2751                resIDs.update(resID)
2752                res = []
2753                resname = []
2754                resID = {}
2755                resIntAct = []
2756                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2757                newChain = True
2758        if atom[0] not in res:  #new residue, get residue no.
2759            if res and int(res[-1]) != int(atom[0])-1:  #a gap in chain - not new chain
2760                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2761                ires = int(res[-1])
2762                for i in range(int(atom[0])-ires-1):
2763                    res.append(str(ires+i+1))
2764                    resname.append('')
2765                    resIntAct.append(sumintact(intact))
2766            res.append(atom[0])
2767            name = '%s-%s%s'%(atom[2],atom[0],atom[1])
2768            resname.append(name)
2769            resID[name] = atom[-1]
2770            if not newChain:
2771                resIntAct.append(sumintact(intact))
2772            intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2773            newChain = False
2774        ibox = iBox[ia]         #box location of atom
2775        tgts = []
2776        for unit in Units:      #assemble list of all possible target atoms
2777            jbox = ibox+unit
2778            if np.all(jbox>=0) and np.all((jbox-nbox[:3])<0):               
2779                tgts += list(Boxes[jbox[0],jbox[1],jbox[2]])
2780        tgts = list(set(tgts))
2781        tgts = [tgt for tgt in tgts if atom[:3] != cartAtoms[tgt][:3]]    #exclude same residue
2782        tgts = [tgt for tgt in tgts if np.sum((XYZ[ia]-XYZ[tgt])**2) < dsmax]
2783        ires = int(atom[0])
2784        if old:
2785            if atom[3].strip() == 'C':
2786                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2787            elif atom[3].strip() == 'N':
2788                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])]
2789            elif atom[3].strip() == 'CA':
2790                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2791        else:
2792            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]]
2793            if atom[3].strip() == 'C':
2794                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) == ires+1)]
2795            elif atom[3].strip() == 'N':
2796                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'C' and int(cartAtoms[tgt][0]) == ires-1)]
2797        for tgt in tgts:
2798            dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
2799            mult = 1.0
2800            if dsqt > 3.25 and not old:
2801                mult = 2.*(3.75-dsqt)
2802            intype = atom[4].strip()+cartAtoms[tgt][4].strip()
2803            if 'S' not in intype:
2804                intact[intype] += mult
2805                jntact[intype] += mult
2806#        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']
2807    resNames += resname
2808    resIDs.update(resID)
2809    resIntAct.append(sumintact(intact))
2810    chainIntAct.append(resIntAct)
2811    chainProb = []
2812    for ich,chn in enumerate(chains):
2813        IntAct = chainIntAct[ich]
2814        nRes = len(IntAct)
2815        Probs = [0.,0.,0.,0.]   #skip 1st 4 residues in chain
2816        for i in range(4,nRes-4):
2817            if resNames[i]:
2818                mtrx = np.zeros(5)
2819                summ = 0.
2820                for j in range(i-4,i+5):
2821                    summ += np.sum(np.array(list(IntAct[j].values())))
2822                    if old:
2823                        mtrx[0] += IntAct[j]['CC']
2824                        mtrx[1] += IntAct[j]['CO']
2825                        mtrx[2] += IntAct[j]['NN']
2826                        mtrx[3] += IntAct[j]['NO']
2827                        mtrx[4] += IntAct[j]['OO']
2828                    else:
2829                        mtrx[0] += IntAct[j]['CC']
2830                        mtrx[1] += IntAct[j]['CN']
2831                        mtrx[2] += IntAct[j]['CO']
2832                        mtrx[3] += IntAct[j]['NN']
2833                        mtrx[4] += IntAct[j]['NO']
2834                mtrx /= summ
2835    #            print i+1,mtrx*summ
2836                if old:
2837                    mtrx -= avg_old
2838                    prob = np.inner(np.inner(mtrx,b1_old),mtrx)
2839                else:
2840                    mtrx -= avg
2841                    prob = np.inner(np.inner(mtrx,b1),mtrx)
2842            else:       #skip the gaps
2843                prob = 0.0
2844            Probs.append(prob)
2845        Probs += 4*[0.,]        #skip last 4 residues in chain
2846        chainProb += Probs
2847    return resNames,chainProb,resIDs
2848   
2849################################################################################
2850##### Texture fitting stuff
2851################################################################################
2852
2853def FitTexture(General,Gangls,refData,keyList,pgbar):
2854    import pytexture as ptx
2855    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2856   
2857    def printSpHarm(textureData,SHtextureSig):
2858        print ('\n Spherical harmonics texture: Order:' + str(textureData['Order']))
2859        names = ['omega','chi','phi']
2860        namstr = '  names :'
2861        ptstr =  '  values:'
2862        sigstr = '  esds  :'
2863        for name in names:
2864            namstr += '%12s'%('Sample '+name)
2865            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2866            if 'Sample '+name in SHtextureSig:
2867                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2868            else:
2869                sigstr += 12*' '
2870        print (namstr)
2871        print (ptstr)
2872        print (sigstr)
2873        print ('\n Texture coefficients:')
2874        SHcoeff = textureData['SH Coeff'][1]
2875        SHkeys = list(SHcoeff.keys())
2876        nCoeff = len(SHcoeff)
2877        nBlock = nCoeff//10+1
2878        iBeg = 0
2879        iFin = min(iBeg+10,nCoeff)
2880        for block in range(nBlock):
2881            namstr = '  names :'
2882            ptstr =  '  values:'
2883            sigstr = '  esds  :'
2884            for name in SHkeys[iBeg:iFin]:
2885                if 'C' in name:
2886                    namstr += '%12s'%(name)
2887                    ptstr += '%12.3f'%(SHcoeff[name])
2888                    if name in SHtextureSig:
2889                        sigstr += '%12.3f'%(SHtextureSig[name])
2890                    else:
2891                        sigstr += 12*' '
2892            print (namstr)
2893            print (ptstr)
2894            print (sigstr)
2895            iBeg += 10
2896            iFin = min(iBeg+10,nCoeff)
2897           
2898    def Dict2Values(parmdict, varylist):
2899        '''Use before call to leastsq to setup list of values for the parameters
2900        in parmdict, as selected by key in varylist'''
2901        return [parmdict[key] for key in varylist] 
2902       
2903    def Values2Dict(parmdict, varylist, values):
2904        ''' Use after call to leastsq to update the parameter dictionary with
2905        values corresponding to keys in varylist'''
2906        parmdict.update(list(zip(varylist,values)))
2907       
2908    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2909        parmDict.update(list(zip(varyList,values)))
2910        Mat = np.empty(0)
2911        sumObs = 0
2912        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2913        for hist in Gangls.keys():
2914            Refs = refData[hist]
2915            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2916            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2917#            wt = 1./np.max(Refs[:,4],.25)
2918            sumObs += np.sum(wt*Refs[:,5])
2919            Refs[:,6] = 1.
2920            H = Refs[:,:3]
2921            phi,beta = G2lat.CrsAng(H,cell,SGData)
2922            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2923            for item in parmDict:
2924                if 'C' in item:
2925                    L,M,N = eval(item.strip('C'))
2926                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2927                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2928                    Lnorm = G2lat.Lnorm(L)
2929                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2930            mat = wt*(Refs[:,5]-Refs[:,6])
2931            Mat = np.concatenate((Mat,mat))
2932        sumD = np.sum(np.abs(Mat))
2933        R = min(100.,100.*sumD/sumObs)
2934        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2935        print (' Residual: %.3f%%'%(R))
2936        return Mat
2937       
2938    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2939        Mat = np.empty(0)
2940        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2941        for hist in Gangls.keys():
2942            mat = np.zeros((len(varyList),len(refData[hist])))
2943            Refs = refData[hist]
2944            H = Refs[:,:3]
2945            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2946#            wt = 1./np.max(Refs[:,4],.25)
2947            phi,beta = G2lat.CrsAng(H,cell,SGData)
2948            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2949            for j,item in enumerate(varyList):
2950                if 'C' in item:
2951                    L,M,N = eval(item.strip('C'))
2952                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2953                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2954                    Lnorm = G2lat.Lnorm(L)
2955                    mat[j] = -wt*Lnorm*Kcl*Ksl
2956                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2957                        try:
2958                            l = varyList.index(itema)
2959                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2960                        except ValueError:
2961                            pass
2962            if len(Mat):
2963                Mat = np.concatenate((Mat,mat.T))
2964            else:
2965                Mat = mat.T
2966        print ('deriv')
2967        return Mat
2968
2969    print (' Fit texture for '+General['Name'])
2970    SGData = General['SGData']
2971    cell = General['Cell'][1:7]
2972    Texture = General['SH Texture']
2973    if not Texture['Order']:
2974        return 'No spherical harmonics coefficients'
2975    varyList = []
2976    parmDict = copy.copy(Texture['SH Coeff'][1])
2977    for item in ['Sample omega','Sample chi','Sample phi']:
2978        parmDict[item] = Texture[item][1]
2979        if Texture[item][0]:
2980            varyList.append(item)
2981    if Texture['SH Coeff'][0]:
2982        varyList += list(Texture['SH Coeff'][1].keys())
2983    while True:
2984        begin = time.time()
2985        values =  np.array(Dict2Values(parmDict, varyList))
2986        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
2987            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
2988        ncyc = int(result[2]['nfev']//2)
2989        if ncyc:
2990            runtime = time.time()-begin   
2991            chisq = np.sum(result[2]['fvec']**2)
2992            Values2Dict(parmDict, varyList, result[0])
2993            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
2994            print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(result[2]['fvec']),len(varyList)))
2995            print ('refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2996            try:
2997                sig = np.sqrt(np.diag(result[1])*GOF)
2998                if np.any(np.isnan(sig)):
2999                    print ('*** Least squares aborted - some invalid esds possible ***')
3000                break                   #refinement succeeded - finish up!
3001            except ValueError:          #result[1] is None on singular matrix
3002                print ('**** Refinement failed - singular matrix ****')
3003                return None
3004        else:
3005            break
3006   
3007    if ncyc:
3008        for parm in parmDict:
3009            if 'C' in parm:
3010                Texture['SH Coeff'][1][parm] = parmDict[parm]
3011            else:
3012                Texture[parm][1] = parmDict[parm] 
3013        sigDict = dict(zip(varyList,sig))
3014        printSpHarm(Texture,sigDict)
3015       
3016    return None
3017   
3018################################################################################
3019##### Fourier & charge flip stuff
3020################################################################################
3021
3022def adjHKLmax(SGData,Hmax):
3023    '''default doc string
3024   
3025    :param type name: description
3026   
3027    :returns: type name: description
3028   
3029    '''
3030    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3031        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
3032        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
3033        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
3034    else:
3035        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
3036        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
3037        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
3038
3039def OmitMap(data,reflDict,pgbar=None):
3040    '''default doc string
3041   
3042    :param type name: description
3043   
3044    :returns: type name: description
3045   
3046    '''
3047    generalData = data['General']
3048    if not generalData['Map']['MapType']:
3049        print ('**** ERROR - Fourier map not defined')
3050        return
3051    mapData = generalData['Map']
3052    dmin = mapData['Resolution']
3053    SGData = generalData['SGData']
3054    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3055    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3056    cell = generalData['Cell'][1:8]       
3057    A = G2lat.cell2A(cell[:6])
3058    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3059    adjHKLmax(SGData,Hmax)
3060    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3061    time0 = time.time()
3062    for iref,ref in enumerate(reflDict['RefList']):
3063        if ref[4] >= dmin:
3064            Fosq,Fcsq,ph = ref[8:11]
3065            Uniq = np.inner(ref[:3],SGMT)
3066            Phi = np.inner(ref[:3],SGT)
3067            for i,hkl in enumerate(Uniq):        #uses uniq
3068                hkl = np.asarray(hkl,dtype='i')
3069                dp = 360.*Phi[i]                #and phi
3070                a = cosd(ph+dp)
3071                b = sind(ph+dp)
3072                phasep = complex(a,b)
3073                phasem = complex(a,-b)
3074                if '2Fo-Fc' in mapData['MapType']:
3075                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
3076                else:
3077                    F = np.sqrt(Fosq)
3078                h,k,l = hkl+Hmax
3079                Fhkl[h,k,l] = F*phasep
3080                h,k,l = -hkl+Hmax
3081                Fhkl[h,k,l] = F*phasem
3082    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3083    M = np.mgrid[0:4,0:4,0:4]
3084    blkIds = np.array(list(zip(M[0].flatten(),M[1].flatten(),M[2].flatten())))
3085    iBeg = blkIds*rho0.shape//4
3086    iFin = (blkIds+1)*rho0.shape//4
3087    rho_omit = np.zeros_like(rho0)
3088    nBlk = 0
3089    for iB,iF in zip(iBeg,iFin):
3090        rho1 = np.copy(rho0)
3091        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
3092        Fnew = fft.ifftshift(fft.ifftn(rho1))
3093        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
3094        phase = Fnew/np.absolute(Fnew)
3095        OFhkl = np.absolute(Fhkl)*phase
3096        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
3097        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]])
3098        nBlk += 1
3099        pgbar.Update(nBlk)
3100    mapData['rho'] = np.real(rho_omit)/cell[6]
3101    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3102    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3103    print ('Omit map time: %.4f no. elements: %d'%(time.time()-time0,Fhkl.size))
3104    return mapData
3105   
3106def FourierMap(data,reflDict):
3107    '''default doc string
3108   
3109    :param type name: description
3110   
3111    :returns: type name: description
3112   
3113    '''
3114    generalData = data['General']
3115    mapData = generalData['Map']
3116    dmin = mapData['Resolution']
3117    SGData = generalData['SGData']
3118    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3119    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3120    cell = generalData['Cell'][1:8]       
3121    A = G2lat.cell2A(cell[:6])
3122    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3123    adjHKLmax(SGData,Hmax)
3124    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3125#    Fhkl[0,0,0] = generalData['F000X']
3126    time0 = time.time()
3127    for iref,ref in enumerate(reflDict['RefList']):
3128        if ref[4] > dmin:
3129            Fosq,Fcsq,ph = ref[8:11]
3130            Uniq = np.inner(ref[:3],SGMT)
3131            Phi = np.inner(ref[:3],SGT)
3132            for i,hkl in enumerate(Uniq):        #uses uniq
3133                hkl = np.asarray(hkl,dtype='i')
3134                dp = 360.*Phi[i]                #and phi
3135                a = cosd(ph+dp)
3136                b = sind(ph+dp)
3137                phasep = complex(a,b)
3138                phasem = complex(a,-b)
3139                if 'Fobs' in mapData['MapType']:
3140                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
3141                    h,k,l = hkl+Hmax
3142                    Fhkl[h,k,l] = F*phasep
3143                    h,k,l = -hkl+Hmax
3144                    Fhkl[h,k,l] = F*phasem
3145                elif 'Fcalc' in mapData['MapType']:
3146                    F = np.sqrt(Fcsq)
3147                    h,k,l = hkl+Hmax
3148                    Fhkl[h,k,l] = F*phasep
3149                    h,k,l = -hkl+Hmax
3150                    Fhkl[h,k,l] = F*phasem
3151                elif 'delt-F' in mapData['MapType']:
3152                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3153                    h,k,l = hkl+Hmax
3154                    Fhkl[h,k,l] = dF*phasep
3155                    h,k,l = -hkl+Hmax
3156                    Fhkl[h,k,l] = dF*phasem
3157                elif '2*Fo-Fc' in mapData['MapType']:
3158                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3159                    h,k,l = hkl+Hmax
3160                    Fhkl[h,k,l] = F*phasep
3161                    h,k,l = -hkl+Hmax
3162                    Fhkl[h,k,l] = F*phasem
3163                elif 'Patterson' in mapData['MapType']:
3164                    h,k,l = hkl+Hmax
3165                    Fhkl[h,k,l] = complex(Fosq,0.)
3166                    h,k,l = -hkl+Hmax
3167                    Fhkl[h,k,l] = complex(Fosq,0.)
3168    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3169    print ('Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size))
3170    mapData['Type'] = reflDict['Type']
3171    mapData['rho'] = np.real(rho)
3172    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3173    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3174   
3175def Fourier4DMap(data,reflDict):
3176    '''default doc string
3177   
3178    :param type name: description
3179   
3180    :returns: type name: description
3181   
3182    '''
3183    generalData = data['General']
3184    map4DData = generalData['4DmapData']
3185    mapData = generalData['Map']
3186    dmin = mapData['Resolution']
3187    SGData = generalData['SGData']
3188    SSGData = generalData['SSGData']
3189    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3190    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3191    cell = generalData['Cell'][1:8]       
3192    A = G2lat.cell2A(cell[:6])
3193    maxM = 4
3194    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
3195    adjHKLmax(SGData,Hmax)
3196    Hmax = np.asarray(Hmax,dtype='i')+1
3197    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3198    time0 = time.time()
3199    for iref,ref in enumerate(reflDict['RefList']):
3200        if ref[5] > dmin:
3201            Fosq,Fcsq,ph = ref[9:12]
3202            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
3203            Uniq = np.inner(ref[:4],SSGMT)
3204            Phi = np.inner(ref[:4],SSGT)
3205            for i,hkl in enumerate(Uniq):        #uses uniq
3206                hkl = np.asarray(hkl,dtype='i')
3207                dp = 360.*Phi[i]                #and phi
3208                a = cosd(ph+dp)
3209                b = sind(ph+dp)
3210                phasep = complex(a,b)
3211                phasem = complex(a,-b)
3212                if 'Fobs' in mapData['MapType']:
3213                    F = np.sqrt(Fosq)
3214                    h,k,l,m = hkl+Hmax
3215                    Fhkl[h,k,l,m] = F*phasep
3216                    h,k,l,m = -hkl+Hmax
3217                    Fhkl[h,k,l,m] = F*phasem
3218                elif 'Fcalc' in mapData['MapType']:
3219                    F = np.sqrt(Fcsq)
3220                    h,k,l,m = hkl+Hmax
3221                    Fhkl[h,k,l,m] = F*phasep
3222                    h,k,l,m = -hkl+Hmax
3223                    Fhkl[h,k,l,m] = F*phasem                   
3224                elif 'delt-F' in mapData['MapType']:
3225                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
3226                    h,k,l,m = hkl+Hmax
3227                    Fhkl[h,k,l,m] = dF*phasep
3228                    h,k,l,m = -hkl+Hmax
3229                    Fhkl[h,k,l,m] = dF*phasem
3230    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
3231    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
3232    map4DData['rho'] = np.real(SSrho)
3233    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3234    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3235    map4DData['Type'] = reflDict['Type']
3236    mapData['Type'] = reflDict['Type']
3237    mapData['rho'] = np.real(rho)
3238    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3239    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3240    print ('Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size))
3241
3242# map printing for testing purposes
3243def printRho(SGLaue,rho,rhoMax):                         
3244    '''default doc string
3245   
3246    :param type name: description
3247   
3248    :returns: type name: description
3249   
3250    '''
3251    dim = len(rho.shape)
3252    if dim == 2:
3253        ix,jy = rho.shape
3254        for j in range(jy):
3255            line = ''
3256            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3257                line += (jy-j)*'  '
3258            for i in range(ix):
3259                r = int(100*rho[i,j]/rhoMax)
3260                line += '%4d'%(r)
3261            print (line+'\n')
3262    else:
3263        ix,jy,kz = rho.shape
3264        for k in range(kz):
3265            print ('k = %d'%k)
3266            for j in range(jy):
3267                line = ''
3268                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3269                    line += (jy-j)*'  '
3270                for i in range(ix):
3271                    r = int(100*rho[i,j,k]/rhoMax)
3272                    line += '%4d'%(r)
3273                print (line+'\n')
3274## keep this
3275               
3276def findOffset(SGData,A,Fhkl):   
3277    '''default doc string
3278   
3279    :param type name: description
3280   
3281    :returns: type name: description
3282   
3283    '''
3284    if SGData['SpGrp'] == 'P 1':
3285        return [0,0,0]   
3286    hklShape = Fhkl.shape
3287    hklHalf = np.array(hklShape)/2
3288    sortHKL = np.argsort(Fhkl.flatten())
3289    Fdict = {}
3290    for hkl in sortHKL:
3291        HKL = np.unravel_index(hkl,hklShape)
3292        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
3293        if F == 0.:
3294            break
3295        Fdict['%.6f'%(np.absolute(F))] = hkl
3296    Flist = np.flipud(np.sort(list(Fdict.keys())))
3297    F = str(1.e6)
3298    i = 0
3299    DH = []
3300    Dphi = []
3301    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3302    for F in Flist:
3303        hkl = np.unravel_index(Fdict[F],hklShape)
3304        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
3305            continue
3306        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
3307        Uniq = np.array(Uniq,dtype='i')
3308        Phi = np.array(Phi)
3309        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
3310        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3311        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
3312        ang0 = np.angle(Fh0,deg=True)/360.
3313        for H,phi in list(zip(Uniq,Phi))[1:]:
3314            ang = (np.angle(Fhkl[int(H[0]),int(H[1]),int(H[2])],deg=True)/360.-phi)
3315            dH = H-hkl
3316            dang = ang-ang0
3317            DH.append(dH)
3318            Dphi.append((dang+.5) % 1.0)
3319        if i > 20 or len(DH) > 30:
3320            break
3321        i += 1
3322    DH = np.array(DH)
3323    print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))
3324    Dphi = np.array(Dphi)
3325    steps = np.array(hklShape)
3326    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
3327    XYZ = np.array(list(zip(X.flatten(),Y.flatten(),Z.flatten())))
3328    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
3329    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3330    hist,bins = np.histogram(Mmap,bins=1000)
3331#    for i,item in enumerate(hist[:10]):
3332#        print item,bins[i]
3333    chisq = np.min(Mmap)
3334    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3335    print (' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2]))
3336#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3337    return DX
3338   
3339def ChargeFlip(data,reflDict,pgbar):
3340    '''default doc string
3341   
3342    :param type name: description
3343   
3344    :returns: type name: description
3345   
3346    '''
3347    generalData = data['General']
3348    mapData = generalData['Map']
3349    flipData = generalData['Flip']
3350    FFtable = {}
3351    if 'None' not in flipData['Norm element']:
3352        normElem = flipData['Norm element'].upper()
3353        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3354        for ff in FFs:
3355            if ff['Symbol'] == normElem:
3356                FFtable.update(ff)
3357    dmin = flipData['Resolution']
3358    SGData = generalData['SGData']
3359    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3360    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3361    cell = generalData['Cell'][1:8]       
3362    A = G2lat.cell2A(cell[:6])
3363    Vol = cell[6]
3364    im = 0
3365    if generalData['Modulated'] == True:
3366        im = 1
3367    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3368    adjHKLmax(SGData,Hmax)
3369    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3370    time0 = time.time()
3371    for iref,ref in enumerate(reflDict['RefList']):
3372        dsp = ref[4+im]
3373        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
3374            continue
3375        if dsp > dmin:
3376            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3377            if FFtable:
3378                SQ = 0.25/dsp**2
3379                ff *= G2el.ScatFac(FFtable,SQ)[0]
3380            if ref[8+im] > 0.:         #use only +ve Fobs**2
3381                E = np.sqrt(ref[8+im])/ff
3382            else:
3383                E = 0.
3384            ph = ref[10]
3385            ph = rn.uniform(0.,360.)
3386            Uniq = np.inner(ref[:3],SGMT)
3387            Phi = np.inner(ref[:3],SGT)
3388            for i,hkl in enumerate(Uniq):        #uses uniq
3389                hkl = np.asarray(hkl,dtype='i')
3390                dp = 360.*Phi[i]                #and phi
3391                a = cosd(ph+dp)
3392                b = sind(ph+dp)
3393                phasep = complex(a,b)
3394                phasem = complex(a,-b)
3395                h,k,l = hkl+Hmax
3396                Ehkl[h,k,l] = E*phasep
3397                h,k,l = -hkl+Hmax
3398                Ehkl[h,k,l] = E*phasem
3399#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3400    testHKL = np.array(flipData['testHKL'])+Hmax
3401    CEhkl = copy.copy(Ehkl)
3402    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3403    Emask = ma.getmask(MEhkl)
3404    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3405    Ncyc = 0
3406    old = np.seterr(all='raise')
3407    twophases = []
3408    while True:       
3409        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3410        CEsig = np.std(CErho)
3411        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3412        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3413        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3414        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3415        phase = CFhkl/np.absolute(CFhkl)
3416        twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
3417        CEhkl = np.absolute(Ehkl)*phase
3418        Ncyc += 1
3419        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3420        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3421        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3422        if Rcf < 5.:
3423            break
3424        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3425        if not GoOn or Ncyc > 10000:
3426            break
3427    np.seterr(**old)
3428    print (' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size))
3429    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
3430    print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))
3431    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3432       
3433    mapData['Rcf'] = Rcf
3434    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3435    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3436    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3437    mapData['Type'] = reflDict['Type']
3438    return mapData,twophases
3439   
3440def findSSOffset(SGData,SSGData,A,Fhklm):   
3441    '''default doc string
3442   
3443    :param type name: description
3444   
3445    :returns: type name: description
3446   
3447    '''
3448    if SGData['SpGrp'] == 'P 1':
3449        return [0,0,0,0]   
3450    hklmShape = Fhklm.shape
3451    hklmHalf = np.array(hklmShape)/2
3452    sortHKLM = np.argsort(Fhklm.flatten())
3453    Fdict = {}
3454    for hklm in sortHKLM:
3455        HKLM = np.unravel_index(hklm,hklmShape)
3456        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
3457        if F == 0.:
3458            break
3459        Fdict['%.6f'%(np.absolute(F))] = hklm
3460    Flist = np.flipud(np.sort(list(Fdict.keys())))
3461    F = str(1.e6)
3462    i = 0
3463    DH = []
3464    Dphi = []
3465    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3466    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3467    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3468    for F in Flist:
3469        hklm = np.unravel_index(Fdict[F],hklmShape)
3470        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
3471            continue
3472        Uniq = np.inner(hklm-hklmHalf,SSGMT)
3473        Phi = np.inner(hklm-hklmHalf,SSGT)
3474        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
3475        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3476        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
3477        ang0 = np.angle(Fh0,deg=True)/360.
3478        for H,phi in list(zip(Uniq,Phi))[1:]:
3479            H = np.array(H,dtype=int)
3480            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
3481            dH = H-hklm
3482            dang = ang-ang0
3483            DH.append(dH)
3484            Dphi.append((dang+.5) % 1.0)
3485        if i > 20 or len(DH) > 30:
3486            break
3487        i += 1
3488    DH = np.array(DH)
3489    print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))
3490    Dphi = np.array(Dphi)
3491    steps = np.array(hklmShape)
3492    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]]
3493    XYZT = np.array(list(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten())))
3494    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
3495    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3496    hist,bins = np.histogram(Mmap,bins=1000)
3497#    for i,item in enumerate(hist[:10]):
3498#        print item,bins[i]
3499    chisq = np.min(Mmap)
3500    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3501    print (' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3]))
3502#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3503    return DX
3504   
3505def SSChargeFlip(data,reflDict,pgbar):
3506    '''default doc string
3507   
3508    :param type name: description
3509   
3510    :returns: type name: description
3511   
3512    '''
3513    generalData = data['General']
3514    mapData = generalData['Map']
3515    map4DData = {}
3516    flipData = generalData['Flip']
3517    FFtable = {}
3518    if 'None' not in flipData['Norm element']:
3519        normElem = flipData['Norm element'].upper()
3520        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3521        for ff in FFs:
3522            if ff['Symbol'] == normElem:
3523                FFtable.update(ff)
3524    dmin = flipData['Resolution']
3525    SGData = generalData['SGData']
3526    SSGData = generalData['SSGData']
3527    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3528    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3529    cell = generalData['Cell'][1:8]       
3530    A = G2lat.cell2A(cell[:6])
3531    Vol = cell[6]
3532    maxM = 4
3533    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3534    adjHKLmax(SGData,Hmax)
3535    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3536    time0 = time.time()
3537    for iref,ref in enumerate(reflDict['RefList']):
3538        dsp = ref[5]
3539        if dsp > dmin:
3540            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3541            if FFtable:
3542                SQ = 0.25/dsp**2
3543                ff *= G2el.ScatFac(FFtable,SQ)[0]
3544            if ref[9] > 0.:         #use only +ve Fobs**2
3545                E = np.sqrt(ref[9])/ff
3546            else:
3547                E = 0.
3548            ph = ref[11]
3549            ph = rn.uniform(0.,360.)
3550            Uniq = np.inner(ref[:4],SSGMT)
3551            Phi = np.inner(ref[:4],SSGT)
3552            for i,hklm in enumerate(Uniq):        #uses uniq
3553                hklm = np.asarray(hklm,dtype='i')
3554                dp = 360.*Phi[i]                #and phi
3555                a = cosd(ph+dp)
3556                b = sind(ph+dp)
3557                phasep = complex(a,b)
3558                phasem = complex(a,-b)
3559                h,k,l,m = hklm+Hmax
3560                Ehkl[h,k,l,m] = E*phasep
3561                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
3562                Ehkl[h,k,l,m] = E*phasem
3563#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3564    CEhkl = copy.copy(Ehkl)
3565    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3566    Emask = ma.getmask(MEhkl)
3567    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3568    Ncyc = 0
3569    old = np.seterr(all='raise')
3570    while True:       
3571        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3572        CEsig = np.std(CErho)
3573        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3574        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3575        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3576        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3577        phase = CFhkl/np.absolute(CFhkl)
3578        CEhkl = np.absolute(Ehkl)*phase
3579        Ncyc += 1
3580        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3581        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3582        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3583        if Rcf < 5.:
3584            break
3585        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3586        if not GoOn or Ncyc > 10000:
3587            break
3588    np.seterr(**old)
3589    print (' Charge flip time: %.4f no. elements: %d'%(time.time()-time0,Ehkl.size))
3590    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
3591    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
3592    print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))
3593    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3594
3595    mapData['Rcf'] = Rcf
3596    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3597    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3598    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3599    mapData['Type'] = reflDict['Type']
3600
3601    map4DData['Rcf'] = Rcf
3602    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))
3603    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3604    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3605    map4DData['Type'] = reflDict['Type']
3606    return mapData,map4DData
3607   
3608def getRho(xyz,mapData):
3609    ''' get scattering density at a point by 8-point interpolation
3610    param xyz: coordinate to be probed
3611    param: mapData: dict of map data
3612   
3613    :returns: density at xyz
3614    '''
3615    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3616    if not len(mapData):
3617        return 0.0
3618    rho = copy.copy(mapData['rho'])     #don't mess up original
3619    if not len(rho):
3620        return 0.0
3621    mapShape = np.array(rho.shape)
3622    mapStep = 1./mapShape
3623    X = np.array(xyz)%1.    #get into unit cell
3624    I = np.array(X*mapShape,dtype='int')
3625    D = X-I*mapStep         #position inside map cell
3626    D12 = D[0]*D[1]
3627    D13 = D[0]*D[2]
3628    D23 = D[1]*D[2]
3629    D123 = np.prod(D)
3630    Rho = rollMap(rho,-I)    #shifts map so point is in corner
3631    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]+  \
3632        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
3633        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
3634        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
3635    return R
3636       
3637def SearchMap(generalData,drawingData,Neg=False):
3638    '''Does a search of a density map for peaks meeting the criterion of peak
3639    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3640    mapData is data['General']['mapData']; the map is also in mapData.
3641
3642    :param generalData: the phase data structure; includes the map
3643    :param drawingData: the drawing data structure
3644    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
3645
3646    :returns: (peaks,mags,dzeros) where
3647
3648        * peaks : ndarray
3649          x,y,z positions of the peaks found in the map
3650        * mags : ndarray
3651          the magnitudes of the peaks
3652        * dzeros : ndarray
3653          the distance of the peaks from  the unit cell origin
3654        * dcent : ndarray
3655          the distance of the peaks from  the unit cell center
3656
3657    '''       
3658    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3659   
3660    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
3661   
3662    def fixSpecialPos(xyz,SGData,Amat):
3663        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3664        X = []
3665        xyzs = [equiv[0] for equiv in equivs]
3666        for x in xyzs:
3667            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3668                X.append(x)
3669        if len(X) > 1:
3670            return np.average(X,axis=0)
3671        else:
3672            return xyz
3673       
3674    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3675        Mag,x0,y0,z0,sig = parms
3676        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
3677#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
3678        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
3679       
3680    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3681        Mag,x0,y0,z0,sig = parms
3682        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3683        return M
3684       
3685    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3686        Mag,x0,y0,z0,sig = parms
3687        dMdv = np.zeros(([5,]+list(rX.shape)))
3688        delt = .01
3689        for i in range(5):
3690            parms[i] -= delt
3691            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3692            parms[i] += 2.*delt
3693            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3694            parms[i] -= delt
3695            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3696        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3697        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3698        dMdv = np.reshape(dMdv,(5,rX.size))
3699        Hess = np.inner(dMdv,dMdv)
3700       
3701        return Vec,Hess
3702       
3703    SGData = generalData['SGData']
3704    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3705    peaks = []
3706    mags = []
3707    dzeros = []
3708    dcent = []
3709    try:
3710        mapData = generalData['Map']
3711        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3712        if Neg:
3713            rho = -copy.copy(mapData['rho'])     #flip +/-
3714        else:
3715            rho = copy.copy(mapData['rho'])     #don't mess up original
3716        mapHalf = np.array(rho.shape)/2
3717        res = mapData['Resolution']
3718        incre = np.array(rho.shape,dtype=np.float)
3719        step = max(1.0,1./res)+1
3720        steps = np.array((3*[step,]),dtype='int32')
3721    except KeyError:
3722        print ('**** ERROR - Fourier map not defined')
3723        return peaks,mags
3724    rhoMask = ma.array(rho,mask=(rho<contLevel))
3725    indices = (-1,0,1)
3726    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3727    for roll in rolls:
3728        if np.any(roll):
3729            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3730    indx = np.transpose(rhoMask.nonzero())
3731    peaks = indx/incre
3732    mags = rhoMask[rhoMask.nonzero()]
3733    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3734        rho = rollMap(rho,ind)
3735        rMM = mapHalf-steps
3736        rMP = mapHalf+steps+1
3737        rhoPeak = rho[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3738        peakInt = np.sum(rhoPeak)*res**3
3739        rX,rY,rZ = np.mgrid[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3740        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
3741        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3742            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3743        x1 = result[0]
3744        if not np.any(x1 < 0):
3745            peak = (np.array(x1[1:4])-ind)/incre
3746        peak = fixSpecialPos(peak,SGData,Amat)
3747        rho = rollMap(rho,-ind)
3748    cent = np.ones(3)*.5     
3749    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3750    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3751    if Neg:     #want negative magnitudes for negative peaks
3752        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3753    else:
3754        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3755   
3756def sortArray(data,pos,reverse=False):
3757    '''data is a list of items
3758    sort by pos in list; reverse if True
3759    '''
3760    T = []
3761    for i,M in enumerate(data):
3762        try:
3763            T.append((M[pos],i))
3764        except IndexError:
3765            return data
3766    D = dict(zip(T,data))
3767    T.sort()
3768    if reverse:
3769        T.reverse()
3770    X = []
3771    for key in T:
3772        X.append(D[key])
3773    return X
3774
3775def PeaksEquiv(data,Ind):
3776    '''Find the equivalent map peaks for those selected. Works on the
3777    contents of data['Map Peaks'].
3778
3779    :param data: the phase data structure
3780    :param list Ind: list of selected peak indices
3781    :returns: augmented list of peaks including those related by symmetry to the
3782      ones in Ind
3783
3784    '''       
3785    def Duplicate(xyz,peaks,Amat):
3786        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3787            return True
3788        return False
3789                           
3790    generalData = data['General']
3791    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3792    SGData = generalData['SGData']
3793    mapPeaks = data['Map Peaks']
3794    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3795    Indx = {}
3796    for ind in Ind:
3797        xyz = np.array(mapPeaks[ind][1:4])
3798        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3799        for jnd,xyz in enumerate(XYZ):       
3800            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3801    Ind = []
3802    for ind in Indx:
3803        if Indx[ind]:
3804            Ind.append(ind)
3805    return Ind
3806               
3807def PeaksUnique(data,Ind):
3808    '''Finds the symmetry unique set of peaks from those selected. Works on the
3809    contents of data['Map Peaks'].
3810
3811    :param data: the phase data structure
3812    :param list Ind: list of selected peak indices
3813    :returns: the list of symmetry unique peaks from among those given in Ind
3814
3815    '''       
3816#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3817
3818    def noDuplicate(xyz,peaks,Amat):
3819        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3820            return False
3821        return True
3822                           
3823    generalData = data['General']
3824    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3825    SGData = generalData['SGData']
3826    mapPeaks = data['Map Peaks']
3827    Indx = {}
3828    XYZ = {}
3829    for ind in Ind:
3830        XYZ[ind] = np.array(mapPeaks[ind][1:4])
3831        Indx[ind] = True
3832    for ind in Ind:
3833        if Indx[ind]:
3834            xyz = XYZ[ind]
3835            for jnd in Ind:
3836                if ind != jnd and Indx[jnd]:                       
3837                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3838                    xyzs = np.array([equiv[0] for equiv in Equiv])
3839                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3840    Ind = []
3841    for ind in Indx:
3842        if Indx[ind]:
3843            Ind.append(ind)
3844    return Ind
3845   
3846################################################################################
3847##### single peak fitting profile fxn stuff
3848################################################################################
3849
3850def getCWsig(ins,pos):
3851    '''get CW peak profile sigma^2
3852   
3853    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3854      as values only
3855    :param float pos: 2-theta of peak
3856    :returns: float getCWsig: peak sigma^2
3857   
3858    '''
3859    tp = tand(pos/2.0)
3860    return ins['U']*tp**2+ins['V']*tp+ins['W']
3861   
3862def getCWsigDeriv(pos):
3863    '''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3864   
3865    :param float pos: 2-theta of peak
3866   
3867    :returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3868   
3869    '''
3870    tp = tand(pos/2.0)
3871    return tp**2,tp,1.0
3872   
3873def getCWgam(ins,pos):
3874    '''get CW peak profile gamma
3875   
3876    :param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3877      as values only
3878    :param float pos: 2-theta of peak
3879    :returns: float getCWgam: peak gamma
3880   
3881    '''
3882    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)+ins['Z']
3883   
3884def getCWgamDeriv(pos):
3885    '''get derivatives of CW peak profile gamma wrt X, Y & Z
3886   
3887    :param float pos: 2-theta of peak
3888   
3889    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3890   
3891    '''
3892    return 1./cosd(pos/2.0),tand(pos/2.0),1.0
3893   
3894def getTOFsig(ins,dsp):
3895    '''get TOF peak profile sigma^2
3896   
3897    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3898      as values only
3899    :param float dsp: d-spacing of peak
3900   
3901    :returns: float getTOFsig: peak sigma^2
3902   
3903    '''
3904    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']*dsp
3905   
3906def getTOFsigDeriv(dsp):
3907    '''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3908   
3909    :param float dsp: d-spacing of peak
3910   
3911    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3912   
3913    '''
3914    return 1.0,dsp**2,dsp**4,dsp
3915   
3916def getTOFgamma(ins,dsp):
3917    '''get TOF peak profile gamma
3918   
3919    :param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3920      as values only
3921    :param float dsp: d-spacing of peak
3922   
3923    :returns: float getTOFgamma: peak gamma
3924   
3925    '''
3926    return ins['Z']+ins['X']*dsp+ins['Y']*dsp**2
3927   
3928def getTOFgammaDeriv(dsp):
3929    '''get derivatives of TOF peak profile gamma wrt X, Y & Z
3930   
3931    :param float dsp: d-spacing of peak
3932   
3933    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3934   
3935    '''
3936    return dsp,dsp**2,1.0
3937   
3938def getTOFbeta(ins,dsp):
3939    '''get TOF peak profile beta
3940   
3941    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3942      as values only
3943    :param float dsp: d-spacing of peak
3944   
3945    :returns: float getTOFbeta: peak beat
3946   
3947    '''
3948    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
3949   
3950def getTOFbetaDeriv(dsp):
3951    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3952   
3953    :param float dsp: d-spacing of peak
3954   
3955    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3956   
3957    '''
3958    return 1.0,1./dsp**4,1./dsp**2
3959   
3960def getTOFalpha(ins,dsp):
3961    '''get TOF peak profile alpha
3962   
3963    :param dict ins: instrument parameters with at least 'alpha'
3964      as values only
3965    :param float dsp: d-spacing of peak
3966   
3967    :returns: flaot getTOFalpha: peak alpha
3968   
3969    '''
3970    return ins['alpha']/dsp
3971   
3972def getTOFalphaDeriv(dsp):
3973    '''get derivatives of TOF peak profile beta wrt alpha
3974   
3975    :param float dsp: d-spacing of peak
3976   
3977    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
3978   
3979    '''
3980    return 1./dsp
3981   
3982def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
3983    '''set starting peak parameters for single peak fits from plot selection or auto selection
3984   
3985    :param dict Parms: instrument parameters dictionary
3986    :param dict Parms2: table lookup for TOF profile coefficients
3987    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
3988    :param float mag: peak top magnitude from pick
3989    :param bool ifQ: True if pos in Q
3990    :param bool useFit: True if use fitted CW Parms values (not defaults)
3991   
3992    :returns: list XY: peak list entry:
3993        for CW: [pos,0,mag,1,sig,0,gam,0]
3994        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3995        NB: mag refinement set by default, all others off
3996   
3997    '''
3998    ind = 0
3999    if useFit:
4000        ind = 1
4001    ins = {}
4002    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an else
4003        for x in ['U','V','W','X','Y','Z']:
4004            ins[x] = Parms[x][ind]
4005        if ifQ:                              #qplot - convert back to 2-theta
4006            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
4007        sig = getCWsig(ins,pos)
4008        gam = getCWgam(ins,pos)           
4009        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
4010    else:
4011        if ifQ:
4012            dsp = 2.*np.pi/pos
4013            pos = Parms['difC']*dsp
4014        else:
4015            dsp = pos/Parms['difC'][1]
4016        if 'Pdabc' in Parms2:
4017            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4018                ins[x] = Parms[x][ind]
4019            Pdabc = Parms2['Pdabc'].T
4020            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
4021            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
4022        else:
4023            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4024                ins[x] = Parms[x][ind]
4025            alp = getTOFalpha(ins,dsp)
4026            bet = getTOFbeta(ins,dsp)
4027        sig = getTOFsig(ins,dsp)
4028        gam = getTOFgamma(ins,dsp)
4029        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4030    return XY
4031   
4032################################################################################
4033##### MC/SA stuff
4034################################################################################
4035
4036#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
4037# Original Author: Travis Oliphant 2002
4038# Bug-fixes in 2006 by Tim Leslie
4039
4040
4041import numpy
4042from numpy import asarray, exp, squeeze, sign, \
4043     all, shape, array, where
4044from numpy import random
4045
4046#__all__ = ['anneal']
4047
4048_double_min = numpy.finfo(float).min
4049_double_max = numpy.finfo(float).max
4050class base_schedule(object):
4051    def __init__(self):
4052        self.dwell = 20
4053        self.lower = 0.
4054        self.upper = 1.
4055        self.Ninit = 50
4056        self.accepted = 0
4057        self.tests = 0
4058        self.feval = 0
4059        self.k = 0
4060        self.T = None
4061
4062    def init(self, **options):
4063        self.__dict__.update(options)
4064        self.lower = asarray(self.lower)
4065        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
4066        self.upper = asarray(self.upper)
4067        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
4068        self.k = 0
4069        self.accepted = 0
4070        self.feval = 0
4071        self.tests = 0
4072
4073    def getstart_temp(self, best_state):
4074        """ Find a matching starting temperature and starting parameters vector
4075        i.e. find x0 such that func(x0) = T0.
4076
4077        :param best_state: _state
4078            A _state object to store the function value and x0 found.
4079
4080        :returns: x0 : array
4081            The starting parameters vector.
4082        """
4083
4084        assert(not self.dims is None)
4085        lrange = self.lower
4086        urange = self.upper
4087        fmax = _double_min
4088        fmin = _double_max
4089        for _ in range(self.Ninit):
4090            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
4091            fval = self.func(x0, *self.args)
4092            self.feval += 1
4093            if fval > fmax:
4094                fmax = fval
4095            if fval < fmin:
4096                fmin = fval
4097                best_state.cost = fval
4098                best_state.x = array(x0)
4099
4100        self.T0 = (fmax-fmin)*1.5
4101        return best_state.x
4102       
4103    def set_range(self,x0,frac):
4104        delrange = frac*(self.upper-self.lower)
4105        self.upper = x0+delrange
4106        self.lower = x0-delrange
4107
4108    def accept_test(self, dE):
4109        T = self.T
4110        self.tests += 1
4111        if dE < 0:
4112            self.accepted += 1
4113            return 1
4114        p = exp(-dE*1.0/T)
4115        if (p > random.uniform(0.0, 1.0)):
4116            self.accepted += 1
4117            return 1
4118        return 0
4119
4120    def update_guess(self, x0):
4121        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
4122
4123    def update_temp(self, x0):
4124        pass
4125
4126class fast_sa(base_schedule):
4127    def init(self, **options):
4128        self.__dict__.update(options)
4129
4130    def update_guess(self, x0):
4131        x0 = asarray(x0)
4132        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4133        T = self.T
4134        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4135        xnew = xc*(self.upper - self.lower)+self.lower
4136        return xnew
4137
4138    def update_temp(self):
4139        self.T = self.T0*exp(-self.c * self.k**(self.quench))
4140        self.k += 1
4141        return
4142
4143class log_sa(base_schedule):        #OK
4144
4145    def init(self,**options):
4146        self.__dict__.update(options)
4147       
4148    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
4149        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4150        T = self.T
4151        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4152        xnew = xc*(self.upper - self.lower)+self.lower
4153        return xnew
4154       
4155    def update_temp(self):
4156        self.k += 1
4157        self.T = self.T0*self.slope**self.k
4158       
4159class _state(object):
4160    def __init__(self):
4161        self.x = None
4162        self.cost = None
4163
4164def makeTsched(data):
4165    if data['Algorithm'] == 'fast':
4166        sched = fast_sa()
4167        sched.quench = data['fast parms'][0]
4168        sched.c = data['fast parms'][1]
4169    elif data['Algorithm'] == 'log':
4170        sched = log_sa()
4171        sched.slope = data['log slope']
4172    sched.T0 = data['Annealing'][0]
4173    if not sched.T0:
4174        sched.T0 = 50.
4175    Tf = data['Annealing'][1]
4176    if not Tf:
4177        Tf = 0.001
4178    Tsched = [sched.T0,]
4179    while Tsched[-1] > Tf:
4180        sched.update_temp()
4181        Tsched.append(sched.T)
4182    return Tsched[1:]
4183   
4184def anneal(func, x0, args=(), schedule='fast', 
4185           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4186           feps=1e-6, quench=1.0, c=1.0,
4187           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4188           ranRange=0.10,autoRan=False,dlg=None):
4189    """Minimize a function using simulated annealing.
4190
4191    Schedule is a schedule class implementing the annealing schedule.
4192    Available ones are 'fast', 'cauchy', 'boltzmann'
4193
4194    :param callable func: f(x, \*args)
4195        Function to be optimized.
4196    :param ndarray x0:
4197        Initial guess.
4198    :param tuple args:
4199        Extra parameters to `func`.
4200    :param base_schedule schedule:
4201        Annealing schedule to use (a class).
4202    :param float T0:
4203        Initial Temperature (estimated as 1.2 times the largest
4204        cost-function deviation over random points in the range).
4205    :param float Tf:
4206        Final goal temperature.
4207    :param int maxeval:
4208        Maximum function evaluations.
4209    :param int maxaccept:
4210        Maximum changes to accept.
4211    :param int maxiter:
4212        Maximum cooling iterations.
4213    :param float feps:
4214        Stopping relative error tolerance for the function value in
4215        last four coolings.
4216    :param float quench,c:
4217        Parameters to alter fast_sa schedule.
4218    :param float/ndarray lower,upper:
4219        Lower and upper bounds on `x`.
4220    :param int dwell:
4221        The number of times to search the space at each temperature.
4222    :param float slope:
4223        Parameter for log schedule
4224    :param bool ranStart=False:
4225        True for set 10% of ranges about x
4226
4227    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4228
4229     * xmin (ndarray): Point giving smallest value found.
4230     * Jmin (float): Minimum value of function found.
4231     * T (float): Final temperature.
4232     * feval (int): Number of function evaluations.
4233     * iters (int): Number of cooling iterations.
4234     * accept (int): Number of tests accepted.
4235     * retval (int): Flag indicating stopping condition:
4236
4237              *  0: Points no longer changing
4238              *  1: Cooled to final temperature
4239              *  2: Maximum function evaluations
4240              *  3: Maximum cooling iterations reached
4241              *  4: Maximum accepted query locations reached
4242              *  5: Final point not the minimum amongst encountered points
4243
4244    *Notes*:
4245    Simulated annealing is a random algorithm which uses no derivative
4246    information from the function being optimized. In practice it has
4247    been more useful in discrete optimization than continuous
4248    optimization, as there are usually better algorithms for continuous
4249    optimization problems.
4250
4251    Some experimentation by trying the difference temperature
4252    schedules and altering their parameters is likely required to
4253    obtain good performance.
4254
4255    The randomness in the algorithm comes from random sampling in numpy.
4256    To obtain the same results you can call numpy.random.seed with the
4257    same seed immediately before calling scipy.optimize.anneal.
4258
4259    We give a brief description of how the three temperature schedules
4260    generate new points and vary their temperature. Temperatures are
4261    only updated with iterations in the outer loop. The inner loop is
4262    over range(dwell), and new points are generated for
4263    every iteration in the inner loop. (Though whether the proposed
4264    new points are accepted is probabilistic.)
4265
4266    For readability, let d denote the dimension of the inputs to func.
4267    Also, let x_old denote the previous state, and k denote the
4268    iteration number of the outer loop. All other variables not
4269    defined below are input variables to scipy.optimize.anneal itself.
4270
4271    In the 'fast' schedule the updates are ::
4272
4273        u ~ Uniform(0, 1, size=d)
4274        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4275        xc = y * (upper - lower)
4276        x_new = x_old + xc
4277
4278        T_new = T0 * exp(-c * k**quench)
4279
4280    """
4281   
4282    ''' Scipy license:
4283        Copyright (c) 2001, 2002 Enthought, Inc.
4284    All rights reserved.
4285   
4286    Copyright (c) 2003-2016 SciPy Developers.
4287    All rights reserved.
4288   
4289    Redistribution and use in source and binary forms, with or without
4290    modification, are permitted provided that the following conditions are met:
4291   
4292      a. Redistributions of source code must retain the above copyright notice,
4293         this list of conditions and the following disclaimer.
4294      b. Redistributions in binary form must reproduce the above copyright
4295         notice, this list of conditions and the following disclaimer in the
4296         documentation and/or other materials provided with the distribution.
4297      c. Neither the name of Enthought nor the names of the SciPy Developers
4298         may be used to endorse or promote products derived from this software
4299         without specific prior written permission.
4300    '''
4301    x0 = asarray(x0)
4302    lower = asarray(lower)
4303    upper = asarray(upper)
4304
4305    schedule = eval(schedule+'_sa()')
4306    #   initialize the schedule
4307    schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4308        c=c, quench=quench, dwell=dwell, slope=slope)
4309
4310    current_state, last_state, best_state = _state(), _state(), _state()
4311    if ranStart:
4312        schedule.set_range(x0,ranRange)
4313    if T0 is None:
4314        x0 = schedule.getstart_temp(best_state)
4315    else:
4316        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4317        best_state.x = None
4318        best_state.cost = numpy.Inf
4319
4320    last_state.x = asarray(x0).copy()
4321    fval = func(x0,*args)
4322    schedule.feval += 1
4323    last_state.cost = fval
4324    if last_state.cost < best_state.cost:
4325        best_state.cost = fval
4326        best_state.x = asarray(x0).copy()
4327    schedule.T = schedule.T0
4328    fqueue = [100, 300, 500, 700]
4329    iters = 1
4330    keepGoing = True
4331    bestn = 0
4332    while keepGoing:
4333        retval = 0
4334        for n in range(dwell):
4335            current_state.x = schedule.update_guess(last_state.x)
4336            current_state.cost = func(current_state.x,*args)
4337            schedule.feval += 1
4338
4339            dE = current_state.cost - last_state.cost
4340            if schedule.accept_test(dE):
4341                last_state.x = current_state.x.copy()
4342                last_state.cost = current_state.cost
4343                if last_state.cost < best_state.cost:
4344                    best_state.x = last_state.x.copy()
4345                    best_state.cost = last_state.cost
4346                    bestn = n
4347                    if best_state.cost < 1.0 and autoRan:
4348                        schedule.set_range(x0,best_state.cost/2.)                       
4349        if dlg:
4350            GoOn = dlg.Update(min(100.,best_state.cost*100),
4351                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4352                    'Best trial:',bestn,  \
4353                    'MC/SA Residual:',best_state.cost*100,'%', \
4354                    ))[0]
4355            if not GoOn:
4356                best_state.x = last_state.x.copy()
4357                best_state.cost = last_state.cost
4358                retval = 5
4359        schedule.update_temp()
4360        iters += 1
4361        # Stopping conditions
4362        # 0) last saved values of f from each cooling step
4363        #     are all very similar (effectively cooled)
4364        # 1) Tf is set and we are below it
4365        # 2) maxeval is set and we are past it
4366        # 3) maxiter is set and we are past it
4367        # 4) maxaccept is set and we are past it
4368        # 5) user canceled run via progress bar
4369
4370        fqueue.append(squeeze(last_state.cost))
4371        fqueue.pop(0)
4372        af = asarray(fqueue)*1.0
4373        if retval == 5:
4374            print (' User terminated run; incomplete MC/SA')
4375            keepGoing = False
4376            break
4377        if all(abs((af-af[0])/af[0]) < feps):
4378            retval = 0
4379            if abs(af[-1]-best_state.cost) > feps*10:
4380                retval = 5
4381                print (" Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1))
4382            break
4383        if (Tf is not None) and (schedule.T < Tf):
4384#            print ' Minimum T reached in %d steps'%(iters-1)
4385            retval = 1
4386            break
4387        if (maxeval is not None) and (schedule.feval > maxeval):
4388            retval = 2
4389            break
4390        if (iters > maxiter):
4391            print  (" Warning: Maximum number of iterations exceeded.")
4392            retval = 3
4393            break
4394        if (maxaccept is not None) and (schedule.accepted > maxaccept):
4395            retval = 4
4396            break
4397
4398    return best_state.x, best_state.cost, schedule.T, \
4399           schedule.feval, iters, schedule.accepted, retval
4400
4401def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4402    outlist = []
4403    timelist = []
4404    nsflist = []
4405    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
4406    for n in range(iCyc):
4407        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False)         #mcsa result,time,rcov
4408        outlist.append(result[0])
4409        timelist.append(result[1])
4410        nsflist.append(result[2])
4411        print (' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1]))
4412    out_q.put(outlist)
4413    out_t.put(timelist)
4414    out_n.put(nsflist)
4415
4416def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4417    import multiprocessing as mp
4418   
4419    out_q = mp.Queue()
4420    out_t = mp.Queue()
4421    out_n = mp.Queue()
4422    procs = []
4423    totsftime = 0.
4424    totnsf = 0
4425    iCyc = np.zeros(nprocs)
4426    for i in range(nCyc):
4427        iCyc[i%nprocs] += 1
4428    for i in range(nprocs):
4429        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4430        procs.append(p)
4431        p.start()
4432    resultlist = []
4433    for i in range(nprocs):
4434        resultlist += out_q.get()
4435        totsftime += np.sum(out_t.get())
4436        totnsf += np.sum(out_n.get())
4437    for p in procs:
4438        p.join()
4439    return resultlist,totsftime,totnsf
4440
4441def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4442    '''default doc string
4443   
4444    :param type name: description
4445   
4446    :returns: type name: description
4447    '''
4448   
4449    class RandomDisplacementBounds(object):
4450        """random displacement with bounds"""
4451        def __init__(self, xmin, xmax, stepsize=0.5):
4452            self.xmin = xmin
4453            self.xmax = xmax
4454            self.stepsize = stepsize
4455   
4456        def __call__(self, x):
4457            """take a random step but ensure the new position is within the bounds"""
4458            while True:
4459                # this could be done in a much more clever way, but it will work for example purposes
4460                steps = self.xmax-self.xmin
4461                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
4462                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4463                    break
4464            return xnew
4465   
4466    global tsum,nsum
4467    tsum = 0.
4468    nsum = 0
4469   
4470    def getMDparms(item,pfx,parmDict,varyList):
4471        parmDict[pfx+'MDaxis'] = item['axis']
4472        parmDict[pfx+'MDval'] = item['Coef'][0]
4473        if item['Coef'][1]:
4474            varyList += [pfx+'MDval',]
4475            limits = item['Coef'][2]
4476            lower.append(limits[0])
4477            upper.append(limits[1])
4478                       
4479    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4480        parmDict[pfx+'Atype'] = item['atType']
4481        aTypes |= set([item['atType'],]) 
4482        pstr = ['Ax','Ay','Az']
4483        XYZ = [0,0,0]
4484        for i in range(3):
4485            name = pfx+pstr[i]
4486            parmDict[name] = item['Pos'][0][i]
4487            XYZ[i] = parmDict[name]
4488            if item['Pos'][1][i]:
4489                varyList += [name,]
4490                limits = item['Pos'][2][i]
4491                lower.append(limits[0])
4492                upper.append(limits[1])
4493        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4494           
4495    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4496        parmDict[mfx+'MolCent'] = item['MolCent']
4497        parmDict[mfx+'RBId'] = item['RBId']
4498        pstr = ['Px','Py','Pz']
4499        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4500        for i in range(3):
4501            name = mfx+pstr[i]
4502            parmDict[name] = item['Pos'