Context Navigation

source: trunk/GSASIImath.py @ 1513

Last change on this file since 1513 was 1513, checked in by toby, 9 years ago
update and rebuild docs
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 107.1 KB

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

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

Download in other formats: