source: trunk/GSASIImath.py @ 2846

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

implement SVD inversion (Moore Penrose invert routine) in Hessian LSQ
clean up commented out code replaced by ValidatedTextCtrl?
set default dM/M to 0.001 (from 0.0001) - more reasonable convergence criterion

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