Context Navigation

source: trunk/GSASIImath_new.py @ 4699

Last change on this file since 4699 was 4699, checked in by toby, 2 years ago
improve refinement warnings
Property svn:eol-style set to `native`
File size: 206.9 KB

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