source: trunk/GSASIImath.py @ 2110

Last change on this file since 2110 was 2110, checked in by vondreele, 6 years ago

add menu item for global setting of wave vary parameters - TBD
split StructureFactorDeriv? over twins & incommensurate
do block refl for the nontwin/powder one

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