source: trunk/GSASIImath.py @ 1513

Last change on this file since 1513 was 1513, checked in by toby, 8 years ago

update and rebuild docs

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