source: trunk/GSASIImath.py @ 3896

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

enormous speed improvement for contour plot - not dependent on map size

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