source: trunk/GSASIImath.py @ 2544

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

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