Context Navigation

source: trunk/GSASIImath.py @ 1830

Last change on this file since 1830 was 1830, checked in by vondreele, 9 years ago
change ftol on texture refinement - stop on real convergence improve spacing for texture display remove up/down buttons from Ddata display; do set focus back to this window when up/down arrows used to scan thru data
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 130.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: