Context Navigation

source: trunk/GSASIImath.py @ 4764

Last change on this file since 4764 was 4764, checked in by toby, 2 years ago
HessianLSQ: force vars to float; small wx 4.1 fixes; scriptable: more docs, TOF peak widths (from Adam)
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 208.6 KB

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