source: trunk/GSASIImath.py @ 3754

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

begin work on incommensurate mag. structure factors
fix problem of magnetic reflection extinctions
replace analytical derivatives of magnetic moments with numerical ones using deltM = 0.0001

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