source: trunk/GSASIImath.py @ 1531

Last change on this file since 1531 was 1531, checked in by vondreele, 8 years ago

fix getVCOV; nonexistant parms return zero.
fix space group output in PrintDistAngle?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 107.1 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2014-10-21 17:28:29 +0000 (Tue, 21 Oct 2014) $
5# $Author: vondreele $
6# $Revision: 1531 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1531 2014-10-21 17:28:29Z 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: 1531 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import GSASIIpwd as G2pwd
34import numpy.fft as fft
35import pypowder as pyd
36
37sind = lambda x: np.sin(x*np.pi/180.)
38cosd = lambda x: np.cos(x*np.pi/180.)
39tand = lambda x: np.tan(x*np.pi/180.)
40asind = lambda x: 180.*np.arcsin(x)/np.pi
41acosd = lambda x: 180.*np.arccos(x)/np.pi
42atand = lambda x: 180.*np.arctan(x)/np.pi
43atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
44
45def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
46   
47    """
48    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
49    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
50   
51    ::
52   
53                    Nobs
54        x = arg min(sum(func(y)**2,axis=0))
55                    y=0
56
57    :param function func: callable method or function
58        should take at least one (possibly length N vector) argument and
59        returns M floating point numbers.
60    :param np.ndarray x0: The starting estimate for the minimization of length N
61    :param function Hess: callable method or function
62        A required function or method to compute the weighted vector and Hessian for func.
63        It must be a symmetric NxN array
64    :param tuple args: Any extra arguments to func are placed in this tuple.
65    :param float ftol: Relative error desired in the sum of squares.
66    :param float xtol: Relative error desired in the approximate solution.
67    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
68        until other limits are met (ftol, xtol)
69    :param bool Print: True for printing results (residuals & times) by cycle
70
71    :returns: (x,cov_x,infodict) where
72
73      * x : ndarray
74        The solution (or the result of the last iteration for an unsuccessful
75        call).
76      * cov_x : ndarray
77        Uses the fjac and ipvt optional outputs to construct an
78        estimate of the jacobian around the solution.  ``None`` if a
79        singular matrix encountered (indicates very flat curvature in
80        some direction).  This matrix must be multiplied by the
81        residual standard deviation to get the covariance of the
82        parameter estimates -- see curve_fit.
83      * infodict : dict
84        a dictionary of optional outputs with the keys:
85
86         * 'fvec' : the function evaluated at the output
87         * 'num cyc':
88         * 'nfev':
89         * 'lamMax':
90         * 'psing':
91           
92    """
93               
94    ifConverged = False
95    deltaChi2 = -10.
96    x0 = np.array(x0, ndmin=1)      #might be redundant?
97    n = len(x0)
98    if type(args) != type(()):
99        args = (args,)
100       
101    icycle = 0
102    One = np.ones((n,n))
103    lam = 0.001
104    lamMax = lam
105    nfev = 0
106    if Print:
107        print ' Hessian refinement on %d variables:'%(n)
108    Lam = np.zeros((n,n))
109    while icycle < maxcyc:
110        time0 = time.time()
111        M = func(x0,*args)
112        nfev += 1
113        chisq0 = np.sum(M**2)
114        Yvec,Amat = Hess(x0,*args)
115        Adiag = np.sqrt(np.diag(Amat))
116        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
117        if np.any(psing):                #hard singularity in matrix
118            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
119        Anorm = np.outer(Adiag,Adiag)
120        Yvec /= Adiag
121        Amat /= Anorm
122        while True:
123            Lam = np.eye(Amat.shape[0])*lam
124            Amatlam = Amat*(One+Lam)
125            try:
126                Xvec = nl.solve(Amatlam,Yvec)
127            except nl.LinAlgError:
128                print 'ouch #1'
129                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
130                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
131            Xvec /= Adiag
132            M2 = func(x0+Xvec,*args)
133            nfev += 1
134            chisq1 = np.sum(M2**2)
135            if chisq1 > chisq0:
136                lam *= 10.
137            else:
138                x0 += Xvec
139                lam /= 10.
140                break
141            if lam > 10.e5:
142                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
143                break
144        lamMax = max(lamMax,lam)
145        deltaChi2 = (chisq0-chisq1)/chisq0
146        if Print:
147            print ' Cycle: %d, Time: %.2fs, Chi**2: %.3g, Lambda: %.3g,  Delta: %.3g'%(
148                icycle,time.time()-time0,chisq1,lam,deltaChi2)
149        if deltaChi2 < ftol:
150            ifConverged = True
151            if Print: print "converged"
152            break
153        icycle += 1
154    M = func(x0,*args)
155    nfev += 1
156    Yvec,Amat = Hess(x0,*args)
157    Adiag = np.sqrt(np.diag(Amat))
158    Anorm = np.outer(Adiag,Adiag)
159    Lam = np.eye(Amat.shape[0])*lam
160    Amatlam = Amat/Anorm  #*(One+Lam)              #don't scale Amat to Marquardt array       
161    try:
162        Bmat = nl.inv(Amatlam)/Anorm  #*(One+Lam)      #don't rescale Bmat to Marquardt array
163        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
164    except nl.LinAlgError:
165        print 'ouch #2 linear algebra error in making v-cov matrix'
166        psing = []
167        if maxcyc:
168            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
169        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
170   
171def getVCov(varyNames,varyList,covMatrix):
172    '''obtain variance-covariance terms for a set of variables. NB: the varyList
173    and covMatrix were saved by the last least squares refinement so they must match.
174   
175    :param list varyNames: variable names to find v-cov matric for
176    :param list varyList: full list of all variables in v-cov matrix
177    :param nparray covMatrix: full variance-covariance matrix from the last
178     least squares refinement
179   
180    :returns: nparray vcov: variance-covariance matrix for the variables given
181     in varyNames
182   
183    '''
184    vcov = np.zeros((len(varyNames),len(varyNames)))
185    for i1,name1 in enumerate(varyNames):
186        for i2,name2 in enumerate(varyNames):
187            try:
188                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
189            except ValueError:
190                vcov[i1][i2] = 0.0
191#                if i1 == i2:
192#                    vcov[i1][i2] = 1e-20
193#                else:
194#                    vcov[i1][i2] = 0.0
195    return vcov
196
197def FindAtomIndexByIDs(atomData,IDs,Draw=True):
198    '''finds the set of atom array indices for a list of atom IDs. Will search
199    either the Atom table or the drawAtom table.
200   
201    :param list atomData: Atom or drawAtom table containting coordinates, etc.
202    :param list IDs: atom IDs to be found
203    :param bool Draw: True if drawAtom table to be searched; False if Atom table
204      is searched
205   
206    :returns: list indx: atom (or drawAtom) indices
207   
208    '''
209    indx = []
210    for i,atom in enumerate(atomData):
211        if Draw and atom[-3] in IDs:
212            indx.append(i)
213        elif atom[-1] in IDs:
214            indx.append(i)
215    return indx
216
217def FillAtomLookUp(atomData):
218    '''create a dictionary of atom indexes with atom IDs as keys
219   
220    :param list atomData: Atom table to be used
221   
222    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
223   
224    '''
225    atomLookUp = {}
226    for iatm,atom in enumerate(atomData):
227        atomLookUp[atom[-1]] = iatm
228    return atomLookUp
229
230def GetAtomsById(atomData,atomLookUp,IdList):
231    '''gets a list of atoms from Atom table that match a set of atom IDs
232   
233    :param list atomData: Atom table to be used
234    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
235    :param list IdList: atom IDs to be found
236   
237    :returns: list atoms: list of atoms found
238   
239    '''
240    atoms = []
241    for id in IdList:
242        atoms.append(atomData[atomLookUp[id]])
243    return atoms
244   
245def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
246    '''gets atom parameters for atoms using atom IDs
247   
248    :param list atomData: Atom table to be used
249    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
250    :param list IdList: atom IDs to be found
251    :param int itemLoc: pointer to desired 1st item in an atom table entry
252    :param int numItems: number of items to be retrieved
253   
254    :returns: type name: description
255   
256    '''
257    Items = []
258    if not isinstance(IdList,list):
259        IdList = [IdList,]
260    for id in IdList:
261        if numItems == 1:
262            Items.append(atomData[atomLookUp[id]][itemLoc])
263        else:
264            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
265    return Items
266   
267def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
268    '''default doc string
269   
270    :param type name: description
271   
272    :returns: type name: description
273   
274    '''
275    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
276    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
277    XYZ = []
278    for ind in indx:
279        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
280        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
281        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
282    return XYZ
283
284def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
285    '''default doc string
286   
287    :param type name: description
288   
289    :returns: type name: description
290   
291    '''
292    for atom in atomData:
293        XYZ = np.inner(Amat,atom[cx:cx+3])
294        if atom[cia] == 'A':
295            UIJ = atom[cia+2:cia+8]
296               
297def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
298    '''default doc string
299   
300    :param type name: description
301   
302    :returns: type name: description
303   
304    '''
305    TLStype,TLS = rbObj['ThermalMotion'][:2]
306    Tmat = np.zeros((3,3))
307    Lmat = np.zeros((3,3))
308    Smat = np.zeros((3,3))
309    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
310        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
311    if 'T' in TLStype:
312        Tmat = G2lat.U6toUij(TLS[:6])
313    if 'L' in TLStype:
314        Lmat = G2lat.U6toUij(TLS[6:12])
315    if 'S' in TLStype:
316        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
317    XYZ = np.inner(Amat,xyz)
318    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
319    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
320    beta = np.inner(np.inner(g,Umat),g)
321    return G2lat.UijtoU6(beta)*gvec   
322       
323def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
324    '''default doc string
325   
326    :param type name: description
327   
328    :returns: type name: description
329   
330    '''
331    cx,ct,cs,cia = atPtrs
332    TLStype,TLS = rbObj['ThermalMotion'][:2]
333    Tmat = np.zeros((3,3))
334    Lmat = np.zeros((3,3))
335    Smat = np.zeros((3,3))
336    G,g = G2lat.A2Gmat(Amat)
337    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
338    if 'T' in TLStype:
339        Tmat = G2lat.U6toUij(TLS[:6])
340    if 'L' in TLStype:
341        Lmat = G2lat.U6toUij(TLS[6:12])
342    if 'S' in TLStype:
343        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
344    for atom in atomData:
345        XYZ = np.inner(Amat,atom[cx:cx+3])
346        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
347        if 'U' in TSLtype:
348            atom[cia+1] = TLS[0]
349            atom[cia] = 'I'
350        else:
351            atom[cia] = 'A'
352            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
353            beta = np.inner(np.inner(g,Umat),g)
354            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
355
356def GetXYZDist(xyz,XYZ,Amat):
357    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
358        and are in crystal coordinates; Amat is crystal to Cart matrix
359   
360    :param type name: description
361   
362    :returns: type name: description
363   
364    '''
365    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
366
367def getAtomXYZ(atoms,cx):
368    '''default doc string
369   
370    :param type name: description
371   
372    :returns: type name: description
373   
374    '''
375    XYZ = []
376    for atom in atoms:
377        XYZ.append(atom[cx:cx+3])
378    return np.array(XYZ)
379
380def RotateRBXYZ(Bmat,Cart,oriQ):
381    '''rotate & transform cartesian coordinates to crystallographic ones
382    no translation applied. To be used for numerical derivatives
383   
384    :param type name: description
385   
386    :returns: type name: description
387   
388    '''
389    ''' returns crystal coordinates for atoms described by RBObj
390    '''
391    XYZ = np.zeros_like(Cart)
392    for i,xyz in enumerate(Cart):
393        XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
394    return XYZ
395
396def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
397    '''default doc string
398   
399    :param type name: description
400   
401    :returns: type name: description
402   
403    '''
404    ''' returns crystal coordinates for atoms described by RBObj
405    '''
406    RBRes = RBData[RBType][RBObj['RBId']]
407    if RBType == 'Vector':
408        vecs = RBRes['rbVect']
409        mags = RBRes['VectMag']
410        Cart = np.zeros_like(vecs[0])
411        for vec,mag in zip(vecs,mags):
412            Cart += vec*mag
413    elif RBType == 'Residue':
414        Cart = np.array(RBRes['rbXYZ'])
415        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
416            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
417            Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
418    XYZ = np.zeros_like(Cart)
419    for i,xyz in enumerate(Cart):
420        XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
421    return XYZ,Cart
422
423def UpdateMCSAxyz(Bmat,MCSA):
424    '''default doc string
425   
426    :param type name: description
427   
428    :returns: type name: description
429   
430    '''
431    xyz = []
432    atTypes = []
433    iatm = 0
434    for model in MCSA['Models'][1:]:        #skip the MD model
435        if model['Type'] == 'Atom':
436            xyz.append(model['Pos'][0])
437            atTypes.append(model['atType'])
438            iatm += 1
439        else:
440            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
441            Pos = np.array(model['Pos'][0])
442            Ori = np.array(model['Ori'][0])
443            Qori = AVdeg2Q(Ori[0],Ori[1:])
444            if model['Type'] == 'Vector':
445                vecs = RBRes['rbVect']
446                mags = RBRes['VectMag']
447                Cart = np.zeros_like(vecs[0])
448                for vec,mag in zip(vecs,mags):
449                    Cart += vec*mag
450            elif model['Type'] == 'Residue':
451                Cart = np.array(RBRes['rbXYZ'])
452                for itor,seq in enumerate(RBRes['rbSeq']):
453                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
454                    Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
455            if model['MolCent'][1]:
456                Cart -= model['MolCent'][0]
457            for i,x in enumerate(Cart):
458                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
459                atType = RBRes['rbTypes'][i]
460                atTypes.append(atType)
461                iatm += 1
462    return np.array(xyz),atTypes
463   
464def SetMolCent(model,RBData):
465    '''default doc string
466   
467    :param type name: description
468   
469    :returns: type name: description
470   
471    '''
472    rideList = []
473    RBRes = RBData[model['Type']][model['RBId']]
474    if model['Type'] == 'Vector':
475        vecs = RBRes['rbVect']
476        mags = RBRes['VectMag']
477        Cart = np.zeros_like(vecs[0])
478        for vec,mag in zip(vecs,mags):
479            Cart += vec*mag
480    elif model['Type'] == 'Residue':
481        Cart = np.array(RBRes['rbXYZ'])
482        for seq in RBRes['rbSeq']:
483            rideList += seq[3]
484    centList = set(range(len(Cart)))-set(rideList)
485    cent = np.zeros(3)
486    for i in centList:
487        cent += Cart[i]
488    model['MolCent'][0] = cent/len(centList)   
489   
490def UpdateRBUIJ(Bmat,Cart,RBObj):
491    '''default doc string
492   
493    :param type name: description
494   
495    :returns: type name: description
496   
497    '''
498    ''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
499    '''
500    TLStype,TLS = RBObj['ThermalMotion'][:2]
501    T = np.zeros(6)
502    L = np.zeros(6)
503    S = np.zeros(8)
504    if 'T' in TLStype:
505        T = TLS[:6]
506    if 'L' in TLStype:
507        L = np.array(TLS[6:12])*(np.pi/180.)**2
508    if 'S' in TLStype:
509        S = np.array(TLS[12:])*(np.pi/180.)
510    g = nl.inv(np.inner(Bmat,Bmat))
511    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
512        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
513    Uout = []
514    Q = RBObj['Orient'][0]
515    for X in Cart:
516        X = prodQVQ(Q,X)
517        if 'U' in TLStype:
518            Uout.append(['I',TLS[0],0,0,0,0,0,0])
519        elif not 'N' in TLStype:
520            U = [0,0,0,0,0,0]
521            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])
522            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])
523            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])
524            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]+  \
525                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
526            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]+  \
527                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
528            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]+  \
529                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
530            Umat = G2lat.U6toUij(U)
531            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
532            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
533        else:
534            Uout.append(['N',])
535    return Uout
536   
537def GetSHCoeff(pId,parmDict,SHkeys):
538    '''default doc string
539   
540    :param type name: description
541   
542    :returns: type name: description
543   
544    '''
545    SHCoeff = {}
546    for shkey in SHkeys:
547        shname = str(pId)+'::'+shkey
548        SHCoeff[shkey] = parmDict[shname]
549    return SHCoeff
550       
551def getMass(generalData):
552    '''Computes mass of unit cell contents
553   
554    :param dict generalData: The General dictionary in Phase
555   
556    :returns: float mass: Crystal unit cell mass in AMU.
557   
558    '''
559    mass = 0.
560    for i,elem in enumerate(generalData['AtomTypes']):
561        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
562    return mass   
563
564def getDensity(generalData):
565    '''calculate crystal structure density
566   
567    :param dict generalData: The General dictionary in Phase
568   
569    :returns: float density: crystal density in gm/cm^3
570   
571    '''
572    mass = getMass(generalData)
573    Volume = generalData['Cell'][7]
574    density = mass/(0.6022137*Volume)
575    return density,Volume/mass
576   
577def getWave(Parms):
578    '''returns wavelength from Instrument parameters dictionary
579   
580    :param dict Parms: Instrument parameters;
581        must contain:
582        Lam: single wavelength
583        or
584        Lam1: Ka1 radiation wavelength
585   
586    :returns: float wave: wavelength
587   
588    '''
589    try:
590        return Parms['Lam'][1]
591    except KeyError:
592        return Parms['Lam1'][1]
593       
594def El2Mass(Elements):
595    '''compute molecular weight from Elements
596   
597    :param dict Elements: elements in molecular formula;
598        each must contain
599        Num: number of atoms in formula
600        Mass: at. wt.
601   
602    :returns: float mass: molecular weight.
603   
604    '''
605    mass = 0
606    for El in Elements:
607        mass += Elements[El]['Num']*Elements[El]['Mass']
608    return mass
609       
610def Den2Vol(Elements,density):
611    '''converts density to molecular volume
612   
613    :param dict Elements: elements in molecular formula;
614        each must contain
615        Num: number of atoms in formula
616        Mass: at. wt.
617    :param float density: material density in gm/cm^3
618   
619    :returns: float volume: molecular volume in A^3
620   
621    '''
622    return El2Mass(Elements)/(density*0.6022137)
623   
624def Vol2Den(Elements,volume):
625    '''converts volume to density
626   
627    :param dict Elements: elements in molecular formula;
628        each must contain
629        Num: number of atoms in formula
630        Mass: at. wt.
631    :param float volume: molecular volume in A^3
632   
633    :returns: float density: material density in gm/cm^3
634   
635    '''
636    return El2Mass(Elements)/(volume*0.6022137)
637   
638def El2EstVol(Elements):
639    '''Estimate volume from molecular formula; assumes atom volume = 10A^3
640   
641    :param dict Elements: elements in molecular formula;
642        each must contain
643        Num: number of atoms in formula
644   
645    :returns: float volume: estimate of molecular volume in A^3
646   
647    '''
648    vol = 0
649    for El in Elements:
650        vol += 10.*Elements[El]['Num']
651    return vol
652   
653def XScattDen(Elements,vol,wave=0.):
654    '''Estimate X-ray scattering density from molecular formula & volume;
655    ignores valence, but includes anomalous effects
656   
657    :param dict Elements: elements in molecular formula;
658        each element must contain
659        Num: number of atoms in formula
660        Z: atomic number
661    :param float vol: molecular volume in A^3
662    :param float wave: optional wavelength in A
663   
664    :returns: float rho: scattering density in 10^10cm^-2;
665        if wave > 0 the includes f' contribution
666    :returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
667   
668    '''
669    rho = 0
670    mu = 0
671    if wave:
672        Xanom = XAnomAbs(Elements,wave)
673    for El in Elements:
674        f0 = Elements[El]['Z']
675        if wave:
676            f0 += Xanom[El][0]
677            mu += Xanom[El][2]*Elements[El]['Num']
678        rho += Elements[El]['Num']*f0
679    return 28.179*rho/vol,0.1*mu/vol
680   
681def wavekE(wavekE):
682    '''Convert wavelength to energy & vise versa
683   
684    :param float waveKe:wavelength in A or energy in kE
685   
686    :returns float waveKe:the other one
687   
688    '''
689    return 12.397639/wavekE
690       
691def XAnomAbs(Elements,wave):
692    kE = wavekE(wave)
693    Xanom = {}
694    for El in Elements:
695        Orbs = G2el.GetXsectionCoeff(El)
696        Xanom[El] = G2el.FPcalc(Orbs, kE)
697    return Xanom
698   
699   
700################################################################################
701##### distance, angle, planes, torsion stuff stuff
702################################################################################
703
704def getSyXYZ(XYZ,ops,SGData):
705    '''default doc string
706   
707    :param type name: description
708   
709    :returns: type name: description
710   
711    '''
712    XYZout = np.zeros_like(XYZ)
713    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
714        if op == '1':
715            XYZout[i] = xyz
716        else:
717            oprs = op.split('+')
718            unit = [0,0,0]
719            if len(oprs)>1:
720                unit = np.array(list(eval(oprs[1])))
721            syop =int(oprs[0])
722            inv = syop/abs(syop)
723            syop *= inv
724            cent = syop/100
725            syop %= 100
726            syop -= 1
727            M,T = SGData['SGOps'][syop]
728            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
729    return XYZout
730   
731def getRestDist(XYZ,Amat):
732    '''default doc string
733   
734    :param type name: description
735   
736    :returns: type name: description
737   
738    '''
739    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
740   
741def getRestDeriv(Func,XYZ,Amat,ops,SGData):
742    '''default doc string
743   
744    :param type name: description
745   
746    :returns: type name: description
747   
748    '''
749    deriv = np.zeros((len(XYZ),3))
750    dx = 0.00001
751    for j,xyz in enumerate(XYZ):
752        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
753            XYZ[j] -= x
754            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
755            XYZ[j] += 2*x
756            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
757            XYZ[j] -= x
758            deriv[j][i] = (d1-d2)/(2*dx)
759    return deriv.flatten()
760
761def getRestAngle(XYZ,Amat):
762    '''default doc string
763   
764    :param type name: description
765   
766    :returns: type name: description
767   
768    '''
769   
770    def calcVec(Ox,Tx,Amat):
771        return np.inner(Amat,(Tx-Ox))
772
773    VecA = calcVec(XYZ[1],XYZ[0],Amat)
774    VecA /= np.sqrt(np.sum(VecA**2))
775    VecB = calcVec(XYZ[1],XYZ[2],Amat)
776    VecB /= np.sqrt(np.sum(VecB**2))
777    edge = VecB-VecA
778    edge = np.sum(edge**2)
779    angle = (2.-edge)/2.
780    angle = max(angle,-1.)
781    return acosd(angle)
782   
783def getRestPlane(XYZ,Amat):
784    '''default doc string
785   
786    :param type name: description
787   
788    :returns: type name: description
789   
790    '''
791    sumXYZ = np.zeros(3)
792    for xyz in XYZ:
793        sumXYZ += xyz
794    sumXYZ /= len(XYZ)
795    XYZ = np.array(XYZ)-sumXYZ
796    XYZ = np.inner(Amat,XYZ).T
797    Zmat = np.zeros((3,3))
798    for i,xyz in enumerate(XYZ):
799        Zmat += np.outer(xyz.T,xyz)
800    Evec,Emat = nl.eig(Zmat)
801    Evec = np.sqrt(Evec)/(len(XYZ)-3)
802    Order = np.argsort(Evec)
803    return Evec[Order[0]]
804   
805def getRestChiral(XYZ,Amat):   
806    '''default doc string
807   
808    :param type name: description
809   
810    :returns: type name: description
811   
812    '''
813    VecA = np.empty((3,3))   
814    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
815    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
816    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
817    return nl.det(VecA)
818   
819def getRestTorsion(XYZ,Amat):
820    '''default doc string
821   
822    :param type name: description
823   
824    :returns: type name: description
825   
826    '''
827    VecA = np.empty((3,3))
828    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
829    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
830    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
831    D = nl.det(VecA)
832    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
833    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
834    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
835    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
836    Ang = 1.0
837    if abs(P12) < 1.0 and abs(P23) < 1.0:
838        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
839    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
840    return TOR
841   
842def calcTorsionEnergy(TOR,Coeff=[]):
843    '''default doc string
844   
845    :param type name: description
846   
847    :returns: type name: description
848   
849    '''
850    sum = 0.
851    if len(Coeff):
852        cof = np.reshape(Coeff,(3,3)).T
853        delt = TOR-cof[1]
854        delt = np.where(delt<-180.,delt+360.,delt)
855        delt = np.where(delt>180.,delt-360.,delt)
856        term = -cof[2]*delt**2
857        val = cof[0]*np.exp(term/1000.0)
858        pMax = cof[0][np.argmin(val)]
859        Eval = np.sum(val)
860        sum = Eval-pMax
861    return sum,Eval
862
863def getTorsionDeriv(XYZ,Amat,Coeff):
864    '''default doc string
865   
866    :param type name: description
867   
868    :returns: type name: description
869   
870    '''
871    deriv = np.zeros((len(XYZ),3))
872    dx = 0.00001
873    for j,xyz in enumerate(XYZ):
874        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
875            XYZ[j] -= x
876            tor = getRestTorsion(XYZ,Amat)
877            p1,d1 = calcTorsionEnergy(tor,Coeff)
878            XYZ[j] += 2*x
879            tor = getRestTorsion(XYZ,Amat)
880            p2,d2 = calcTorsionEnergy(tor,Coeff)           
881            XYZ[j] -= x
882            deriv[j][i] = (p2-p1)/(2*dx)
883    return deriv.flatten()
884
885def getRestRama(XYZ,Amat):
886    '''Computes a pair of torsion angles in a 5 atom string
887   
888    :param nparray XYZ: crystallographic coordinates of 5 atoms
889    :param nparray Amat: crystal to cartesian transformation matrix
890   
891    :returns: list (phi,psi) two torsion angles in degrees
892   
893    '''
894    phi = getRestTorsion(XYZ[:5],Amat)
895    psi = getRestTorsion(XYZ[1:],Amat)
896    return phi,psi
897   
898def calcRamaEnergy(phi,psi,Coeff=[]):
899    '''Computes pseudo potential energy from a pair of torsion angles and a
900    numerical description of the potential energy surface. Used to create
901    penalty function in LS refinement:     
902    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
903
904    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
905   
906    :param float phi: first torsion angle (:math:`\\phi`)
907    :param float psi: second torsion angle (:math:`\\psi`)
908    :param list Coeff: pseudo potential coefficients
909   
910    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
911      sum is used for penalty function.
912   
913    '''
914    sum = 0.
915    Eval = 0.
916    if len(Coeff):
917        cof = Coeff.T
918        dPhi = phi-cof[1]
919        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
920        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
921        dPsi = psi-cof[2]
922        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
923        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
924        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
925        val = cof[0]*np.exp(val/1000.)
926        pMax = cof[0][np.argmin(val)]
927        Eval = np.sum(val)
928        sum = Eval-pMax
929    return sum,Eval
930
931def getRamaDeriv(XYZ,Amat,Coeff):
932    '''Computes numerical derivatives of torsion angle pair pseudo potential
933    with respect of crystallographic atom coordinates of the 5 atom sequence
934   
935    :param nparray XYZ: crystallographic coordinates of 5 atoms
936    :param nparray Amat: crystal to cartesian transformation matrix
937    :param list Coeff: pseudo potential coefficients
938   
939    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
940     crystallographic xyz coordinates.
941   
942    '''
943    deriv = np.zeros((len(XYZ),3))
944    dx = 0.00001
945    for j,xyz in enumerate(XYZ):
946        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
947            XYZ[j] -= x
948            phi,psi = getRestRama(XYZ,Amat)
949            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
950            XYZ[j] += 2*x
951            phi,psi = getRestRama(XYZ,Amat)
952            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
953            XYZ[j] -= x
954            deriv[j][i] = (p2-p1)/(2*dx)
955    return deriv.flatten()
956
957def getRestPolefig(ODFln,SamSym,Grid):
958    '''default doc string
959   
960    :param type name: description
961   
962    :returns: type name: description
963   
964    '''
965    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
966    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
967    R = np.where(R <= 1.,2.*atand(R),0.0)
968    Z = np.zeros_like(R)
969    Z = G2lat.polfcal(ODFln,SamSym,R,P)
970    Z = np.reshape(Z,(Grid,Grid))
971    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
972
973def getRestPolefigDerv(HKL,Grid,SHCoeff):
974    '''default doc string
975   
976    :param type name: description
977   
978    :returns: type name: description
979   
980    '''
981    pass
982       
983def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
984    '''default doc string
985   
986    :param type name: description
987   
988    :returns: type name: description
989   
990    '''
991    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
992        TxT = inv*(np.inner(M,Tx)+T)+C+U
993        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
994       
995    inv = Top/abs(Top)
996    cent = abs(Top)/100
997    op = abs(Top)%100-1
998    M,T = SGData['SGOps'][op]
999    C = SGData['SGCen'][cent]
1000    dx = .00001
1001    deriv = np.zeros(6)
1002    for i in [0,1,2]:
1003        Oxyz[i] -= dx
1004        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1005        Oxyz[i] += 2*dx
1006        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1007        Oxyz[i] -= dx
1008        Txyz[i] -= dx
1009        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1010        Txyz[i] += 2*dx
1011        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1012        Txyz[i] -= dx
1013    return deriv
1014   
1015def getAngSig(VA,VB,Amat,SGData,covData={}):
1016    '''default doc string
1017   
1018    :param type name: description
1019   
1020    :returns: type name: description
1021   
1022    '''
1023    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
1024        TxT = inv*(np.inner(M,Tx)+T)+C
1025        TxT = G2spc.MoveToUnitCell(TxT)+U
1026        return np.inner(Amat,(TxT-Ox))
1027       
1028    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1029        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1030        VecA /= np.sqrt(np.sum(VecA**2))
1031        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1032        VecB /= np.sqrt(np.sum(VecB**2))
1033        edge = VecB-VecA
1034        edge = np.sum(edge**2)
1035        angle = (2.-edge)/2.
1036        angle = max(angle,-1.)
1037        return acosd(angle)
1038       
1039    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1040    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1041    invA = invB = 1
1042    invA = TopA/abs(TopA)
1043    invB = TopB/abs(TopB)
1044    centA = abs(TopA)/100
1045    centB = abs(TopB)/100
1046    opA = abs(TopA)%100-1
1047    opB = abs(TopB)%100-1
1048    MA,TA = SGData['SGOps'][opA]
1049    MB,TB = SGData['SGOps'][opB]
1050    CA = SGData['SGCen'][centA]
1051    CB = SGData['SGCen'][centB]
1052    if 'covMatrix' in covData:
1053        covMatrix = covData['covMatrix']
1054        varyList = covData['varyList']
1055        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1056        dx = .00001
1057        dadx = np.zeros(9)
1058        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1059        for i in [0,1,2]:
1060            OxA[i] -= dx
1061            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1062            OxA[i] += 2*dx
1063            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1064            OxA[i] -= dx
1065           
1066            TxA[i] -= dx
1067            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1068            TxA[i] += 2*dx
1069            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1070            TxA[i] -= dx
1071           
1072            TxB[i] -= dx
1073            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1074            TxB[i] += 2*dx
1075            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1076            TxB[i] -= dx
1077           
1078        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1079        if sigAng < 0.01:
1080            sigAng = 0.0
1081        return Ang,sigAng
1082    else:
1083        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1084
1085def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1086    '''default doc string
1087   
1088    :param type name: description
1089   
1090    :returns: type name: description
1091   
1092    '''
1093    def calcDist(Atoms,SyOps,Amat):
1094        XYZ = []
1095        for i,atom in enumerate(Atoms):
1096            Inv,M,T,C,U = SyOps[i]
1097            XYZ.append(np.array(atom[1:4]))
1098            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1099            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1100        V1 = XYZ[1]-XYZ[0]
1101        return np.sqrt(np.sum(V1**2))
1102       
1103    Inv = []
1104    SyOps = []
1105    names = []
1106    for i,atom in enumerate(Oatoms):
1107        names += atom[-1]
1108        Op,unit = Atoms[i][-1]
1109        inv = Op/abs(Op)
1110        m,t = SGData['SGOps'][abs(Op)%100-1]
1111        c = SGData['SGCen'][abs(Op)/100]
1112        SyOps.append([inv,m,t,c,unit])
1113    Dist = calcDist(Oatoms,SyOps,Amat)
1114   
1115    sig = -0.001
1116    if 'covMatrix' in covData:
1117        parmNames = []
1118        dx = .00001
1119        dadx = np.zeros(6)
1120        for i in range(6):
1121            ia = i/3
1122            ix = i%3
1123            Oatoms[ia][ix+1] += dx
1124            a0 = calcDist(Oatoms,SyOps,Amat)
1125            Oatoms[ia][ix+1] -= 2*dx
1126            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1127        covMatrix = covData['covMatrix']
1128        varyList = covData['varyList']
1129        DistVcov = getVCov(names,varyList,covMatrix)
1130        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1131        if sig < 0.001:
1132            sig = -0.001
1133   
1134    return Dist,sig
1135
1136def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1137    '''default doc string
1138   
1139    :param type name: description
1140   
1141    :returns: type name: description
1142   
1143    '''
1144
1145    def calcAngle(Atoms,SyOps,Amat):
1146        XYZ = []
1147        for i,atom in enumerate(Atoms):
1148            Inv,M,T,C,U = SyOps[i]
1149            XYZ.append(np.array(atom[1:4]))
1150            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1151            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1152        V1 = XYZ[1]-XYZ[0]
1153        V1 /= np.sqrt(np.sum(V1**2))
1154        V2 = XYZ[1]-XYZ[2]
1155        V2 /= np.sqrt(np.sum(V2**2))
1156        V3 = V2-V1
1157        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1158        return acosd(cang)
1159
1160    Inv = []
1161    SyOps = []
1162    names = []
1163    for i,atom in enumerate(Oatoms):
1164        names += atom[-1]
1165        Op,unit = Atoms[i][-1]
1166        inv = Op/abs(Op)
1167        m,t = SGData['SGOps'][abs(Op)%100-1]
1168        c = SGData['SGCen'][abs(Op)/100]
1169        SyOps.append([inv,m,t,c,unit])
1170    Angle = calcAngle(Oatoms,SyOps,Amat)
1171   
1172    sig = -0.01
1173    if 'covMatrix' in covData:
1174        parmNames = []
1175        dx = .00001
1176        dadx = np.zeros(9)
1177        for i in range(9):
1178            ia = i/3
1179            ix = i%3
1180            Oatoms[ia][ix+1] += dx
1181            a0 = calcAngle(Oatoms,SyOps,Amat)
1182            Oatoms[ia][ix+1] -= 2*dx
1183            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1184        covMatrix = covData['covMatrix']
1185        varyList = covData['varyList']
1186        AngVcov = getVCov(names,varyList,covMatrix)
1187        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1188        if sig < 0.01:
1189            sig = -0.01
1190   
1191    return Angle,sig
1192
1193def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1194    '''default doc string
1195   
1196    :param type name: description
1197   
1198    :returns: type name: description
1199   
1200    '''
1201
1202    def calcTorsion(Atoms,SyOps,Amat):
1203       
1204        XYZ = []
1205        for i,atom in enumerate(Atoms):
1206            Inv,M,T,C,U = SyOps[i]
1207            XYZ.append(np.array(atom[1:4]))
1208            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1209            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1210        V1 = XYZ[1]-XYZ[0]
1211        V2 = XYZ[2]-XYZ[1]
1212        V3 = XYZ[3]-XYZ[2]
1213        V1 /= np.sqrt(np.sum(V1**2))
1214        V2 /= np.sqrt(np.sum(V2**2))
1215        V3 /= np.sqrt(np.sum(V3**2))
1216        M = np.array([V1,V2,V3])
1217        D = nl.det(M)
1218        Ang = 1.0
1219        P12 = np.dot(V1,V2)
1220        P13 = np.dot(V1,V3)
1221        P23 = np.dot(V2,V3)
1222        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1223        return Tors
1224           
1225    Inv = []
1226    SyOps = []
1227    names = []
1228    for i,atom in enumerate(Oatoms):
1229        names += atom[-1]
1230        Op,unit = Atoms[i][-1]
1231        inv = Op/abs(Op)
1232        m,t = SGData['SGOps'][abs(Op)%100-1]
1233        c = SGData['SGCen'][abs(Op)/100]
1234        SyOps.append([inv,m,t,c,unit])
1235    Tors = calcTorsion(Oatoms,SyOps,Amat)
1236   
1237    sig = -0.01
1238    if 'covMatrix' in covData:
1239        parmNames = []
1240        dx = .00001
1241        dadx = np.zeros(12)
1242        for i in range(12):
1243            ia = i/3
1244            ix = i%3
1245            Oatoms[ia][ix+1] -= dx
1246            a0 = calcTorsion(Oatoms,SyOps,Amat)
1247            Oatoms[ia][ix+1] += 2*dx
1248            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1249            Oatoms[ia][ix+1] -= dx           
1250        covMatrix = covData['covMatrix']
1251        varyList = covData['varyList']
1252        TorVcov = getVCov(names,varyList,covMatrix)
1253        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1254        if sig < 0.01:
1255            sig = -0.01
1256   
1257    return Tors,sig
1258       
1259def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1260    '''default doc string
1261   
1262    :param type name: description
1263   
1264    :returns: type name: description
1265   
1266    '''
1267
1268    def calcDist(Atoms,SyOps,Amat):
1269        XYZ = []
1270        for i,atom in enumerate(Atoms):
1271            Inv,M,T,C,U = SyOps[i]
1272            XYZ.append(np.array(atom[1:4]))
1273            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1274            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1275        V1 = XYZ[1]-XYZ[0]
1276        return np.sqrt(np.sum(V1**2))
1277       
1278    def calcAngle(Atoms,SyOps,Amat):
1279        XYZ = []
1280        for i,atom in enumerate(Atoms):
1281            Inv,M,T,C,U = SyOps[i]
1282            XYZ.append(np.array(atom[1:4]))
1283            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1284            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1285        V1 = XYZ[1]-XYZ[0]
1286        V1 /= np.sqrt(np.sum(V1**2))
1287        V2 = XYZ[1]-XYZ[2]
1288        V2 /= np.sqrt(np.sum(V2**2))
1289        V3 = V2-V1
1290        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1291        return acosd(cang)
1292
1293    def calcTorsion(Atoms,SyOps,Amat):
1294       
1295        XYZ = []
1296        for i,atom in enumerate(Atoms):
1297            Inv,M,T,C,U = SyOps[i]
1298            XYZ.append(np.array(atom[1:4]))
1299            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1300            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1301        V1 = XYZ[1]-XYZ[0]
1302        V2 = XYZ[2]-XYZ[1]
1303        V3 = XYZ[3]-XYZ[2]
1304        V1 /= np.sqrt(np.sum(V1**2))
1305        V2 /= np.sqrt(np.sum(V2**2))
1306        V3 /= np.sqrt(np.sum(V3**2))
1307        M = np.array([V1,V2,V3])
1308        D = nl.det(M)
1309        Ang = 1.0
1310        P12 = np.dot(V1,V2)
1311        P13 = np.dot(V1,V3)
1312        P23 = np.dot(V2,V3)
1313        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1314        return Tors
1315           
1316    Inv = []
1317    SyOps = []
1318    names = []
1319    for i,atom in enumerate(Oatoms):
1320        names += atom[-1]
1321        Op,unit = Atoms[i][-1]
1322        inv = Op/abs(Op)
1323        m,t = SGData['SGOps'][abs(Op)%100-1]
1324        c = SGData['SGCen'][abs(Op)/100]
1325        SyOps.append([inv,m,t,c,unit])
1326    M = len(Oatoms)
1327    if M == 2:
1328        Val = calcDist(Oatoms,SyOps,Amat)
1329    elif M == 3:
1330        Val = calcAngle(Oatoms,SyOps,Amat)
1331    else:
1332        Val = calcTorsion(Oatoms,SyOps,Amat)
1333   
1334    sigVals = [-0.001,-0.01,-0.01]
1335    sig = sigVals[M-3]
1336    if 'covMatrix' in covData:
1337        parmNames = []
1338        dx = .00001
1339        N = M*3
1340        dadx = np.zeros(N)
1341        for i in range(N):
1342            ia = i/3
1343            ix = i%3
1344            Oatoms[ia][ix+1] += dx
1345            if M == 2:
1346                a0 = calcDist(Oatoms,SyOps,Amat)
1347            elif M == 3:
1348                a0 = calcAngle(Oatoms,SyOps,Amat)
1349            else:
1350                a0 = calcTorsion(Oatoms,SyOps,Amat)
1351            Oatoms[ia][ix+1] -= 2*dx
1352            if M == 2:
1353                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1354            elif M == 3:
1355                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1356            else:
1357                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1358        covMatrix = covData['covMatrix']
1359        varyList = covData['varyList']
1360        Vcov = getVCov(names,varyList,covMatrix)
1361        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1362        if sig < sigVals[M-3]:
1363            sig = sigVals[M-3]
1364   
1365    return Val,sig
1366       
1367def ValEsd(value,esd=0,nTZ=False):
1368    '''Format a floating point number with a given level of precision or
1369    with in crystallographic format with a "esd", as value(esd). If esd is
1370    negative the number is formatted with the level of significant figures
1371    appropriate if abs(esd) were the esd, but the esd is not included.
1372    if the esd is zero, approximately 6 significant figures are printed.
1373    nTZ=True causes "extra" zeros to be removed after the decimal place.
1374    for example:
1375
1376      * "1.235(3)" for value=1.2346 & esd=0.003
1377      * "1.235(3)e4" for value=12346. & esd=30
1378      * "1.235(3)e6" for value=0.12346e7 & esd=3000
1379      * "1.235" for value=1.2346 & esd=-0.003
1380      * "1.240" for value=1.2395 & esd=-0.003
1381      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1382      * "1.23460" for value=1.2346 & esd=0.0
1383
1384    :param float value: number to be formatted
1385    :param float esd: uncertainty or if esd < 0, specifies level of
1386      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1387    :param bool nTZ: True to remove trailing zeros (default is False)
1388    :returns: value(esd) or value as a string
1389
1390    '''
1391    # Note: this routine is Python 3 compatible -- I think
1392    if math.isnan(value): # invalid value, bail out
1393        return '?'
1394    if math.isnan(esd): # invalid esd, treat as zero
1395        esd = 0
1396        esdoff = 5
1397    elif esd != 0:
1398        # transform the esd to a one or two digit integer
1399        l = math.log10(abs(esd)) % 1
1400        # cut off of 19 1.9==>(19) but 1.95==>(2) N.B. log10(1.95) = 0.2900...
1401        if l < 0.290034611362518: l+= 1.       
1402        intesd = int(round(10**l)) # esd as integer
1403        # determine the number of digits offset for the esd
1404        esdoff = int(round(math.log10(intesd*1./abs(esd))))
1405    else:
1406        esdoff = 5
1407    valoff = 0
1408    if abs(value) < abs(esdoff): # value is effectively zero
1409        pass
1410    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1411        # where the digit offset is to the left of the decimal place or where too many
1412        # digits are needed
1413        if abs(value) > 1:
1414            valoff = int(math.log10(abs(value)))
1415        elif abs(value) > 0:
1416            valoff = int(math.log10(abs(value))-0.9999999)
1417        else:
1418            valoff = 0
1419    if esd != 0:
1420        if valoff+esdoff < 0:
1421            valoff = esdoff = 0
1422        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1423    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1424        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1425    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1426        if abs(value) > 0:           
1427            extra = -math.log10(abs(value))
1428        else:
1429            extra = 0
1430        if extra > 0: extra += 1
1431        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1432    if esd > 0:
1433        out += ("({:d})").format(intesd)  # add the esd
1434    elif nTZ and '.' in out:
1435        out = out.rstrip('0')  # strip zeros to right of decimal
1436        out = out.rstrip('.')  # and decimal place when not needed
1437    if valoff != 0:
1438        out += ("e{:d}").format(valoff) # add an exponent, when needed
1439    return out
1440
1441################################################################################
1442##### Fourier & charge flip stuff
1443################################################################################
1444
1445def adjHKLmax(SGData,Hmax):
1446    '''default doc string
1447   
1448    :param type name: description
1449   
1450    :returns: type name: description
1451   
1452    '''
1453    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1454        Hmax[0] = ((Hmax[0]+3)/6)*6
1455        Hmax[1] = ((Hmax[1]+3)/6)*6
1456        Hmax[2] = ((Hmax[2]+1)/4)*4
1457    else:
1458        Hmax[0] = ((Hmax[0]+2)/4)*4
1459        Hmax[1] = ((Hmax[1]+2)/4)*4
1460        Hmax[2] = ((Hmax[2]+2)/4)*4
1461
1462def OmitMap(data,reflDict,pgbar=None):
1463    '''default doc string
1464   
1465    :param type name: description
1466   
1467    :returns: type name: description
1468   
1469    '''
1470    generalData = data['General']
1471    if not generalData['Map']['MapType']:
1472        print '**** ERROR - Fourier map not defined'
1473        return
1474    mapData = generalData['Map']
1475    dmin = mapData['Resolution']
1476    SGData = generalData['SGData']
1477    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1478    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1479    cell = generalData['Cell'][1:8]       
1480    A = G2lat.cell2A(cell[:6])
1481    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1482    adjHKLmax(SGData,Hmax)
1483    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1484    time0 = time.time()
1485    for iref,ref in enumerate(reflDict['RefList']):
1486        if ref[4] >= dmin:
1487            Fosq,Fcsq,ph = ref[8:11]
1488            Uniq = np.inner(ref[:3],SGMT)
1489            Phi = np.inner(ref[:3],SGT)
1490            for i,hkl in enumerate(Uniq):        #uses uniq
1491                hkl = np.asarray(hkl,dtype='i')
1492                dp = 360.*Phi[i]                #and phi
1493                a = cosd(ph+dp)
1494                b = sind(ph+dp)
1495                phasep = complex(a,b)
1496                phasem = complex(a,-b)
1497                Fo = np.sqrt(Fosq)
1498                if '2Fo-Fc' in mapData['MapType']:
1499                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1500                else:
1501                    F = np.sqrt(Fosq)
1502                h,k,l = hkl+Hmax
1503                Fhkl[h,k,l] = F*phasep
1504                h,k,l = -hkl+Hmax
1505                Fhkl[h,k,l] = F*phasem
1506    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1507    M = np.mgrid[0:4,0:4,0:4]
1508    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
1509    iBeg = blkIds*rho0.shape/4
1510    iFin = (blkIds+1)*rho0.shape/4
1511    rho_omit = np.zeros_like(rho0)
1512    nBlk = 0
1513    for iB,iF in zip(iBeg,iFin):
1514        rho1 = np.copy(rho0)
1515        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
1516        Fnew = fft.ifftshift(fft.ifftn(rho1))
1517        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
1518        phase = Fnew/np.absolute(Fnew)
1519        OFhkl = np.absolute(Fhkl)*phase
1520        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
1521        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]])
1522        nBlk += 1
1523        pgbar.Update(nBlk)
1524    mapData['rho'] = np.real(rho_omit)/cell[6]
1525    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1526    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1527    return mapData
1528
1529def FourierMap(data,reflDict):
1530    '''default doc string
1531   
1532    :param type name: description
1533   
1534    :returns: type name: description
1535   
1536    '''
1537    generalData = data['General']
1538    if not generalData['Map']['MapType']:
1539        print '**** ERROR - Fourier map not defined'
1540        return
1541    mapData = generalData['Map']
1542    dmin = mapData['Resolution']
1543    SGData = generalData['SGData']
1544    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1545    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1546    cell = generalData['Cell'][1:8]       
1547    A = G2lat.cell2A(cell[:6])
1548    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1549    adjHKLmax(SGData,Hmax)
1550    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1551#    Fhkl[0,0,0] = generalData['F000X']
1552    time0 = time.time()
1553    for iref,ref in enumerate(reflDict['RefList']):
1554        if ref[4] >= dmin:
1555            Fosq,Fcsq,ph = ref[8:11]
1556            Uniq = np.inner(ref[:3],SGMT)
1557            Phi = np.inner(ref[:3],SGT)
1558            for i,hkl in enumerate(Uniq):        #uses uniq
1559                hkl = np.asarray(hkl,dtype='i')
1560                dp = 360.*Phi[i]                #and phi
1561                a = cosd(ph+dp)
1562                b = sind(ph+dp)
1563                phasep = complex(a,b)
1564                phasem = complex(a,-b)
1565                if 'Fobs' in mapData['MapType']:
1566                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1567                    h,k,l = hkl+Hmax
1568                    Fhkl[h,k,l] = F*phasep
1569                    h,k,l = -hkl+Hmax
1570                    Fhkl[h,k,l] = F*phasem
1571                elif 'Fcalc' in mapData['MapType']:
1572                    F = np.sqrt(Fcsq)
1573                    h,k,l = hkl+Hmax
1574                    Fhkl[h,k,l] = F*phasep
1575                    h,k,l = -hkl+Hmax
1576                    Fhkl[h,k,l] = F*phasem
1577                elif 'delt-F' in mapData['MapType']:
1578                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1579                    h,k,l = hkl+Hmax
1580                    Fhkl[h,k,l] = dF*phasep
1581                    h,k,l = -hkl+Hmax
1582                    Fhkl[h,k,l] = dF*phasem
1583                elif '2*Fo-Fc' in mapData['MapType']:
1584                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1585                    h,k,l = hkl+Hmax
1586                    Fhkl[h,k,l] = F*phasep
1587                    h,k,l = -hkl+Hmax
1588                    Fhkl[h,k,l] = F*phasem
1589                elif 'Patterson' in mapData['MapType']:
1590                    h,k,l = hkl+Hmax
1591                    Fhkl[h,k,l] = complex(Fosq,0.)
1592                    h,k,l = -hkl+Hmax
1593                    Fhkl[h,k,l] = complex(Fosq,0.)
1594    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1595    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1596    mapData['rho'] = np.real(rho)
1597    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1598    return mapData
1599   
1600# map printing for testing purposes
1601def printRho(SGLaue,rho,rhoMax):                         
1602    '''default doc string
1603   
1604    :param type name: description
1605   
1606    :returns: type name: description
1607   
1608    '''
1609    dim = len(rho.shape)
1610    if dim == 2:
1611        ix,jy = rho.shape
1612        for j in range(jy):
1613            line = ''
1614            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1615                line += (jy-j)*'  '
1616            for i in range(ix):
1617                r = int(100*rho[i,j]/rhoMax)
1618                line += '%4d'%(r)
1619            print line+'\n'
1620    else:
1621        ix,jy,kz = rho.shape
1622        for k in range(kz):
1623            print 'k = ',k
1624            for j in range(jy):
1625                line = ''
1626                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1627                    line += (jy-j)*'  '
1628                for i in range(ix):
1629                    r = int(100*rho[i,j,k]/rhoMax)
1630                    line += '%4d'%(r)
1631                print line+'\n'
1632## keep this
1633               
1634def findOffset(SGData,A,Fhkl):   
1635    '''default doc string
1636   
1637    :param type name: description
1638   
1639    :returns: type name: description
1640   
1641    '''
1642    if SGData['SpGrp'] == 'P 1':
1643        return [0,0,0]   
1644    hklShape = Fhkl.shape
1645    hklHalf = np.array(hklShape)/2
1646    sortHKL = np.argsort(Fhkl.flatten())
1647    Fdict = {}
1648    for hkl in sortHKL:
1649        HKL = np.unravel_index(hkl,hklShape)
1650        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1651        if F == 0.:
1652            break
1653        Fdict['%.6f'%(np.absolute(F))] = hkl
1654    Flist = np.flipud(np.sort(Fdict.keys()))
1655    F = str(1.e6)
1656    i = 0
1657    DH = []
1658    Dphi = []
1659    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1660    for F in Flist:
1661        hkl = np.unravel_index(Fdict[F],hklShape)
1662        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1663        Uniq = np.array(Uniq,dtype='i')
1664        Phi = np.array(Phi)
1665        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1666        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1667        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1668        ang0 = np.angle(Fh0,deg=True)/360.
1669        for H,phi in zip(Uniq,Phi)[1:]:
1670            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1671            dH = H-hkl
1672            dang = ang-ang0
1673            if np.any(np.abs(dH)- Hmax > 0):    #keep low order DHs
1674                continue
1675            DH.append(dH)
1676            Dphi.append((dang+.5) % 1.0)
1677        if i > 20 or len(DH) > 30:
1678            break
1679        i += 1
1680    DH = np.array(DH)
1681    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1682    Dphi = np.array(Dphi)
1683    steps = np.array(hklShape)
1684    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1685    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1686    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1687    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1688    hist,bins = np.histogram(Mmap,bins=1000)
1689#    for i,item in enumerate(hist[:10]):
1690#        print item,bins[i]
1691    chisq = np.min(Mmap)
1692    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1693    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1694#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1695    return DX
1696   
1697def ChargeFlip(data,reflDict,pgbar):
1698    '''default doc string
1699   
1700    :param type name: description
1701   
1702    :returns: type name: description
1703   
1704    '''
1705    generalData = data['General']
1706    mapData = generalData['Map']
1707    flipData = generalData['Flip']
1708    FFtable = {}
1709    if 'None' not in flipData['Norm element']:
1710        normElem = flipData['Norm element'].upper()
1711        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1712        for ff in FFs:
1713            if ff['Symbol'] == normElem:
1714                FFtable.update(ff)
1715    dmin = flipData['Resolution']
1716    SGData = generalData['SGData']
1717    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1718    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1719    cell = generalData['Cell'][1:8]       
1720    A = G2lat.cell2A(cell[:6])
1721    Vol = cell[6]
1722    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1723    adjHKLmax(SGData,Hmax)
1724    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1725    time0 = time.time()
1726    for iref,ref in enumerate(reflDict['RefList']):
1727        dsp = ref[4]
1728        if dsp >= dmin:
1729            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1730            if FFtable:
1731                SQ = 0.25/dsp**2
1732                ff *= G2el.ScatFac(FFtable,SQ)[0]
1733            if ref[8] > 0.:         #use only +ve Fobs**2
1734                E = np.sqrt(ref[8])/ff
1735            else:
1736                E = 0.
1737            ph = ref[10]
1738            ph = rn.uniform(0.,360.)
1739            Uniq = np.inner(ref[:3],SGMT)
1740            Phi = np.inner(ref[:3],SGT)
1741            for i,hkl in enumerate(Uniq):        #uses uniq
1742                hkl = np.asarray(hkl,dtype='i')
1743                dp = 360.*Phi[i]                #and phi
1744                a = cosd(ph+dp)
1745                b = sind(ph+dp)
1746                phasep = complex(a,b)
1747                phasem = complex(a,-b)
1748                h,k,l = hkl+Hmax
1749                Ehkl[h,k,l] = E*phasep
1750                h,k,l = -hkl+Hmax       #Friedel pair refl.
1751                Ehkl[h,k,l] = E*phasem
1752#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1753    CEhkl = copy.copy(Ehkl)
1754    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1755    Emask = ma.getmask(MEhkl)
1756    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1757    Ncyc = 0
1758    old = np.seterr(all='raise')
1759    while True:       
1760        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1761        CEsig = np.std(CErho)
1762        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1763        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
1764        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1765        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1766        phase = CFhkl/np.absolute(CFhkl)
1767        CEhkl = np.absolute(Ehkl)*phase
1768        Ncyc += 1
1769        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1770        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1771        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1772        if Rcf < 5.:
1773            break
1774        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1775        if not GoOn or Ncyc > 10000:
1776            break
1777    np.seterr(**old)
1778    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1779    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1780    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1781    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
1782       
1783    mapData['Rcf'] = Rcf
1784    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1785    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1786    mapData['rollMap'] = [0,0,0]
1787    return mapData
1788   
1789def SearchMap(data):
1790    '''Does a search of a density map for peaks meeting the criterion of peak
1791    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
1792    mapData is data['General']['mapData']; the map is also in mapData.
1793
1794    :param data: the phase data structure
1795
1796    :returns: (peaks,mags,dzeros) where
1797
1798        * peaks : ndarray
1799          x,y,z positions of the peaks found in the map
1800        * mags : ndarray
1801          the magnitudes of the peaks
1802        * dzeros : ndarray
1803          the distance of the peaks from  the unit cell origin
1804
1805    '''       
1806    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1807   
1808    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1809   
1810#    def noDuplicate(xyz,peaks,Amat):
1811#        XYZ = np.inner(Amat,xyz)
1812#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1813#            print ' Peak',xyz,' <0.5A from another peak'
1814#            return False
1815#        return True
1816#                           
1817    def fixSpecialPos(xyz,SGData,Amat):
1818        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1819        X = []
1820        xyzs = [equiv[0] for equiv in equivs]
1821        for x in xyzs:
1822            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1823                X.append(x)
1824        if len(X) > 1:
1825            return np.average(X,axis=0)
1826        else:
1827            return xyz
1828       
1829    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1830        Mag,x0,y0,z0,sig = parms
1831        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1832#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1833        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
1834       
1835    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1836        Mag,x0,y0,z0,sig = parms
1837        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1838        return M
1839       
1840    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1841        Mag,x0,y0,z0,sig = parms
1842        dMdv = np.zeros(([5,]+list(rX.shape)))
1843        delt = .01
1844        for i in range(5):
1845            parms[i] -= delt
1846            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1847            parms[i] += 2.*delt
1848            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1849            parms[i] -= delt
1850            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1851        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1852        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1853        dMdv = np.reshape(dMdv,(5,rX.size))
1854        Hess = np.inner(dMdv,dMdv)
1855       
1856        return Vec,Hess
1857       
1858    generalData = data['General']
1859    phaseName = generalData['Name']
1860    SGData = generalData['SGData']
1861    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1862    drawingData = data['Drawing']
1863    peaks = []
1864    mags = []
1865    dzeros = []
1866    try:
1867        mapData = generalData['Map']
1868        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1869        rho = copy.copy(mapData['rho'])     #don't mess up original
1870        mapHalf = np.array(rho.shape)/2
1871        res = mapData['Resolution']
1872        incre = np.array(rho.shape,dtype=np.float)
1873        step = max(1.0,1./res)+1
1874        steps = np.array(3*[step,])
1875    except KeyError:
1876        print '**** ERROR - Fourier map not defined'
1877        return peaks,mags
1878    rhoMask = ma.array(rho,mask=(rho<contLevel))
1879    indices = (-1,0,1)
1880    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1881    for roll in rolls:
1882        if np.any(roll):
1883            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1884    indx = np.transpose(rhoMask.nonzero())
1885    peaks = indx/incre
1886    mags = rhoMask[rhoMask.nonzero()]
1887    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1888        rho = rollMap(rho,ind)
1889        rMM = mapHalf-steps
1890        rMP = mapHalf+steps+1
1891        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1892        peakInt = np.sum(rhoPeak)*res**3
1893        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1894        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1895        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1896            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1897        x1 = result[0]
1898        if not np.any(x1 < 0):
1899            mag = x1[0]
1900            peak = (np.array(x1[1:4])-ind)/incre
1901        peak = fixSpecialPos(peak,SGData,Amat)
1902        rho = rollMap(rho,-ind)       
1903    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1904    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1905   
1906def sortArray(data,pos,reverse=False):
1907    '''data is a list of items
1908    sort by pos in list; reverse if True
1909    '''
1910    T = []
1911    for i,M in enumerate(data):
1912        T.append((M[pos],i))
1913    D = dict(zip(T,data))
1914    T.sort()
1915    if reverse:
1916        T.reverse()
1917    X = []
1918    for key in T:
1919        X.append(D[key])
1920    return X
1921
1922def PeaksEquiv(data,Ind):
1923    '''Find the equivalent map peaks for those selected. Works on the
1924    contents of data['Map Peaks'].
1925
1926    :param data: the phase data structure
1927    :param list Ind: list of selected peak indices
1928    :returns: augmented list of peaks including those related by symmetry to the
1929      ones in Ind
1930
1931    '''       
1932    def Duplicate(xyz,peaks,Amat):
1933        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1934            return True
1935        return False
1936                           
1937    generalData = data['General']
1938    cell = generalData['Cell'][1:7]
1939    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1940    A = G2lat.cell2A(cell)
1941    SGData = generalData['SGData']
1942    mapPeaks = data['Map Peaks']
1943    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1944    Indx = {}
1945    for ind in Ind:
1946        xyz = np.array(mapPeaks[ind][1:4])
1947        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1948        for jnd,xyz in enumerate(XYZ):       
1949            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1950    Ind = []
1951    for ind in Indx:
1952        if Indx[ind]:
1953            Ind.append(ind)
1954    return Ind
1955               
1956def PeaksUnique(data,Ind):
1957    '''Finds the symmetry unique set of peaks from those selected. Works on the
1958    contents of data['Map Peaks'].
1959
1960    :param data: the phase data structure
1961    :param list Ind: list of selected peak indices
1962    :returns: the list of symmetry unique peaks from among those given in Ind
1963
1964    '''       
1965#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1966
1967    def noDuplicate(xyz,peaks,Amat):
1968        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=1.0) for peak in peaks]:
1969            return False
1970        return True
1971                           
1972    generalData = data['General']
1973    cell = generalData['Cell'][1:7]
1974    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1975    A = G2lat.cell2A(cell)
1976    SGData = generalData['SGData']
1977    mapPeaks = data['Map Peaks']
1978    Indx = {}
1979    XYZ = {}
1980    for ind in Ind:
1981        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1982        Indx[ind] = True
1983    for ind in Ind:
1984        if Indx[ind]:
1985            xyz = XYZ[ind]
1986            for jnd in Ind:
1987                if ind != jnd and Indx[jnd]:                       
1988                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1989                    xyzs = np.array([equiv[0] for equiv in Equiv])
1990                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1991    Ind = []
1992    for ind in Indx:
1993        if Indx[ind]:
1994            Ind.append(ind)
1995    return Ind
1996   
1997################################################################################
1998##### single peak fitting profile fxn stuff
1999################################################################################
2000
2001def getCWsig(ins,pos):
2002    '''get CW peak profile sigma
2003   
2004    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
2005      as values only
2006    :param float pos: 2-theta of peak
2007    :returns: float getCWsig: peak sigma
2008   
2009    '''
2010    tp = tand(pos/2.0)
2011    return ins['U']*tp**2+ins['V']*tp+ins['W']
2012   
2013def getCWsigDeriv(pos):
2014    '''get derivatives of CW peak profile sigma wrt U,V, & W
2015   
2016    :param float pos: 2-theta of peak
2017   
2018    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
2019   
2020    '''
2021    tp = tand(pos/2.0)
2022    return tp**2,tp,1.0
2023   
2024def getCWgam(ins,pos):
2025    '''get CW peak profile gamma
2026   
2027    :param dict ins: instrument parameters with at least 'X' & 'Y'
2028      as values only
2029    :param float pos: 2-theta of peak
2030    :returns: float getCWgam: peak gamma
2031   
2032    '''
2033    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2034   
2035def getCWgamDeriv(pos):
2036    '''get derivatives of CW peak profile gamma wrt X & Y
2037   
2038    :param float pos: 2-theta of peak
2039   
2040    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
2041   
2042    '''
2043    return 1./cosd(pos/2.0),tand(pos/2.0)
2044   
2045def getTOFsig(ins,dsp):
2046    '''get TOF peak profile sigma
2047   
2048    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
2049      as values only
2050    :param float dsp: d-spacing of peak
2051   
2052    :returns: float getTOFsig: peak sigma
2053   
2054    '''
2055    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
2056   
2057def getTOFsigDeriv(dsp):
2058    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
2059   
2060    :param float dsp: d-spacing of peak
2061   
2062    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
2063   
2064    '''
2065    return 1.0,dsp**2,dsp**4,1./dsp**2
2066   
2067def getTOFgamma(ins,dsp):
2068    '''get TOF peak profile gamma
2069   
2070    :param dict ins: instrument parameters with at least 'X' & 'Y'
2071      as values only
2072    :param float dsp: d-spacing of peak
2073   
2074    :returns: float getTOFgamma: peak gamma
2075   
2076    '''
2077    return ins['X']*dsp+ins['Y']*dsp**2
2078   
2079def getTOFgammaDeriv(dsp):
2080    '''get derivatives of TOF peak profile gamma wrt X & Y
2081   
2082    :param float dsp: d-spacing of peak
2083   
2084    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
2085   
2086    '''
2087    return dsp,dsp**2
2088   
2089def getTOFbeta(ins,dsp):
2090    '''get TOF peak profile beta
2091   
2092    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
2093      as values only
2094    :param float dsp: d-spacing of peak
2095   
2096    :returns: float getTOFbeta: peak beat
2097   
2098    '''
2099    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
2100   
2101def getTOFbetaDeriv(dsp):
2102    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
2103   
2104    :param float dsp: d-spacing of peak
2105   
2106    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
2107   
2108    '''
2109    return 1.0,1./dsp**4,1./dsp**2
2110   
2111def getTOFalpha(ins,dsp):
2112    '''get TOF peak profile alpha
2113   
2114    :param dict ins: instrument parameters with at least 'alpha'
2115      as values only
2116    :param float dsp: d-spacing of peak
2117   
2118    :returns: flaot getTOFalpha: peak alpha
2119   
2120    '''
2121    return ins['alpha']/dsp
2122   
2123def getTOFalphaDeriv(dsp):
2124    '''get derivatives of TOF peak profile beta wrt alpha
2125   
2126    :param float dsp: d-spacing of peak
2127   
2128    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
2129   
2130    '''
2131    return 1./dsp
2132   
2133def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2134    '''set starting peak parameters for single peak fits from plot selection or auto selection
2135   
2136    :param dict Parms: instrument parameters dictionary
2137    :param dict Parms2: table lookup for TOF profile coefficients
2138    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
2139    :param float mag: peak top magnitude from pick
2140    :param bool ifQ: True if pos in Q
2141    :param bool useFit: True if use fitted CW Parms values (not defaults)
2142   
2143    :returns: list XY: peak list entry:
2144        for CW: [pos,0,mag,1,sig,0,gam,0]
2145        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2146        NB: mag refinement set by default, all others off
2147   
2148    '''
2149    ind = 0
2150    if useFit:
2151        ind = 1
2152    ins = {}
2153    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2154        for x in ['U','V','W','X','Y']:
2155            ins[x] = Parms[x][ind]
2156        if ifQ:                              #qplot - convert back to 2-theta
2157            pos = 2.0*asind(pos*wave/(4*math.pi))
2158        sig = getCWsig(ins,pos)
2159        gam = getCWgam(ins,pos)           
2160        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2161    else:
2162        if ifQ:
2163            dsp = 2.*np.pi/pos
2164            pos = Parms['difC']*dsp
2165        else:
2166            dsp = pos/Parms['difC'][1]
2167        if 'Pdabc' in Parms2:
2168            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
2169                ins[x] = Parms[x][ind]
2170            Pdabc = Parms2['Pdabc'].T
2171            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2172            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2173        else:
2174            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
2175                ins[x] = Parms[x][ind]
2176            alp = getTOFalpha(ins,dsp)
2177            bet = getTOFbeta(ins,dsp)
2178        sig = getTOFsig(ins,dsp)
2179        gam = getTOFgamma(ins,dsp)
2180        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2181    return XY
2182   
2183################################################################################
2184##### MC/SA stuff
2185################################################################################
2186
2187#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2188# Original Author: Travis Oliphant 2002
2189# Bug-fixes in 2006 by Tim Leslie
2190
2191
2192import numpy
2193from numpy import asarray, tan, exp, ones, squeeze, sign, \
2194     all, log, sqrt, pi, shape, array, minimum, where
2195from numpy import random
2196
2197#__all__ = ['anneal']
2198
2199_double_min = numpy.finfo(float).min
2200_double_max = numpy.finfo(float).max
2201class base_schedule(object):
2202    def __init__(self):
2203        self.dwell = 20
2204        self.learn_rate = 0.5
2205        self.lower = -10
2206        self.upper = 10
2207        self.Ninit = 50
2208        self.accepted = 0
2209        self.tests = 0
2210        self.feval = 0
2211        self.k = 0
2212        self.T = None
2213
2214    def init(self, **options):
2215        self.__dict__.update(options)
2216        self.lower = asarray(self.lower)
2217        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2218        self.upper = asarray(self.upper)
2219        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2220        self.k = 0
2221        self.accepted = 0
2222        self.feval = 0
2223        self.tests = 0
2224
2225    def getstart_temp(self, best_state):
2226        """ Find a matching starting temperature and starting parameters vector
2227        i.e. find x0 such that func(x0) = T0.
2228
2229        Parameters
2230        ----------
2231        best_state : _state
2232            A _state object to store the function value and x0 found.
2233
2234        returns
2235        -------
2236        x0 : array
2237            The starting parameters vector.
2238        """
2239
2240        assert(not self.dims is None)
2241        lrange = self.lower
2242        urange = self.upper
2243        fmax = _double_min
2244        fmin = _double_max
2245        for _ in range(self.Ninit):
2246            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2247            fval = self.func(x0, *self.args)
2248            self.feval += 1
2249            if fval > fmax:
2250                fmax = fval
2251            if fval < fmin:
2252                fmin = fval
2253                best_state.cost = fval
2254                best_state.x = array(x0)
2255
2256        self.T0 = (fmax-fmin)*1.5
2257        return best_state.x
2258       
2259    def set_range(self,x0,frac):
2260        delrange = frac*(self.upper-self.lower)
2261        self.upper = x0+delrange
2262        self.lower = x0-delrange
2263
2264    def accept_test(self, dE):
2265        T = self.T
2266        self.tests += 1
2267        if dE < 0:
2268            self.accepted += 1
2269            return 1
2270        p = exp(-dE*1.0/self.boltzmann/T)
2271        if (p > random.uniform(0.0, 1.0)):
2272            self.accepted += 1
2273            return 1
2274        return 0
2275
2276    def update_guess(self, x0):
2277        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2278
2279    def update_temp(self, x0):
2280        pass
2281
2282
2283#  A schedule due to Lester Ingber modified to use bounds - OK
2284class fast_sa(base_schedule):
2285    def init(self, **options):
2286        self.__dict__.update(options)
2287        if self.m is None:
2288            self.m = 1.0
2289        if self.n is None:
2290            self.n = 1.0
2291        self.c = self.m * exp(-self.n * self.quench)
2292
2293    def update_guess(self, x0):
2294        x0 = asarray(x0)
2295        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2296        T = self.T
2297        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2298        xnew = xc*(self.upper - self.lower)+self.lower
2299        return xnew
2300#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2301#        xc = y*(self.upper - self.lower)
2302#        xnew = x0 + xc
2303#        return xnew
2304
2305    def update_temp(self):
2306        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2307        self.k += 1
2308        return
2309
2310class cauchy_sa(base_schedule):     #modified to use bounds - not good
2311    def update_guess(self, x0):
2312        x0 = asarray(x0)
2313        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2314        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2315        xnew = xc*(self.upper - self.lower)+self.lower
2316        return xnew
2317#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2318#        xc = self.learn_rate * self.T * tan(numbers)
2319#        xnew = x0 + xc
2320#        return xnew
2321
2322    def update_temp(self):
2323        self.T = self.T0/(1+self.k)
2324        self.k += 1
2325        return
2326
2327class boltzmann_sa(base_schedule):
2328#    def update_guess(self, x0):
2329#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2330#        x0 = asarray(x0)
2331#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2332#
2333#        xnew = x0 + xc*std*self.learn_rate
2334#        return xnew
2335
2336    def update_temp(self):
2337        self.k += 1
2338        self.T = self.T0 / log(self.k+1.0)
2339        return
2340
2341class log_sa(base_schedule):        #OK
2342
2343    def init(self,**options):
2344        self.__dict__.update(options)
2345       
2346    def update_guess(self,x0):     #same as default
2347        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2348       
2349    def update_temp(self):
2350        self.k += 1
2351        self.T = self.T0*self.slope**self.k
2352       
2353class _state(object):
2354    def __init__(self):
2355        self.x = None
2356        self.cost = None
2357
2358# TODO:
2359#     allow for general annealing temperature profile
2360#     in that case use update given by alpha and omega and
2361#     variation of all previous updates and temperature?
2362
2363# Simulated annealing
2364
2365def anneal(func, x0, args=(), schedule='fast', full_output=0,
2366           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2367           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2368           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2369           ranRange=0.10,autoRan=False,dlg=None):
2370    """Minimize a function using simulated annealing.
2371
2372    Schedule is a schedule class implementing the annealing schedule.
2373    Available ones are 'fast', 'cauchy', 'boltzmann'
2374
2375    :param callable func: f(x, \*args)
2376        Function to be optimized.
2377    :param ndarray x0:
2378        Initial guess.
2379    :param tuple args:
2380        Extra parameters to `func`.
2381    :param base_schedule schedule:
2382        Annealing schedule to use (a class).
2383    :param bool full_output:
2384        Whether to return optional outputs.
2385    :param float T0:
2386        Initial Temperature (estimated as 1.2 times the largest
2387        cost-function deviation over random points in the range).
2388    :param float Tf:
2389        Final goal temperature.
2390    :param int maxeval:
2391        Maximum function evaluations.
2392    :param int maxaccept:
2393        Maximum changes to accept.
2394    :param int maxiter:
2395        Maximum cooling iterations.
2396    :param float learn_rate:
2397        Scale constant for adjusting guesses.
2398    :param float boltzmann:
2399        Boltzmann constant in acceptance test
2400        (increase for less stringent test at each temperature).
2401    :param float feps:
2402        Stopping relative error tolerance for the function value in
2403        last four coolings.
2404    :param float quench,m,n:
2405        Parameters to alter fast_sa schedule.
2406    :param float/ndarray lower,upper:
2407        Lower and upper bounds on `x`.
2408    :param int dwell:
2409        The number of times to search the space at each temperature.
2410    :param float slope:
2411        Parameter for log schedule
2412    :param bool ranStart=False:
2413        True for set 10% of ranges about x
2414
2415    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2416
2417     * xmin (ndarray): Point giving smallest value found.
2418     * Jmin (float): Minimum value of function found.
2419     * T (float): Final temperature.
2420     * feval (int): Number of function evaluations.
2421     * iters (int): Number of cooling iterations.
2422     * accept (int): Number of tests accepted.
2423     * retval (int): Flag indicating stopping condition:
2424
2425              *  0: Points no longer changing
2426              *  1: Cooled to final temperature
2427              *  2: Maximum function evaluations
2428              *  3: Maximum cooling iterations reached
2429              *  4: Maximum accepted query locations reached
2430              *  5: Final point not the minimum amongst encountered points
2431
2432    *Notes*:
2433    Simulated annealing is a random algorithm which uses no derivative
2434    information from the function being optimized. In practice it has
2435    been more useful in discrete optimization than continuous
2436    optimization, as there are usually better algorithms for continuous
2437    optimization problems.
2438
2439    Some experimentation by trying the difference temperature
2440    schedules and altering their parameters is likely required to
2441    obtain good performance.
2442
2443    The randomness in the algorithm comes from random sampling in numpy.
2444    To obtain the same results you can call numpy.random.seed with the
2445    same seed immediately before calling scipy.optimize.anneal.
2446
2447    We give a brief description of how the three temperature schedules
2448    generate new points and vary their temperature. Temperatures are
2449    only updated with iterations in the outer loop. The inner loop is
2450    over xrange(dwell), and new points are generated for
2451    every iteration in the inner loop. (Though whether the proposed
2452    new points are accepted is probabilistic.)
2453
2454    For readability, let d denote the dimension of the inputs to func.
2455    Also, let x_old denote the previous state, and k denote the
2456    iteration number of the outer loop. All other variables not
2457    defined below are input variables to scipy.optimize.anneal itself.
2458
2459    In the 'fast' schedule the updates are ::
2460
2461        u ~ Uniform(0, 1, size=d)
2462        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2463        xc = y * (upper - lower)
2464        x_new = x_old + xc
2465
2466        c = n * exp(-n * quench)
2467        T_new = T0 * exp(-c * k**quench)
2468
2469
2470    In the 'cauchy' schedule the updates are ::
2471
2472        u ~ Uniform(-pi/2, pi/2, size=d)
2473        xc = learn_rate * T * tan(u)
2474        x_new = x_old + xc
2475
2476        T_new = T0 / (1+k)
2477
2478    In the 'boltzmann' schedule the updates are ::
2479
2480        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2481        y ~ Normal(0, std, size=d)
2482        x_new = x_old + learn_rate * y
2483
2484        T_new = T0 / log(1+k)
2485
2486    """
2487    x0 = asarray(x0)
2488    lower = asarray(lower)
2489    upper = asarray(upper)
2490
2491    schedule = eval(schedule+'_sa()')
2492    #   initialize the schedule
2493    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2494                  learn_rate=learn_rate, lower=lower, upper=upper,
2495                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2496
2497    current_state, last_state, best_state = _state(), _state(), _state()
2498    if ranStart:
2499        schedule.set_range(x0,ranRange)
2500    if T0 is None:
2501        x0 = schedule.getstart_temp(best_state)
2502    else:
2503        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2504        best_state.x = None
2505        best_state.cost = numpy.Inf
2506
2507    last_state.x = asarray(x0).copy()
2508    fval = func(x0,*args)
2509    schedule.feval += 1
2510    last_state.cost = fval
2511    if last_state.cost < best_state.cost:
2512        best_state.cost = fval
2513        best_state.x = asarray(x0).copy()
2514    schedule.T = schedule.T0
2515    fqueue = [100, 300, 500, 700]
2516    iters = 1
2517    keepGoing = True
2518    bestn = 0
2519    while keepGoing:
2520        retval = 0
2521        for n in xrange(dwell):
2522            current_state.x = schedule.update_guess(last_state.x)
2523            current_state.cost = func(current_state.x,*args)
2524            schedule.feval += 1
2525
2526            dE = current_state.cost - last_state.cost
2527            if schedule.accept_test(dE):
2528                last_state.x = current_state.x.copy()
2529                last_state.cost = current_state.cost
2530                if last_state.cost < best_state.cost:
2531                    best_state.x = last_state.x.copy()
2532                    best_state.cost = last_state.cost
2533                    bestn = n
2534                    if best_state.cost < 1.0 and autoRan:
2535                        schedule.set_range(x0,best_state.cost/2.)                       
2536        if dlg:
2537            GoOn = dlg.Update(min(100.,best_state.cost*100),
2538                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2539                    'Best trial:',bestn,  \
2540                    'MC/SA Residual:',best_state.cost*100,'%', \
2541                    ))[0]
2542            if not GoOn:
2543                best_state.x = last_state.x.copy()
2544                best_state.cost = last_state.cost
2545                retval = 5
2546        schedule.update_temp()
2547        iters += 1
2548        # Stopping conditions
2549        # 0) last saved values of f from each cooling step
2550        #     are all very similar (effectively cooled)
2551        # 1) Tf is set and we are below it
2552        # 2) maxeval is set and we are past it
2553        # 3) maxiter is set and we are past it
2554        # 4) maxaccept is set and we are past it
2555        # 5) user canceled run via progress bar
2556
2557        fqueue.append(squeeze(last_state.cost))
2558        fqueue.pop(0)
2559        af = asarray(fqueue)*1.0
2560        if retval == 5:
2561            print ' User terminated run; incomplete MC/SA'
2562            keepGoing = False
2563            break
2564        if all(abs((af-af[0])/af[0]) < feps):
2565            retval = 0
2566            if abs(af[-1]-best_state.cost) > feps*10:
2567                retval = 5
2568#                print "Warning: Cooled to %f at %s but this is not" \
2569#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2570#                      + " the smallest point found."
2571            break
2572        if (Tf is not None) and (schedule.T < Tf):
2573            retval = 1
2574            break
2575        if (maxeval is not None) and (schedule.feval > maxeval):
2576            retval = 2
2577            break
2578        if (iters > maxiter):
2579            print "Warning: Maximum number of iterations exceeded."
2580            retval = 3
2581            break
2582        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2583            retval = 4
2584            break
2585
2586    if full_output:
2587        return best_state.x, best_state.cost, schedule.T, \
2588               schedule.feval, iters, schedule.accepted, retval
2589    else:
2590        return best_state.x, retval
2591
2592def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
2593    outlist = []
2594    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
2595    for n in range(iCyc):
2596        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2597        outlist.append(result[0])
2598        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2599    out_q.put(outlist)
2600
2601def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2602    import multiprocessing as mp
2603   
2604    nprocs = mp.cpu_count()
2605    out_q = mp.Queue()
2606    procs = []
2607    iCyc = np.zeros(nprocs)
2608    for i in range(nCyc):
2609        iCyc[i%nprocs] += 1
2610    for i in range(nprocs):
2611        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
2612        procs.append(p)
2613        p.start()
2614    resultlist = []
2615    for i in range(nprocs):
2616        resultlist += out_q.get()
2617    for p in procs:
2618        p.join()
2619    return resultlist
2620
2621def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2622    '''default doc string
2623   
2624    :param type name: description
2625   
2626    :returns: type name: description
2627   
2628    '''
2629   
2630    twopi = 2.0*np.pi
2631    global tsum
2632    tsum = 0.
2633   
2634    def getMDparms(item,pfx,parmDict,varyList):
2635        parmDict[pfx+'MDaxis'] = item['axis']
2636        parmDict[pfx+'MDval'] = item['Coef'][0]
2637        if item['Coef'][1]:
2638            varyList += [pfx+'MDval',]
2639            limits = item['Coef'][2]
2640            lower.append(limits[0])
2641            upper.append(limits[1])
2642                       
2643    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2644        parmDict[pfx+'Atype'] = item['atType']
2645        aTypes |= set([item['atType'],]) 
2646        pstr = ['Ax','Ay','Az']
2647        XYZ = [0,0,0]
2648        for i in range(3):
2649            name = pfx+pstr[i]
2650            parmDict[name] = item['Pos'][0][i]
2651            XYZ[i] = parmDict[name]
2652            if item['Pos'][1][i]:
2653                varyList += [name,]
2654                limits = item['Pos'][2][i]
2655                lower.append(limits[0])
2656                upper.append(limits[1])
2657        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2658           
2659    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2660        parmDict[mfx+'MolCent'] = item['MolCent']
2661        parmDict[mfx+'RBId'] = item['RBId']
2662        pstr = ['Px','Py','Pz']
2663        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
2664        for i in range(3):
2665            name = mfx+pstr[i]
2666            parmDict[name] = item['Pos'][0][i]
2667            if item['Pos'][1][i]:
2668                varyList += [name,]
2669                limits = item['Pos'][2][i]
2670                lower.append(limits[0])
2671                upper.append(limits[1])
2672        AV = item['Ori'][0]
2673        A = AV[0]
2674        V = AV[1:]
2675        for i in range(4):
2676            name = mfx+ostr[i]
2677            if i:
2678                parmDict[name] = V[i-1]
2679            else:
2680                parmDict[name] = A
2681            if item['Ovar'] == 'AV':
2682                varyList += [name,]
2683                limits = item['Ori'][2][i]
2684                lower.append(limits[0])
2685                upper.append(limits[1])
2686            elif item['Ovar'] == 'A' and not i:
2687                varyList += [name,]
2688                limits = item['Ori'][2][i]
2689                lower.append(limits[0])
2690                upper.append(limits[1])
2691        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2692            for i in range(len(item['Tor'][0])):
2693                name = mfx+'Tor'+str(i)
2694                parmDict[name] = item['Tor'][0][i]
2695                if item['Tor'][1][i]:
2696                    varyList += [name,]
2697                    limits = item['Tor'][2][i]
2698                    lower.append(limits[0])
2699                    upper.append(limits[1])
2700        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2701        aTypes |= set(atypes)
2702        atNo += len(atypes)
2703        return atNo
2704               
2705    def GetAtomM(Xdata,SGData):
2706        Mdata = []
2707        for xyz in Xdata:
2708            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2709        return np.array(Mdata)
2710       
2711    def GetAtomT(RBdata,parmDict):
2712        'Needs a doc string'
2713        atNo = parmDict['atNo']
2714        nfixAt = parmDict['nfixAt']
2715        Tdata = atNo*[' ',]
2716        for iatm in range(nfixAt):
2717            parm = ':'+str(iatm)+':Atype'
2718            if parm in parmDict:
2719                Tdata[iatm] = aTypes.index(parmDict[parm])
2720        iatm = nfixAt
2721        for iObj in range(parmDict['nObj']):
2722            pfx = str(iObj)+':'
2723            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2724                if parmDict[pfx+'Type'] == 'Vector':
2725                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2726                    nAtm = len(RBRes['rbVect'][0])
2727                else:       #Residue
2728                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2729                    nAtm = len(RBRes['rbXYZ'])
2730                for i in range(nAtm):
2731                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2732                    iatm += 1
2733            elif parmDict[pfx+'Type'] == 'Atom':
2734                atNo = parmDict[pfx+'atNo']
2735                parm = pfx+'Atype'              #remove extra ':'
2736                if parm in parmDict:
2737                    Tdata[atNo] = aTypes.index(parmDict[parm])
2738                iatm += 1
2739            else:
2740                continue        #skips March Dollase
2741        return Tdata
2742       
2743    def GetAtomX(RBdata,parmDict):
2744        'Needs a doc string'
2745        Bmat = parmDict['Bmat']
2746        atNo = parmDict['atNo']
2747        nfixAt = parmDict['nfixAt']
2748        Xdata = np.zeros((3,atNo))
2749        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2750        for iatm in range(nfixAt):
2751            for key in keys:
2752                parm = ':'+str(iatm)+key
2753                if parm in parmDict:
2754                    keys[key][iatm] = parmDict[parm]
2755        iatm = nfixAt
2756        for iObj in range(parmDict['nObj']):
2757            pfx = str(iObj)+':'
2758            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2759                if parmDict[pfx+'Type'] == 'Vector':
2760                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2761                    vecs = RBRes['rbVect']
2762                    mags = RBRes['VectMag']
2763                    Cart = np.zeros_like(vecs[0])
2764                    for vec,mag in zip(vecs,mags):
2765                        Cart += vec*mag
2766                elif parmDict[pfx+'Type'] == 'Residue':
2767                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2768                    Cart = np.array(RBRes['rbXYZ'])
2769                    for itor,seq in enumerate(RBRes['rbSeq']):
2770                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
2771                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
2772                if parmDict[pfx+'MolCent'][1]:
2773                    Cart -= parmDict[pfx+'MolCent'][0]
2774                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2775                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2776                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
2777                iatm += len(Cart)
2778            elif parmDict[pfx+'Type'] == 'Atom':
2779                atNo = parmDict[pfx+'atNo']
2780                for key in keys:
2781                    parm = pfx+key[1:]              #remove extra ':'
2782                    if parm in parmDict:
2783                        keys[key][atNo] = parmDict[parm]
2784                iatm += 1
2785            else:
2786                continue        #skips March Dollase
2787        return Xdata.T
2788       
2789    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
2790        allX = np.inner(Xdata,SGM)+SGT
2791        allT = np.repeat(Tdata,allX.shape[1])
2792        allM = np.repeat(Mdata,allX.shape[1])
2793        allX = np.reshape(allX,(-1,3))
2794        return allT,allM,allX
2795
2796    def getAllX(Xdata,SGM,SGT):
2797        allX = np.inner(Xdata,SGM)+SGT
2798        allX = np.reshape(allX,(-1,3))
2799        return allX
2800       
2801    def normQuaternions(RBdata,parmDict,varyList,values):
2802        for iObj in range(parmDict['nObj']):
2803            pfx = str(iObj)+':'
2804            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2805                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2806                A,V = Q2AVdeg(Qori)
2807                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
2808                    if i:
2809                        parmDict[pfx+name] = V[i-1]
2810                    else:
2811                        parmDict[pfx+name] = A
2812       
2813    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
2814        ''' Compute structure factors for all h,k,l for phase
2815        input:
2816            refList: [ref] where each ref = h,k,l,m,d,...
2817            rcov:   array[nref,nref] covariance terms between Fo^2 values
2818            ifInv:  bool True if centrosymmetric
2819            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
2820            RBdata: [dict] rigid body dictionary
2821            varyList: [list] names of varied parameters in MC/SA (not used here)           
2822            ParmDict: [dict] problem parameters
2823        puts result F^2 in each ref[5] in refList
2824        returns:
2825            delt-F*rcov*delt-F/sum(Fo^2)^2
2826        '''       
2827        global tsum
2828        t0 = time.time()
2829        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
2830        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
2831        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
2832        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
2833        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
2834        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
2835        refList[5] = Aterm
2836        if not ifInv:
2837            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
2838        if len(cosTable):        #apply MD correction
2839            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
2840        sumFcsq = np.sum(refList[5])
2841        scale = parmDict['sumFosq']/sumFcsq
2842        refList[5] *= scale
2843        refList[6] = refList[4]-refList[5]
2844        M = np.inner(refList[6],np.inner(rcov,refList[6]))
2845        tsum += (time.time()-t0)
2846        return M/np.sum(refList[4]**2)
2847
2848    sq8ln2 = np.sqrt(8*np.log(2))
2849    sq2pi = np.sqrt(2*np.pi)
2850    sq4pi = np.sqrt(4*np.pi)
2851    generalData = data['General']
2852    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2853    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
2854    SGData = generalData['SGData']
2855    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
2856    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
2857    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
2858    fixAtoms = data['Atoms']                       #if any
2859    cx,ct,cs = generalData['AtomPtrs'][:3]
2860    aTypes = set([])
2861    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2862    varyList = []
2863    atNo = 0
2864    for atm in fixAtoms:
2865        pfx = ':'+str(atNo)+':'
2866        parmDict[pfx+'Atype'] = atm[ct]
2867        aTypes |= set([atm[ct],])
2868        pstr = ['Ax','Ay','Az']
2869        parmDict[pfx+'Amul'] = atm[cs+1]
2870        for i in range(3):
2871            name = pfx+pstr[i]
2872            parmDict[name] = atm[cx+i]
2873        atNo += 1
2874    parmDict['nfixAt'] = len(fixAtoms)       
2875    MCSA = generalData['MCSA controls']
2876    reflName = MCSA['Data source']
2877    phaseName = generalData['Name']
2878    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2879    upper = []
2880    lower = []
2881    MDvec = np.zeros(3)
2882    for i,item in enumerate(MCSAObjs):
2883        mfx = str(i)+':'
2884        parmDict[mfx+'Type'] = item['Type']
2885        if item['Type'] == 'MD':
2886            getMDparms(item,mfx,parmDict,varyList)
2887            MDvec = np.array(item['axis'])
2888        elif item['Type'] == 'Atom':
2889            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2890            parmDict[mfx+'atNo'] = atNo
2891            atNo += 1
2892        elif item['Type'] in ['Residue','Vector']:
2893            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2894    parmDict['atNo'] = atNo                 #total no. of atoms
2895    parmDict['nObj'] = len(MCSAObjs)
2896    aTypes = list(aTypes)
2897    Tdata = GetAtomT(RBdata,parmDict)
2898    Xdata = GetAtomX(RBdata,parmDict)
2899    Mdata = GetAtomM(Xdata,SGData)
2900    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
2901    FFtables = G2el.GetFFtable(aTypes)
2902    refs = []
2903    allFF = []
2904    cosTable = []
2905    sumFosq = 0
2906    if 'PWDR' in reflName:
2907        for ref in reflData:
2908            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2909            if d >= MCSA['dmin']:
2910                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2911                SQ = 0.25/d**2
2912                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2913                refs.append([h,k,l,m,f*m,pos,sig])
2914                sumFosq += f*m
2915                Heqv = np.inner(np.array([h,k,l]),SGMT)
2916                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
2917        nRef = len(refs)
2918        cosTable = np.array(cosTable)**2
2919        rcov = np.zeros((nRef,nRef))
2920        for iref,refI in enumerate(refs):
2921            rcov[iref][iref] = 1./(sq4pi*refI[6])
2922            for jref,refJ in enumerate(refs[:iref]):
2923                t1 = refI[6]**2+refJ[6]**2
2924                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2925                if t2 > 10.:
2926                    rcov[iref][jref] = 0.
2927                else:
2928                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2929        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2930        Rdiag = np.sqrt(np.diag(rcov))
2931        Rnorm = np.outer(Rdiag,Rdiag)
2932        rcov /= Rnorm
2933    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
2934        vNames = []
2935        pfx = str(data['pId'])+'::PWLref:'
2936        for iref,refI in enumerate(reflData):           #Pawley reflection set
2937            h,k,l,m,d,v,f,s = refI
2938            if d >= MCSA['dmin'] and v:       #skip unrefined ones
2939                vNames.append(pfx+str(iref))
2940                SQ = 0.25/d**2
2941                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2942                refs.append([h,k,l,m,f*m,iref,0.])
2943                sumFosq += f*m
2944                Heqv = np.inner(np.array([h,k,l]),SGMT)
2945                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
2946        cosTable = np.array(cosTable)**2
2947        nRef = len(refs)
2948#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
2949        if covData['freshCOV'] or  MCSA['newDmin']:
2950            vList = covData['varyList']
2951            covMatrix = covData['covMatrix']
2952            rcov = getVCov(vNames,vList,covMatrix)
2953            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2954            Rdiag = np.sqrt(np.diag(rcov))
2955            Rdiag = np.where(Rdiag,Rdiag,1.0)
2956            Rnorm = np.outer(Rdiag,Rdiag)
2957            rcov /= Rnorm
2958            MCSA['rcov'] = rcov
2959            covData['freshCOV'] = False
2960            MCSA['newDmin'] = False
2961        else:
2962            rcov = MCSA['rcov']
2963    elif 'HKLF' in reflName:
2964        for ref in reflData:
2965            [h,k,l,m,d],f = ref[:5],ref[6]
2966            if d >= MCSA['dmin']:
2967                SQ = 0.25/d**2
2968                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2969                refs.append([h,k,l,m,f*m,0.,0.])
2970                sumFosq += f*m
2971        nRef = len(refs)
2972        rcov = np.identity(len(refs))
2973    allFF = np.array(allFF).T
2974    refs = np.array(refs).T
2975    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2976    print ' Number of parameters varied: %d'%(len(varyList))
2977    parmDict['sumFosq'] = sumFosq
2978    x0 = [parmDict[val] for val in varyList]
2979    ifInv = SGData['SGInv']
2980    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
2981        schedule=MCSA['Algorithm'], full_output=True,
2982        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2983        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
2984        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2985        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
2986        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
2987    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
2988#    for ref in refs.T:
2989#        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])
2990#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
2991    Result = [False,False,results[1],results[2],]+list(results[0])
2992    Result.append(varyList)
2993    return Result,tsum
2994
2995       
2996################################################################################
2997##### Quaternion stuff
2998################################################################################
2999
3000def prodQQ(QA,QB):
3001    ''' Grassman quaternion product
3002        QA,QB quaternions; q=r+ai+bj+ck
3003    '''
3004    D = np.zeros(4)
3005    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
3006    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
3007    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
3008    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
3009   
3010#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
3011#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
3012   
3013    return D
3014   
3015def normQ(QA):
3016    ''' get length of quaternion & normalize it
3017        q=r+ai+bj+ck
3018    '''
3019    n = np.sqrt(np.sum(np.array(QA)**2))
3020    return QA/n
3021   
3022def invQ(Q):
3023    '''
3024        get inverse of quaternion
3025        q=r+ai+bj+ck; q* = r-ai-bj-ck
3026    '''
3027    return Q*np.array([1,-1,-1,-1])
3028   
3029def prodQVQ(Q,V):
3030    """
3031    compute the quaternion vector rotation qvq-1 = v'
3032    q=r+ai+bj+ck
3033    """
3034    T2 = Q[0]*Q[1]
3035    T3 = Q[0]*Q[2]
3036    T4 = Q[0]*Q[3]
3037    T5 = -Q[1]*Q[1]
3038    T6 = Q[1]*Q[2]
3039    T7 = Q[1]*Q[3]
3040    T8 = -Q[2]*Q[2]
3041    T9 = Q[2]*Q[3]
3042    T10 = -Q[3]*Q[3]
3043    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3044    VP = 2.*np.inner(V,M)
3045    return VP+V
3046   
3047def Q2Mat(Q):
3048    ''' make rotation matrix from quaternion
3049        q=r+ai+bj+ck
3050    '''
3051    QN = normQ(Q)
3052    aa = QN[0]**2
3053    ab = QN[0]*QN[1]
3054    ac = QN[0]*QN[2]
3055    ad = QN[0]*QN[3]
3056    bb = QN[1]**2
3057    bc = QN[1]*QN[2]
3058    bd = QN[1]*QN[3]
3059    cc = QN[2]**2
3060    cd = QN[2]*QN[3]
3061    dd = QN[3]**2
3062    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3063        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3064        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3065    return np.array(M)
3066   
3067def AV2Q(A,V):
3068    ''' convert angle (radians) & vector to quaternion
3069        q=r+ai+bj+ck
3070    '''
3071    Q = np.zeros(4)
3072    d = nl.norm(np.array(V))
3073    if d:
3074        V /= d
3075        if not A:       #==0.
3076            A = 2.*np.pi
3077        p = A/2.
3078        Q[0] = np.cos(p)
3079        Q[1:4] = V*np.sin(p)
3080    else:
3081        Q[3] = 1.
3082    return Q
3083   
3084def AVdeg2Q(A,V):
3085    ''' convert angle (degrees) & vector to quaternion
3086        q=r+ai+bj+ck
3087    '''
3088    Q = np.zeros(4)
3089    d = nl.norm(np.array(V))
3090    if not A:       #== 0.!
3091        A = 360.
3092    if d:
3093        V /= d
3094        p = A/2.
3095        Q[0] = cosd(p)
3096        Q[1:4] = V*sind(p)
3097    else:
3098        Q[3] = 1.
3099    return Q
3100   
3101def Q2AVdeg(Q):
3102    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3103        q=r+ai+bj+ck
3104    '''
3105    A = 2.*acosd(Q[0])
3106    V = np.array(Q[1:])
3107    V /= sind(A/2.)
3108    return A,V
3109   
3110def Q2AV(Q):
3111    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3112        q=r+ai+bj+ck
3113    '''
3114    A = 2.*np.arccos(Q[0])
3115    V = np.array(Q[1:])
3116    V /= np.sin(A/2.)
3117    return A,V
3118   
3119def randomQ(r0,r1,r2,r3):
3120    ''' create random quaternion from 4 random numbers in range (-1,1)
3121    '''
3122    sum = 0
3123    Q = np.array(4)
3124    Q[0] = r0
3125    sum += Q[0]**2
3126    Q[1] = np.sqrt(1.-sum)*r1
3127    sum += Q[1]**2
3128    Q[2] = np.sqrt(1.-sum)*r2
3129    sum += Q[2]**2
3130    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3131    return Q
3132   
3133def randomAVdeg(r0,r1,r2,r3):
3134    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3135    '''
3136    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3137   
3138def makeQuat(A,B,C):
3139    ''' Make quaternion from rotation of A vector to B vector about C axis
3140
3141        :param np.array A,B,C: Cartesian 3-vectors
3142        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3143    '''
3144
3145    V1 = np.cross(A,C)
3146    V2 = np.cross(B,C)
3147    if nl.norm(V1)*nl.norm(V2):
3148        V1 /= nl.norm(V1)
3149        V2 /= nl.norm(V2)
3150        V3 = np.cross(V1,V2)
3151    else:
3152        V3 = np.zeros(3)
3153    Q = np.array([0.,0.,0.,1.])
3154    D = 0.
3155    if nl.norm(V3):
3156        V3 /= nl.norm(V3)
3157        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3158        D = np.arccos(D1)/2.0
3159        V1 = C-V3
3160        V2 = C+V3
3161        DM = nl.norm(V1)
3162        DP = nl.norm(V2)
3163        S = np.sin(D)
3164        Q[0] = np.cos(D)
3165        Q[1:] = V3*S
3166        D *= 2.
3167        if DM > DP:
3168            D *= -1.
3169    return Q,D
3170   
3171def annealtests():
3172    from numpy import cos
3173    # minimum expected at ~-0.195
3174    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3175    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3176    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3177    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3178
3179    # minimum expected at ~[-0.195, -0.1]
3180    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3181    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')
3182    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')
3183    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')
3184
3185
3186if __name__ == '__main__':
3187    annealtests()
Note: See TracBrowser for help on using the repository browser.