Context Navigation

source: trunk/GSASIImath.py @ 4025

Last change on this file since 4025 was 4025, checked in by toby, 4 years ago
more print filtering
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 191.2 KB

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