source: trunk/GSASIImath.py @ 2326

Last change on this file since 2326 was 2326, checked in by vondreele, 7 years ago

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