source: trunk/GSASIImath.py @ 2062

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

exclude non variable SS parms from constraint lists
define ZigZag? & Block position modulations; eliminate Sawtooth (it's a variant of ZigZag?)
fix SS names in constraint lists
implement ZigZag? & Block position wave plots
implement ZigZag? & Block atom motion in structure plots
add movie making option (hidden - no file output for it yet)
fix LS I/O for ZigZag? & Block waves

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