source: trunk/GSASIImath.py @ 2328

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

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