Context Navigation

source: trunk/GSASIImath.py @ 4213

Last change on this file since 4213 was 4213, checked in by toby, 3 years ago
multiple small changes to allow docs build under 3.x
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 191.2 KB

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