Context Navigation

source: trunk/GSASIImath.py @ 4185

Last change on this file since 4185 was 4185, checked in by vondreele, 4 years ago
mag math again
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 191.1 KB

Line
1	# -- coding: utf-8 --
2	#GSASIImath - major mathematics routines
3	########### SVN repository information ###################
4	# $Date: 2019-10-21 20:17:29 +0000 (Mon, 21 Oct 2019) $
5	# $Author: vondreele $
6	# $Revision: 4185 $
7	# $URL: trunk/GSASIImath.py $
8	# $Id: GSASIImath.py 4185 2019-10-21 20:17:29Z vondreele $
9	########### SVN repository information ###################
10	'''
11	GSASIImath: computation module
12	================================
13
14	Routines for least-squares minimization and other stuff
15
16	'''
17	from __future__ import division, print_function
18	import random as rn
19	import numpy as np
20	import numpy.linalg as nl
21	import numpy.ma as ma
22	import time
23	import math
24	import copy
25	import GSASIIpath
26	GSASIIpath.SetVersionNumber("$Revision: 4185 $")
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(glTau,XYZ,modQ,MSSdata,SGData,SSGData):
1363	'''
1364	this needs to make magnetic moment modulations & magnitudes as
1365	fxn of gTau points
1366	'''
1367	Am = np.array(MSSdata[3:]).T #atoms x waves x cos pos mods
1368	Bm = np.array(MSSdata[:3]).T #...sin pos mods
1369	nWaves = Am.shape[1]
1370	SGMT = np.array([ops[0] for ops in SGData['SGOps']]) #not .T!!
1371	Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
1372	SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1373	if SGData['SGInv']:
1374	SGMT = np.vstack((SGMT,-SGMT))
1375	Sinv =np.vstack((Sinv,-Sinv))
1376	SGT = np.vstack((SGT,-SGT))
1377	SGMT = np.vstack([SGMT for cen in SGData['SGCen']])
1378	Sinv = np.vstack([Sinv for cen in SGData['SGCen']])
1379	SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
1380	if SGData['SGGray']:
1381	SGMT = np.vstack((SGMT,SGMT))
1382	Sinv =np.vstack((Sinv,Sinv))
1383	SGT = np.vstack((SGT,SGT+.5))%1.
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]*(glTau[:,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	MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)/2. #cos term
1394	MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)/2. #sin term
1395	MmodA = np.sum(SGMT[nxs,:,nxs,:,:]MmodA[:,:,:,nxs,:],axis=-1)SGData['MagMom'][nxs,:,nxs,nxs]
1396	MmodB = np.sum(SGMT[nxs,:,nxs,:,:]MmodB[:,:,:,nxs,:],axis=-1)SGData['MagMom'][nxs,:,nxs,nxs]
1397	return MmodA,MmodB #Ntau,Nops,Natm,,Mxyz; cos & sin parts; sum matches drawn atom moments
1398
1399	def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1400	'''
1401	H: array nRefBlk x ops X hklt
1402	HP: array nRefBlk x ops X hklt proj to hkl
1403	nWaves: list number of waves for frac, pos, uij & mag
1404	Fmod: array 2 x atoms x waves (sin,cos terms)
1405	Xmod: array atoms X 3 X ngl
1406	Umod: array atoms x 3x3 x ngl
1407	glTau,glWt: arrays Gauss-Lorentzian pos & wts
1408	'''
1409
1410	if nWaves[2]: #uij (adp) waves
1411	if len(HP.shape) > 2:
1412	HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1413	else:
1414	HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1415	else:
1416	HbH = 1.0
1417	HdotX = np.inner(HP,Xmod) #refBlk x ops x atoms X ngl
1418	if len(H.shape) > 2:
1419	D = H[:,:,3:]glTau[nxs,nxs,:] #me*tau; refBlk x ops X ngl
1420	HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1421	else:
1422	D = H[:,3:]glTau[nxs,:] #me*tau; refBlk x ops X ngl
1423	HdotXD = twopi*(HdotX+D[:,nxs,:])
1424	cosHA = np.sum(FmodHbHnp.cos(HdotXD)*glWt,axis=-1) #real part; refBlk X ops x atoms; sum for G-L integration
1425	sinHA = np.sum(FmodHbHnp.sin(HdotXD)*glWt,axis=-1) #imag part; ditto
1426	return np.array([cosHA,sinHA]) # 2 x refBlk x SGops x atoms
1427
1428	#def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
1429	# '''
1430	# H: array nRefBlk x ops X hklt
1431	# HP: array nRefBlk x ops X hklt proj to hkl
1432	# nWaves: list number of waves for frac, pos, uij & mag
1433	# Fmod: array 2 x atoms x waves (sin,cos terms)
1434	# Xmod: array atoms X 3 X ngl
1435	# Umod: array atoms x 3x3 x ngl
1436	# Mmod: array atoms x 3x3 x ngl
1437	# glTau,glWt: arrays Gauss-Lorentzian pos & wts
1438	# '''
1439	#
1440	# if nWaves[2]: #uij (adp) waves
1441	# if len(HP.shape) > 2:
1442	# HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1443	# else:
1444	# HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1445	# else:
1446	# HbH = 1.0
1447	# HdotX = np.inner(HP,Xmod) #refBlk x ops x atoms X ngl
1448	# if len(H.shape) > 2:
1449	# D = H[:,:,3:]glTau[nxs,nxs,:] #me*tau; refBlk x ops X ngl
1450	# HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1451	# else:
1452	# D = H[:,3:]glTau[nxs,:] #me*tau; refBlk x ops X ngl
1453	# HdotXD = twopi*(HdotX+D[:,nxs,:])
1454	# M = np.swapaxes(Mmod,1,2)
1455	# 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
1456	# sinHA = np.sum(M[nxs,nxs,:,:,:](FmodHbHnp.sin(HdotXD)[:,:,:,nxs,:]glWt),axis=-1) #imag part; ditto
1457	# return np.array([cosHA,sinHA]) # 2 x refBlk x SGops x atoms
1458	#
1459	def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1460	'''
1461	H: array nRefBlk x tw x ops X hklt
1462	HP: array nRefBlk x tw x ops X hklt proj to hkl
1463	Fmod: array 2 x atoms x waves (sin,cos terms)
1464	Xmod: array atoms X ngl X 3
1465	Umod: array atoms x ngl x 3x3
1466	glTau,glWt: arrays Gauss-Lorentzian pos & wts
1467	'''
1468
1469	if nWaves[2]:
1470	if len(HP.shape) > 3: #Blocks of reflections
1471	HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1472	else: #single reflections
1473	HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1474	else:
1475	HbH = 1.0
1476	HdotX = np.inner(HP,Xmod) #refBlk x tw x ops x atoms X ngl
1477	if len(H.shape) > 3:
1478	D = glTauH[:,:,:,3:,nxs] #me*tau; refBlk x tw x ops X ngl
1479	HdotXD = twopi*(HdotX+D[:,:,:,nxs,:])
1480	else:
1481	D = HglTau[nxs,:] #me*tau; refBlk x ops X ngl
1482	HdotXD = twopi*(HdotX+D[:,nxs,:])
1483	cosHA = np.sum(FmodHbHnp.cos(HdotXD)*glWt,axis=-1) #real part; refBlk X ops x atoms; sum for G-L integration
1484	sinHA = np.sum(FmodHbHnp.sin(HdotXD)*glWt,axis=-1) #imag part; ditto
1485	return np.array([cosHA,sinHA]) # 2 x refBlk x SGops x atoms
1486
1487	def makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast):
1488	'''
1489	Only for Fourier waves for fraction, position & adp (probably not used for magnetism)
1490	FSSdata: array 2 x atoms x waves (sin,cos terms)
1491	XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1492	USSdata: array 2x6 x atoms X waves (sin,cos terms)
1493	Mast: array orthogonalization matrix for Uij
1494	'''
1495	glTau,glWt = pwd.pygauleg(0.,1.,ngl) #get Gauss-Legendre intervals & weights
1496	waveShapes = [FSSdata.T.shape,XSSdata.T.shape,USSdata.T.shape]
1497	Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms
1498	Bf = np.array(FSSdata[1]).T #cos frac mods...
1499	Ax = np.array(XSSdata[:3]).T #...cos pos mods x waves x atoms
1500	Bx = np.array(XSSdata[3:]).T #...cos pos mods
1501	Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods
1502	Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods
1503	nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]]
1504	StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl)) #atoms x waves x 3 x ngl
1505	CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1506	ZtauXt = np.zeros((Ax.shape[0],2,3,ngl)) #atoms x Tminmax x 3 x ngl
1507	ZtauXx = np.zeros((Ax.shape[0],3,ngl)) #atoms x XYZmax x ngl
1508	for iatm in range(Ax.shape[0]):
1509	nx = 0
1510	if 'ZigZag' in waveTypes[iatm]:
1511	nx = 1
1512	elif 'Block' in waveTypes[iatm]:
1513	nx = 1
1514	tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau #Xwaves x ngl
1515	if nx:
1516	StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]np.sin(twopitauX)[nxs,:,nxs,:] #atoms X waves X 3(xyz) X ngl
1517	CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]np.cos(twopitauX)[nxs,:,nxs,:] #ditto
1518	else:
1519	StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]np.sin(twopitauX)[nxs,:,nxs,:] #atoms X waves X 3(xyz) X ngl
1520	CtauX[iatm] = np.ones_like(Bx)[iatm,:,:,nxs]np.cos(twopitauX)[nxs,:,nxs,:] #ditto
1521	if nWaves[0]:
1522	tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau #Fwaves x ngl
1523	StauF = np.ones_like(Af)[:,:,nxs]np.sin(twopitauF)[nxs,:,:] #also dFmod/dAf
1524	CtauF = np.ones_like(Bf)[:,:,nxs]np.cos(twopitauF)[nxs,:,:] #also dFmod/dBf
1525	else:
1526	StauF = 1.0
1527	CtauF = 1.0
1528	if nWaves[2]:
1529	tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau #Uwaves x ngl
1530	StauU = np.ones_like(Au)[:,:,:,:,nxs]np.sin(twopitauU)[nxs,:,nxs,nxs,:] #also dUmodA/dAu
1531	CtauU = np.ones_like(Bu)[:,:,:,:,nxs]np.cos(twopitauU)[nxs,:,nxs,nxs,:] #also dUmodB/dBu
1532	UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x ngl
1533	UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
1534	#derivs need to be ops x atoms x waves x 6uij; ops x atoms x waves x ngl x 6uij before sum
1535	StauU = np.rollaxis(np.rollaxis(np.swapaxes(StauU,2,4),-1),-1)
1536	CtauU = np.rollaxis(np.rollaxis(np.swapaxes(CtauU,2,4),-1),-1)
1537	else:
1538	StauU = 1.0
1539	CtauU = 1.0
1540	UmodA = 0.
1541	UmodB = 0.
1542	return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
1543
1544	def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
1545	'''
1546	Compute Fourier modulation derivatives
1547	H: array ops X hklt proj to hkl
1548	HP: array ops X hklt proj to hkl
1549	Hij: array 2pi^2[a^2h^2 b^2k^2 c^2l^2 abhk achl bc*kl] of projected hklm to hkl space
1550	'''
1551
1552	Mf = [H.shape[0],]+list(waveShapes[0]) #=[ops,atoms,waves,2] (sin+cos frac mods)
1553	dGdMfC = np.zeros(Mf)
1554	dGdMfS = np.zeros(Mf)
1555	Mx = [H.shape[0],]+list(waveShapes[1]) #=[ops,atoms,waves,6] (sin+cos pos mods)
1556	dGdMxC = np.zeros(Mx)
1557	dGdMxS = np.zeros(Mx)
1558	Mu = [H.shape[0],]+list(waveShapes[2]) #=[ops,atoms,waves,12] (sin+cos Uij mods)
1559	dGdMuC = np.zeros(Mu)
1560	dGdMuS = np.zeros(Mu)
1561
1562	D = twopiH[:,3][:,nxs]glTau[nxs,:] #metau; ops X ngl
1563	HdotX = twopi*np.inner(HP,Xmod) #ops x atoms X ngl
1564	HdotXD = HdotX+D[:,nxs,:]
1565	if nWaves[2]:
1566	Umod = np.swapaxes((UmodAB),2,4) #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
1567	HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1) #ops x atoms x waves x ngl
1568	HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1) #ops x atoms x waves x ngl
1569	HbH = np.exp(-np.sum(HuH,axis=-2)) # ops x atoms x ngl; sum waves - OK vs Modulation version
1570	# part1 = -np.exp(-HuH)*Fmod[nxs,:,nxs,:] #ops x atoms x waves x ngl
1571	part1 = -np.exp(-HuH)*Fmod #ops x atoms x waves x ngl
1572	dUdAu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,:,:,:,:] #ops x atoms x waves x ngl x 6sinUij
1573	dUdBu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,:,:,:,:] #ops x atoms x waves x ngl x 6cosUij
1574	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
1575	dGdMuCb = np.sum(part1[:,:,:,:,nxs]dUdBunp.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1576	dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1) #ops x atoms x waves x 12uij
1577	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
1578	dGdMuSb = np.sum(part1[:,:,:,:,nxs]dUdBunp.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1579	dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1) #ops x atoms x waves x 12uij
1580	else:
1581	HbH = np.ones_like(HdotX)
1582	dHdXA = twopiHP[:,nxs,nxs,nxs,:]np.swapaxes(SCtauX[0],-1,-2)[nxs,:,:,:,:] #ops x atoms x sine waves x ngl x xyz
1583	dHdXB = twopiHP[:,nxs,nxs,nxs,:]np.swapaxes(SCtauX[1],-1,-2)[nxs,:,:,:,:] #ditto - cos waves
1584	# ops x atoms x waves x 2xyz - real part - good
1585	# dGdMxCa = -np.sum((Fmod[nxs,:,:]HbH)[:,:,nxs,:,nxs](dHdXAnp.sin(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1586	# dGdMxCb = -np.sum((Fmod[nxs,:,:]HbH)[:,:,nxs,:,nxs](dHdXBnp.sin(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1587	dGdMxCa = -np.sum((FmodHbH)[:,:,nxs,:,nxs](dHdXAnp.sin(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1588	dGdMxCb = -np.sum((FmodHbH)[:,:,nxs,:,nxs](dHdXBnp.sin(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1589	dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
1590	# ops x atoms x waves x 2xyz - imag part - good
1591	# dGdMxSa = np.sum((Fmod[nxs,:,:]HbH)[:,:,nxs,:,nxs](dHdXAnp.cos(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1592	# dGdMxSb = np.sum((Fmod[nxs,:,:]HbH)[:,:,nxs,:,nxs](dHdXBnp.cos(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1593	dGdMxSa = np.sum((FmodHbH)[:,:,nxs,:,nxs](dHdXAnp.cos(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1594	dGdMxSb = np.sum((FmodHbH)[:,:,nxs,:,nxs](dHdXBnp.cos(HdotXD)[:,:,nxs,:,nxs])glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1595	dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
1596	return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS]
1597
1598	def posFourier(tau,psin,pcos):
1599	A = np.array([ps[:,nxs]np.sin(2np.pi(i+1)tau) for i,ps in enumerate(psin)])
1600	B = np.array([pc[:,nxs]np.cos(2np.pi(i+1)tau) for i,pc in enumerate(pcos)])
1601	return np.sum(A,axis=0)+np.sum(B,axis=0)
1602
1603	def posZigZag(T,Tmm,Xmax):
1604	DT = Tmm[1]-Tmm[0]
1605	Su = 2.*Xmax/DT
1606	Sd = 2.*Xmax/(1.-DT)
1607	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])
1608	return A
1609
1610	#def posZigZagDerv(T,Tmm,Xmax):
1611	# DT = Tmm[1]-Tmm[0]
1612	# Su = 2.*Xmax/DT
1613	# Sd = 2.*Xmax/(1.-DT)
1614	# dAdT = np.zeros((2,3,len(T)))
1615	# 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
1616	# 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
1617	# 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])
1618	# return dAdT,dAdX
1619
1620	def posBlock(T,Tmm,Xmax):
1621	A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax,Xmax) for t in T])
1622	return A
1623
1624	#def posBlockDerv(T,Tmm,Xmax):
1625	# dAdT = np.zeros((2,3,len(T)))
1626	# ind = np.searchsorted(T,Tmm)
1627	# dAdT[0,:,ind[0]] = -Xmax/len(T)
1628	# dAdT[1,:,ind[1]] = Xmax/len(T)
1629	# dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t <= Tmm[1],-1.,1.) for t in T]) #OK
1630	# return dAdT,dAdX
1631
1632	def fracCrenel(tau,Toff,Twid):
1633	Tau = (tau-Toff)%1.
1634	A = np.where(Tau<Twid,1.,0.)
1635	return A
1636
1637	def fracFourier(tau,fsin,fcos):
1638	if len(fsin) == 1:
1639	A = np.array([fsin[0]np.sin(2.np.pi*tau)])
1640	B = np.array([fcos[0]np.cos(2.np.pi*tau)])
1641	else:
1642	A = np.array([fs[:,nxs]np.sin(2.np.pi(i+1)tau) for i,fs in enumerate(fsin)])
1643	B = np.array([fc[:,nxs]np.cos(2.np.pi(i+1)tau) for i,fc in enumerate(fcos)])
1644	return np.sum(A,axis=0)+np.sum(B,axis=0)
1645
1646	def ApplyModulation(data,tau):
1647	'''Applies modulation to drawing atom positions & Uijs for given tau
1648	'''
1649	generalData = data['General']
1650	cell = generalData['Cell'][1:7]
1651	G,g = G2lat.cell2Gmat(cell)
1652	SGData = generalData['SGData']
1653	SSGData = generalData['SSGData']
1654	cx,ct,cs,cia = generalData['AtomPtrs']
1655	drawingData = data['Drawing']
1656	modul = generalData['SuperVec'][0]
1657	dcx,dct,dcs,dci = drawingData['atomPtrs']
1658	atoms = data['Atoms']
1659	drawAtoms = drawingData['Atoms']
1660	Fade = np.ones(len(drawAtoms))
1661	for atom in atoms:
1662	atxyz = np.array(atom[cx:cx+3])
1663	atuij = np.array(atom[cia+2:cia+8])
1664	Sfrac = atom[-1]['SS1']['Sfrac']
1665	Spos = atom[-1]['SS1']['Spos']
1666	Sadp = atom[-1]['SS1']['Sadp']
1667	if generalData['Type'] == 'magnetic':
1668	Smag = atom[-1]['SS1']['Smag']
1669	atmom = np.array(atom[cx+4:cx+7])
1670	indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
1671	for ind in indx:
1672	drawatom = drawAtoms[ind]
1673	opr = drawatom[dcs-1]
1674	sop,ssop,icent,cent,unit = G2spc.OpsfromStringOps(opr,SGData,SSGData)
1675	drxyz = (np.inner(sop[0],atxyz)+sop[1]+cent)*icent+np.array(unit)
1676	tauT = G2spc.getTauT(tau,sop,ssop,drxyz,modul)[-1]
1677	tauT *= icent #invert wave on -1
1678	# print(tau,tauT,opr,G2spc.MT2text(sop).replace(' ',''),G2spc.SSMT2text(ssop).replace(' ',''))
1679	wave = np.zeros(3)
1680	uwave = np.zeros(6)
1681	mom = np.zeros(3)
1682	if len(Sfrac):
1683	scof = []
1684	ccof = []
1685	waveType = Sfrac[0]
1686	for i,sfrac in enumerate(Sfrac[1:]):
1687	if not i and 'Crenel' in waveType:
1688	Fade[ind] += fracCrenel(tauT,sfrac[0][0],sfrac[0][1])
1689	else:
1690	scof.append(sfrac[0][0])
1691	ccof.append(sfrac[0][1])
1692	if len(scof):
1693	Fade[ind] += np.sum(fracFourier(tauT,scof,ccof))
1694	if len(Spos):
1695	scof = []
1696	ccof = []
1697	waveType = Spos[0]
1698	for i,spos in enumerate(Spos[1:]):
1699	if waveType in ['ZigZag','Block'] and not i:
1700	Tminmax = spos[0][:2]
1701	XYZmax = np.array(spos[0][2:5])
1702	if waveType == 'Block':
1703	wave = np.array(posBlock([tauT,],Tminmax,XYZmax))[0]
1704	elif waveType == 'ZigZag':
1705	wave = np.array(posZigZag([tauT,],Tminmax,XYZmax))[0]
1706	else:
1707	scof.append(spos[0][:3])
1708	ccof.append(spos[0][3:])
1709	if len(scof):
1710	wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1711	if generalData['Type'] == 'magnetic' and len(Smag):
1712	scof = []
1713	ccof = []
1714	waveType = Smag[0]
1715	for i,spos in enumerate(Smag[1:]):
1716	scof.append(spos[0][:3])
1717	ccof.append(spos[0][3:])
1718	if len(scof):
1719	mom += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1720	if len(Sadp):
1721	scof = []
1722	ccof = []
1723	waveType = Sadp[0]
1724	for i,sadp in enumerate(Sadp[1:]):
1725	scof.append(sadp[0][:6])
1726	ccof.append(sadp[0][6:])
1727	ures = posFourier(tauT,np.array(scof),np.array(ccof))
1728	if np.any(ures):
1729	uwave += np.sum(ures,axis=1)
1730	if atom[cia] == 'A':
1731	X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
1732	drawatom[dcx:dcx+3] = X
1733	drawatom[dci-6:dci] = U
1734	else:
1735	X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
1736	drawatom[dcx:dcx+3] = X
1737	if generalData['Type'] == 'magnetic':
1738	M = G2spc.ApplyStringOpsMom(opr,SGData,atmom+mom)
1739	drawatom[dcx+3:dcx+6] = M
1740	return drawAtoms,Fade
1741
1742	# gauleg.py Gauss Legendre numerical quadrature, x and w computation
1743	# integrate from a to b using n evaluations of the function f(x)
1744	# usage: from gauleg import gaulegf
1745	# x,w = gaulegf( a, b, n)
1746	# area = 0.0
1747	# for i in range(1,n+1): # yes, 1..n
1748	# area += w[i]*f(x[i])
1749
1750	def gaulegf(a, b, n):
1751	x = range(n+1) # x[0] unused
1752	w = range(n+1) # w[0] unused
1753	eps = 3.0E-14
1754	m = (n+1)/2
1755	xm = 0.5*(b+a)
1756	xl = 0.5*(b-a)
1757	for i in range(1,m+1):
1758	z = math.cos(3.141592654*(i-0.25)/(n+0.5))
1759	while True:
1760	p1 = 1.0
1761	p2 = 0.0
1762	for j in range(1,n+1):
1763	p3 = p2
1764	p2 = p1
1765	p1 = ((2.0j-1.0)zp2-(j-1.0)p3)/j
1766
1767	pp = n(zp1-p2)/(z*z-1.0)
1768	z1 = z
1769	z = z1 - p1/pp
1770	if abs(z-z1) <= eps:
1771	break
1772
1773	x[i] = xm - xl*z
1774	x[n+1-i] = xm + xl*z
1775	w[i] = 2.0xl/((1.0-zz)pppp)
1776	w[n+1-i] = w[i]
1777	return np.array(x), np.array(w)
1778	# end gaulegf
1779
1780
1781	def BessJn(nmax,x):
1782	''' compute Bessel function J(n,x) from scipy routine & recurrance relation
1783	returns sequence of J(n,x) for n in range [-nmax...0...nmax]
1784
1785	:param integer nmax: maximul order for Jn(x)
1786	:param float x: argument for Jn(x)
1787
1788	:returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
1789
1790	'''
1791	import scipy.special as sp
1792	bessJn = np.zeros(2*nmax+1)
1793	bessJn[nmax] = sp.j0(x)
1794	bessJn[nmax+1] = sp.j1(x)
1795	bessJn[nmax-1] = -bessJn[nmax+1]
1796	for i in range(2,nmax+1):
1797	bessJn[i+nmax] = 2(i-1)bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
1798	bessJn[nmax-i] = bessJn[i+nmax](-1)*i
1799	return bessJn
1800
1801	def BessIn(nmax,x):
1802	''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
1803	returns sequence of I(n,x) for n in range [-nmax...0...nmax]
1804
1805	:param integer nmax: maximul order for In(x)
1806	:param float x: argument for In(x)
1807
1808	:returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
1809
1810	'''
1811	import scipy.special as sp
1812	bessIn = np.zeros(2*nmax+1)
1813	bessIn[nmax] = sp.i0(x)
1814	bessIn[nmax+1] = sp.i1(x)
1815	bessIn[nmax-1] = bessIn[nmax+1]
1816	for i in range(2,nmax+1):
1817	bessIn[i+nmax] = bessIn[nmax+i-2]-2(i-1)bessIn[nmax+i-1]/x
1818	bessIn[nmax-i] = bessIn[i+nmax]
1819	return bessIn
1820
1821
1822	################################################################################
1823	##### distance, angle, planes, torsion stuff
1824	################################################################################
1825
1826	def CalcDist(distance_dict, distance_atoms, parmDict):
1827	if not len(parmDict):
1828	return 0.
1829	pId = distance_dict['pId']
1830	A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1831	Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1832	Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1833	Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1834	inv = 1
1835	symNo = distance_dict['symNo']
1836	if symNo < 0:
1837	inv = -1
1838	symNo *= -1
1839	cen = symNo//100
1840	op = symNo%100-1
1841	M,T = distance_dict['SGData']['SGOps'][op]
1842	D = T*inv+distance_dict['SGData']['SGCen'][cen]
1843	D += distance_dict['cellNo']
1844	Txyz = np.inner(M*inv,Txyz)+D
1845	dist = np.sqrt(np.sum(np.inner(Amat,(Txyz-Oxyz))**2))
1846	# GSASIIpath.IPyBreak()
1847	return dist
1848
1849	def CalcDistDeriv(distance_dict, distance_atoms, parmDict):
1850	if not len(parmDict):
1851	return None
1852	pId = distance_dict['pId']
1853	A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1854	Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1855	Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1856	Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1857	symNo = distance_dict['symNo']
1858	Tunit = distance_dict['cellNo']
1859	SGData = distance_dict['SGData']
1860	deriv = getDistDerv(Oxyz,Txyz,Amat,Tunit,symNo,SGData)
1861	return deriv
1862
1863	def CalcAngle(angle_dict, angle_atoms, parmDict):
1864	if not len(parmDict):
1865	return 0.
1866	pId = angle_dict['pId']
1867	A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1868	Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1869	Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1870	Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1871	Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1872	ABxyz = [Axyz,Bxyz]
1873	symNo = angle_dict['symNo']
1874	vec = np.zeros((2,3))
1875	for i in range(2):
1876	inv = 1
1877	if symNo[i] < 0:
1878	inv = -1
1879	cen = inv*symNo[i]//100
1880	op = inv*symNo[i]%100-1
1881	M,T = angle_dict['SGData']['SGOps'][op]
1882	D = T*inv+angle_dict['SGData']['SGCen'][cen]
1883	D += angle_dict['cellNo'][i]
1884	ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
1885	vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
1886	dist = np.sqrt(np.sum(vec[i]**2))
1887	if not dist:
1888	return 0.
1889	vec[i] /= dist
1890	angle = acosd(np.sum(vec[0]*vec[1]))
1891	# GSASIIpath.IPyBreak()
1892	return angle
1893
1894	def CalcAngleDeriv(angle_dict, angle_atoms, parmDict):
1895	if not len(parmDict):
1896	return None
1897	pId = angle_dict['pId']
1898	A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1899	Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1900	Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1901	Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1902	Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1903	symNo = angle_dict['symNo']
1904	Tunit = angle_dict['cellNo']
1905	SGData = angle_dict['SGData']
1906	deriv = getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData)
1907	return deriv
1908
1909	def getSyXYZ(XYZ,ops,SGData):
1910	'''default doc
1911
1912	:param type name: description
1913
1914	:returns: type name: description
1915
1916	'''
1917	XYZout = np.zeros_like(XYZ)
1918	for i,[xyz,op] in enumerate(zip(XYZ,ops)):
1919	if op == '1':
1920	XYZout[i] = xyz
1921	else:
1922	oprs = op.split('+')
1923	unit = [0,0,0]
1924	if len(oprs)>1:
1925	unit = np.array(list(eval(oprs[1])))
1926	syop =int(oprs[0])
1927	inv = syop//abs(syop)
1928	syop *= inv
1929	cent = syop//100
1930	syop %= 100
1931	syop -= 1
1932	M,T = SGData['SGOps'][syop]
1933	XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
1934	return XYZout
1935
1936	def getRestDist(XYZ,Amat):
1937	'''default doc string
1938
1939	:param type name: description
1940
1941	:returns: type name: description
1942
1943	'''
1944	return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1945
1946	def getRestDeriv(Func,XYZ,Amat,ops,SGData):
1947	'''default doc string
1948
1949	:param type name: description
1950
1951	:returns: type name: description
1952
1953	'''
1954	deriv = np.zeros((len(XYZ),3))
1955	dx = 0.00001
1956	for j,xyz in enumerate(XYZ):
1957	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1958	XYZ[j] -= x
1959	d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1960	XYZ[j] += 2*x
1961	d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1962	XYZ[j] -= x
1963	deriv[j][i] = (d1-d2)/(2*dx)
1964	return deriv.flatten()
1965
1966	def getRestAngle(XYZ,Amat):
1967	'''default doc string
1968
1969	:param type name: description
1970
1971	:returns: type name: description
1972
1973	'''
1974
1975	def calcVec(Ox,Tx,Amat):
1976	return np.inner(Amat,(Tx-Ox))
1977
1978	VecA = calcVec(XYZ[1],XYZ[0],Amat)
1979	VecA /= np.sqrt(np.sum(VecA**2))
1980	VecB = calcVec(XYZ[1],XYZ[2],Amat)
1981	VecB /= np.sqrt(np.sum(VecB**2))
1982	edge = VecB-VecA
1983	edge = np.sum(edge**2)
1984	angle = (2.-edge)/2.
1985	angle = max(angle,-1.)
1986	return acosd(angle)
1987
1988	def getRestPlane(XYZ,Amat):
1989	'''default doc string
1990
1991	:param type name: description
1992
1993	:returns: type name: description
1994
1995	'''
1996	sumXYZ = np.zeros(3)
1997	for xyz in XYZ:
1998	sumXYZ += xyz
1999	sumXYZ /= len(XYZ)
2000	XYZ = np.array(XYZ)-sumXYZ
2001	XYZ = np.inner(Amat,XYZ).T
2002	Zmat = np.zeros((3,3))
2003	for i,xyz in enumerate(XYZ):
2004	Zmat += np.outer(xyz.T,xyz)
2005	Evec,Emat = nl.eig(Zmat)
2006	Evec = np.sqrt(Evec)/(len(XYZ)-3)
2007	Order = np.argsort(Evec)
2008	return Evec[Order[0]]
2009
2010	def getRestChiral(XYZ,Amat):
2011	'''default doc string
2012
2013	:param type name: description
2014
2015	:returns: type name: description
2016
2017	'''
2018	VecA = np.empty((3,3))
2019	VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
2020	VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
2021	VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
2022	return nl.det(VecA)
2023
2024	def getRestTorsion(XYZ,Amat):
2025	'''default doc string
2026
2027	:param type name: description
2028
2029	:returns: type name: description
2030
2031	'''
2032	VecA = np.empty((3,3))
2033	VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
2034	VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
2035	VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
2036	D = nl.det(VecA)
2037	Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
2038	P12 = np.sum(VecA[0]VecA[1])/(Mag[0]Mag[1])
2039	P13 = np.sum(VecA[0]VecA[2])/(Mag[0]Mag[2])
2040	P23 = np.sum(VecA[1]VecA[2])/(Mag[1]Mag[2])
2041	Ang = 1.0
2042	if abs(P12) < 1.0 and abs(P23) < 1.0:
2043	Ang = (P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23**2))
2044	TOR = (acosd(Ang)*D/abs(D)+720.)%360.
2045	return TOR
2046
2047	def calcTorsionEnergy(TOR,Coeff=[]):
2048	'''default doc string
2049
2050	:param type name: description
2051
2052	:returns: type name: description
2053
2054	'''
2055	sum = 0.
2056	if len(Coeff):
2057	cof = np.reshape(Coeff,(3,3)).T
2058	delt = TOR-cof[1]
2059	delt = np.where(delt<-180.,delt+360.,delt)
2060	delt = np.where(delt>180.,delt-360.,delt)
2061	term = -cof[2]delt*2
2062	val = cof[0]*np.exp(term/1000.0)
2063	pMax = cof[0][np.argmin(val)]
2064	Eval = np.sum(val)
2065	sum = Eval-pMax
2066	return sum,Eval
2067
2068	def getTorsionDeriv(XYZ,Amat,Coeff):
2069	'''default doc string
2070
2071	:param type name: description
2072
2073	:returns: type name: description
2074
2075	'''
2076	deriv = np.zeros((len(XYZ),3))
2077	dx = 0.00001
2078	for j,xyz in enumerate(XYZ):
2079	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2080	XYZ[j] -= x
2081	tor = getRestTorsion(XYZ,Amat)
2082	p1,d1 = calcTorsionEnergy(tor,Coeff)
2083	XYZ[j] += 2*x
2084	tor = getRestTorsion(XYZ,Amat)
2085	p2,d2 = calcTorsionEnergy(tor,Coeff)
2086	XYZ[j] -= x
2087	deriv[j][i] = (p2-p1)/(2*dx)
2088	return deriv.flatten()
2089
2090	def getRestRama(XYZ,Amat):
2091	'''Computes a pair of torsion angles in a 5 atom string
2092
2093	:param nparray XYZ: crystallographic coordinates of 5 atoms
2094	:param nparray Amat: crystal to cartesian transformation matrix
2095
2096	:returns: list (phi,psi) two torsion angles in degrees
2097
2098	'''
2099	phi = getRestTorsion(XYZ[:5],Amat)
2100	psi = getRestTorsion(XYZ[1:],Amat)
2101	return phi,psi
2102
2103	def calcRamaEnergy(phi,psi,Coeff=[]):
2104	'''Computes pseudo potential energy from a pair of torsion angles and a
2105	numerical description of the potential energy surface. Used to create
2106	penalty function in LS refinement:
2107	:math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
2108
2109	where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4](\\psi-C[2])^2 - 2(\\phi-C[1])*(\\psi-C[2])`
2110
2111	:param float phi: first torsion angle (:math:`\\phi`)
2112	:param float psi: second torsion angle (:math:`\\psi`)
2113	:param list Coeff: pseudo potential coefficients
2114
2115	:returns: list (sum,Eval): pseudo-potential difference from minimum & value;
2116	sum is used for penalty function.
2117
2118	'''
2119	sum = 0.
2120	Eval = 0.
2121	if len(Coeff):
2122	cof = Coeff.T
2123	dPhi = phi-cof[1]
2124	dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
2125	dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
2126	dPsi = psi-cof[2]
2127	dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
2128	dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
2129	val = -cof[3]dPhi2-cof[4]dPsi*2-2.0cof[5]dPhidPsi
2130	val = cof[0]*np.exp(val/1000.)
2131	pMax = cof[0][np.argmin(val)]
2132	Eval = np.sum(val)
2133	sum = Eval-pMax
2134	return sum,Eval
2135
2136	def getRamaDeriv(XYZ,Amat,Coeff):
2137	'''Computes numerical derivatives of torsion angle pair pseudo potential
2138	with respect of crystallographic atom coordinates of the 5 atom sequence
2139
2140	:param nparray XYZ: crystallographic coordinates of 5 atoms
2141	:param nparray Amat: crystal to cartesian transformation matrix
2142	:param list Coeff: pseudo potential coefficients
2143
2144	:returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
2145	crystallographic xyz coordinates.
2146
2147	'''
2148	deriv = np.zeros((len(XYZ),3))
2149	dx = 0.00001
2150	for j,xyz in enumerate(XYZ):
2151	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2152	XYZ[j] -= x
2153	phi,psi = getRestRama(XYZ,Amat)
2154	p1,d1 = calcRamaEnergy(phi,psi,Coeff)
2155	XYZ[j] += 2*x
2156	phi,psi = getRestRama(XYZ,Amat)
2157	p2,d2 = calcRamaEnergy(phi,psi,Coeff)
2158	XYZ[j] -= x
2159	deriv[j][i] = (p2-p1)/(2*dx)
2160	return deriv.flatten()
2161
2162	def getRestPolefig(ODFln,SamSym,Grid):
2163	'''default doc string
2164
2165	:param type name: description
2166
2167	:returns: type name: description
2168
2169	'''
2170	X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
2171	R,P = np.sqrt(X2+Y2).flatten(),atan2d(Y,X).flatten()
2172	R = np.where(R <= 1.,2.*atand(R),0.0)
2173	Z = np.zeros_like(R)
2174	Z = G2lat.polfcal(ODFln,SamSym,R,P)
2175	Z = np.reshape(Z,(Grid,Grid))
2176	return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
2177
2178	def getRestPolefigDerv(HKL,Grid,SHCoeff):
2179	'''default doc string
2180
2181	:param type name: description
2182
2183	:returns: type name: description
2184
2185	'''
2186	pass
2187
2188	def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
2189	'''default doc string
2190
2191	:param type name: description
2192
2193	:returns: type name: description
2194
2195	'''
2196	def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
2197	TxT = inv*(np.inner(M,Tx)+T)+C+U
2198	return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
2199
2200	inv = Top/abs(Top)
2201	cent = abs(Top)//100
2202	op = abs(Top)%100-1
2203	M,T = SGData['SGOps'][op]
2204	C = SGData['SGCen'][cent]
2205	dx = .00001
2206	deriv = np.zeros(6)
2207	for i in [0,1,2]:
2208	Oxyz[i] -= dx
2209	d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2210	Oxyz[i] += 2*dx
2211	deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2212	Oxyz[i] -= dx
2213	Txyz[i] -= dx
2214	d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2215	Txyz[i] += 2*dx
2216	deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2217	Txyz[i] -= dx
2218	return deriv
2219
2220	def getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData):
2221
2222	def calcAngle(Oxyz,ABxyz,Amat,Tunit,symNo,SGData):
2223	vec = np.zeros((2,3))
2224	for i in range(2):
2225	inv = 1
2226	if symNo[i] < 0:
2227	inv = -1
2228	cen = inv*symNo[i]//100
2229	op = inv*symNo[i]%100-1
2230	M,T = SGData['SGOps'][op]
2231	D = T*inv+SGData['SGCen'][cen]
2232	D += Tunit[i]
2233	ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
2234	vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
2235	dist = np.sqrt(np.sum(vec[i]**2))
2236	if not dist:
2237	return 0.
2238	vec[i] /= dist
2239	angle = acosd(np.sum(vec[0]*vec[1]))
2240	# GSASIIpath.IPyBreak()
2241	return angle
2242
2243	dx = .00001
2244	deriv = np.zeros(9)
2245	for i in [0,1,2]:
2246	Oxyz[i] -= dx
2247	a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2248	Oxyz[i] += 2*dx
2249	deriv[i] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2250	Oxyz[i] -= dx
2251	Axyz[i] -= dx
2252	a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2253	Axyz[i] += 2*dx
2254	deriv[i+3] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2255	Axyz[i] -= dx
2256	Bxyz[i] -= dx
2257	a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2258	Bxyz[i] += 2*dx
2259	deriv[i+6] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2260	Bxyz[i] -= dx
2261	return deriv
2262
2263	def getAngSig(VA,VB,Amat,SGData,covData={}):
2264	'''default doc string
2265
2266	:param type name: description
2267
2268	:returns: type name: description
2269
2270	'''
2271	def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
2272	TxT = inv*(np.inner(M,Tx)+T)+C+U
2273	return np.inner(Amat,(TxT-Ox))
2274
2275	def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
2276	VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
2277	VecA /= np.sqrt(np.sum(VecA**2))
2278	VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
2279	VecB /= np.sqrt(np.sum(VecB**2))
2280	edge = VecB-VecA
2281	edge = np.sum(edge**2)
2282	angle = (2.-edge)/2.
2283	angle = max(angle,-1.)
2284	return acosd(angle)
2285
2286	OxAN,OxA,TxAN,TxA,unitA,TopA = VA
2287	OxBN,OxB,TxBN,TxB,unitB,TopB = VB
2288	invA = invB = 1
2289	invA = TopA//abs(TopA)
2290	invB = TopB//abs(TopB)
2291	centA = abs(TopA)//100
2292	centB = abs(TopB)//100
2293	opA = abs(TopA)%100-1
2294	opB = abs(TopB)%100-1
2295	MA,TA = SGData['SGOps'][opA]
2296	MB,TB = SGData['SGOps'][opB]
2297	CA = SGData['SGCen'][centA]
2298	CB = SGData['SGCen'][centB]
2299	if 'covMatrix' in covData:
2300	covMatrix = covData['covMatrix']
2301	varyList = covData['varyList']
2302	AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
2303	dx = .00001
2304	dadx = np.zeros(9)
2305	Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2306	for i in [0,1,2]:
2307	OxA[i] -= dx
2308	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2309	OxA[i] += 2*dx
2310	dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2311	OxA[i] -= dx
2312
2313	TxA[i] -= dx
2314	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2315	TxA[i] += 2*dx
2316	dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2317	TxA[i] -= dx
2318
2319	TxB[i] -= dx
2320	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2321	TxB[i] += 2*dx
2322	dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2323	TxB[i] -= dx
2324
2325	sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2326	if sigAng < 0.01:
2327	sigAng = 0.0
2328	return Ang,sigAng
2329	else:
2330	return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
2331
2332	def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
2333	'''default doc string
2334
2335	:param type name: description
2336
2337	:returns: type name: description
2338
2339	'''
2340	def calcDist(Atoms,SyOps,Amat):
2341	XYZ = []
2342	for i,atom in enumerate(Atoms):
2343	Inv,M,T,C,U = SyOps[i]
2344	XYZ.append(np.array(atom[1:4]))
2345	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2346	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2347	V1 = XYZ[1]-XYZ[0]
2348	return np.sqrt(np.sum(V1**2))
2349
2350	SyOps = []
2351	names = []
2352	for i,atom in enumerate(Oatoms):
2353	names += atom[-1]
2354	Op,unit = Atoms[i][-1]
2355	inv = Op//abs(Op)
2356	m,t = SGData['SGOps'][abs(Op)%100-1]
2357	c = SGData['SGCen'][abs(Op)//100]
2358	SyOps.append([inv,m,t,c,unit])
2359	Dist = calcDist(Oatoms,SyOps,Amat)
2360
2361	sig = -0.001
2362	if 'covMatrix' in covData:
2363	dx = .00001
2364	dadx = np.zeros(6)
2365	for i in range(6):
2366	ia = i//3
2367	ix = i%3
2368	Oatoms[ia][ix+1] += dx
2369	a0 = calcDist(Oatoms,SyOps,Amat)
2370	Oatoms[ia][ix+1] -= 2*dx
2371	dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2372	covMatrix = covData['covMatrix']
2373	varyList = covData['varyList']
2374	DistVcov = getVCov(names,varyList,covMatrix)
2375	sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
2376	if sig < 0.001:
2377	sig = -0.001
2378
2379	return Dist,sig
2380
2381	def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
2382	'''default doc string
2383
2384	:param type name: description
2385
2386	:returns: type name: description
2387
2388	'''
2389
2390	def calcAngle(Atoms,SyOps,Amat):
2391	XYZ = []
2392	for i,atom in enumerate(Atoms):
2393	Inv,M,T,C,U = SyOps[i]
2394	XYZ.append(np.array(atom[1:4]))
2395	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2396	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2397	V1 = XYZ[1]-XYZ[0]
2398	V1 /= np.sqrt(np.sum(V1**2))
2399	V2 = XYZ[1]-XYZ[2]
2400	V2 /= np.sqrt(np.sum(V2**2))
2401	V3 = V2-V1
2402	cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2403	return acosd(cang)
2404
2405	SyOps = []
2406	names = []
2407	for i,atom in enumerate(Oatoms):
2408	names += atom[-1]
2409	Op,unit = Atoms[i][-1]
2410	inv = Op//abs(Op)
2411	m,t = SGData['SGOps'][abs(Op)%100-1]
2412	c = SGData['SGCen'][abs(Op)//100]
2413	SyOps.append([inv,m,t,c,unit])
2414	Angle = calcAngle(Oatoms,SyOps,Amat)
2415
2416	sig = -0.01
2417	if 'covMatrix' in covData:
2418	dx = .00001
2419	dadx = np.zeros(9)
2420	for i in range(9):
2421	ia = i//3
2422	ix = i%3
2423	Oatoms[ia][ix+1] += dx
2424	a0 = calcAngle(Oatoms,SyOps,Amat)
2425	Oatoms[ia][ix+1] -= 2*dx
2426	dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2427	covMatrix = covData['covMatrix']
2428	varyList = covData['varyList']
2429	AngVcov = getVCov(names,varyList,covMatrix)
2430	sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2431	if sig < 0.01:
2432	sig = -0.01
2433
2434	return Angle,sig
2435
2436	def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
2437	'''default doc string
2438
2439	:param type name: description
2440
2441	:returns: type name: description
2442
2443	'''
2444
2445	def calcTorsion(Atoms,SyOps,Amat):
2446
2447	XYZ = []
2448	for i,atom in enumerate(Atoms):
2449	Inv,M,T,C,U = SyOps[i]
2450	XYZ.append(np.array(atom[1:4]))
2451	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2452	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2453	V1 = XYZ[1]-XYZ[0]
2454	V2 = XYZ[2]-XYZ[1]
2455	V3 = XYZ[3]-XYZ[2]
2456	V1 /= np.sqrt(np.sum(V1**2))
2457	V2 /= np.sqrt(np.sum(V2**2))
2458	V3 /= np.sqrt(np.sum(V3**2))
2459	M = np.array([V1,V2,V3])
2460	D = nl.det(M)
2461	P12 = np.dot(V1,V2)
2462	P13 = np.dot(V1,V3)
2463	P23 = np.dot(V2,V3)
2464	Tors = acosd((P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23*2)))D/abs(D)
2465	return Tors
2466
2467	SyOps = []
2468	names = []
2469	for i,atom in enumerate(Oatoms):
2470	names += atom[-1]
2471	Op,unit = Atoms[i][-1]
2472	inv = Op//abs(Op)
2473	m,t = SGData['SGOps'][abs(Op)%100-1]
2474	c = SGData['SGCen'][abs(Op)//100]
2475	SyOps.append([inv,m,t,c,unit])
2476	Tors = calcTorsion(Oatoms,SyOps,Amat)
2477
2478	sig = -0.01
2479	if 'covMatrix' in covData:
2480	dx = .00001
2481	dadx = np.zeros(12)
2482	for i in range(12):
2483	ia = i//3
2484	ix = i%3
2485	Oatoms[ia][ix+1] -= dx
2486	a0 = calcTorsion(Oatoms,SyOps,Amat)
2487	Oatoms[ia][ix+1] += 2*dx
2488	dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2489	Oatoms[ia][ix+1] -= dx
2490	covMatrix = covData['covMatrix']
2491	varyList = covData['varyList']
2492	TorVcov = getVCov(names,varyList,covMatrix)
2493	sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
2494	if sig < 0.01:
2495	sig = -0.01
2496
2497	return Tors,sig
2498
2499	def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
2500	'''default doc string
2501
2502	:param type name: description
2503
2504	:returns: type name: description
2505
2506	'''
2507
2508	def calcDist(Atoms,SyOps,Amat):
2509	XYZ = []
2510	for i,atom in enumerate(Atoms):
2511	Inv,M,T,C,U = SyOps[i]
2512	XYZ.append(np.array(atom[1:4]))
2513	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2514	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2515	V1 = XYZ[1]-XYZ[0]
2516	return np.sqrt(np.sum(V1**2))
2517
2518	def calcAngle(Atoms,SyOps,Amat):
2519	XYZ = []
2520	for i,atom in enumerate(Atoms):
2521	Inv,M,T,C,U = SyOps[i]
2522	XYZ.append(np.array(atom[1:4]))
2523	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2524	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2525	V1 = XYZ[1]-XYZ[0]
2526	V1 /= np.sqrt(np.sum(V1**2))
2527	V2 = XYZ[1]-XYZ[2]
2528	V2 /= np.sqrt(np.sum(V2**2))
2529	V3 = V2-V1
2530	cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2531	return acosd(cang)
2532
2533	def calcTorsion(Atoms,SyOps,Amat):
2534
2535	XYZ = []
2536	for i,atom in enumerate(Atoms):
2537	Inv,M,T,C,U = SyOps[i]
2538	XYZ.append(np.array(atom[1:4]))
2539	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2540	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2541	V1 = XYZ[1]-XYZ[0]
2542	V2 = XYZ[2]-XYZ[1]
2543	V3 = XYZ[3]-XYZ[2]
2544	V1 /= np.sqrt(np.sum(V1**2))
2545	V2 /= np.sqrt(np.sum(V2**2))
2546	V3 /= np.sqrt(np.sum(V3**2))
2547	M = np.array([V1,V2,V3])
2548	D = nl.det(M)
2549	P12 = np.dot(V1,V2)
2550	P13 = np.dot(V1,V3)
2551	P23 = np.dot(V2,V3)
2552	Tors = acosd((P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23*2)))D/abs(D)
2553	return Tors
2554
2555	SyOps = []
2556	names = []
2557	for i,atom in enumerate(Oatoms):
2558	names += atom[-1]
2559	Op,unit = Atoms[i][-1]
2560	inv = Op//abs(Op)
2561	m,t = SGData['SGOps'][abs(Op)%100-1]
2562	c = SGData['SGCen'][abs(Op)//100]
2563	SyOps.append([inv,m,t,c,unit])
2564	M = len(Oatoms)
2565	if M == 2:
2566	Val = calcDist(Oatoms,SyOps,Amat)
2567	elif M == 3:
2568	Val = calcAngle(Oatoms,SyOps,Amat)
2569	else:
2570	Val = calcTorsion(Oatoms,SyOps,Amat)
2571
2572	sigVals = [-0.001,-0.01,-0.01]
2573	sig = sigVals[M-3]
2574	if 'covMatrix' in covData:
2575	dx = .00001
2576	N = M*3
2577	dadx = np.zeros(N)
2578	for i in range(N):
2579	ia = i//3
2580	ix = i%3
2581	Oatoms[ia][ix+1] += dx
2582	if M == 2:
2583	a0 = calcDist(Oatoms,SyOps,Amat)
2584	elif M == 3:
2585	a0 = calcAngle(Oatoms,SyOps,Amat)
2586	else:
2587	a0 = calcTorsion(Oatoms,SyOps,Amat)
2588	Oatoms[ia][ix+1] -= 2*dx
2589	if M == 2:
2590	dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2591	elif M == 3:
2592	dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2593	else:
2594	dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2595	covMatrix = covData['covMatrix']
2596	varyList = covData['varyList']
2597	Vcov = getVCov(names,varyList,covMatrix)
2598	sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
2599	if sig < sigVals[M-3]:
2600	sig = sigVals[M-3]
2601
2602	return Val,sig
2603
2604	def ValEsd(value,esd=0,nTZ=False):
2605	'''Format a floating point number with a given level of precision or
2606	with in crystallographic format with a "esd", as value(esd). If esd is
2607	negative the number is formatted with the level of significant figures
2608	appropriate if abs(esd) were the esd, but the esd is not included.
2609	if the esd is zero, approximately 6 significant figures are printed.
2610	nTZ=True causes "extra" zeros to be removed after the decimal place.
2611	for example:
2612
2613	* "1.235(3)" for value=1.2346 & esd=0.003
2614	* "1.235(3)e4" for value=12346. & esd=30
2615	* "1.235(3)e6" for value=0.12346e7 & esd=3000
2616	* "1.235" for value=1.2346 & esd=-0.003
2617	* "1.240" for value=1.2395 & esd=-0.003
2618	* "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
2619	* "1.23460" for value=1.2346 & esd=0.0
2620
2621	:param float value: number to be formatted
2622	:param float esd: uncertainty or if esd < 0, specifies level of
2623	precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
2624	:param bool nTZ: True to remove trailing zeros (default is False)
2625	:returns: value(esd) or value as a string
2626
2627	'''
2628	# Note: this routine is Python 3 compatible -- I think
2629	cutoff = 3.16228 #=(sqrt(10); same as old GSAS was 1.95
2630	if math.isnan(value): # invalid value, bail out
2631	return '?'
2632	if math.isnan(esd): # invalid esd, treat as zero
2633	esd = 0
2634	esdoff = 5
2635	# if esd < 1.e-5:
2636	# esd = 0
2637	# esdoff = 5
2638	elif esd != 0:
2639	# transform the esd to a one or two digit integer
2640	l = math.log10(abs(esd)) % 1.
2641	if l < math.log10(cutoff): l+= 1.
2642	intesd = int(round(10**l)) # esd as integer
2643	# determine the number of digits offset for the esd
2644	esdoff = int(round(math.log10(intesd*1./abs(esd))))
2645	else:
2646	esdoff = 5
2647	valoff = 0
2648	if abs(value) < abs(esdoff): # value is effectively zero
2649	pass
2650	elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
2651	# where the digit offset is to the left of the decimal place or where too many
2652	# digits are needed
2653	if abs(value) > 1:
2654	valoff = int(math.log10(abs(value)))
2655	elif abs(value) > 0:
2656	valoff = int(math.log10(abs(value))-0.9999999)
2657	else:
2658	valoff = 0
2659	if esd != 0:
2660	if valoff+esdoff < 0:
2661	valoff = esdoff = 0
2662	out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
2663	elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
2664	out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
2665	else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
2666	if abs(value) > 0:
2667	extra = -math.log10(abs(value))
2668	else:
2669	extra = 0
2670	if extra > 0: extra += 1
2671	out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
2672	if esd > 0:
2673	out += ("({:d})").format(intesd) # add the esd
2674	elif nTZ and '.' in out:
2675	out = out.rstrip('0') # strip zeros to right of decimal
2676	out = out.rstrip('.') # and decimal place when not needed
2677	if valoff != 0:
2678	out += ("e{:d}").format(valoff) # add an exponent, when needed
2679	return out
2680
2681	###############################################################################
2682	##### Protein validation - "ERRATV2" analysis
2683	###############################################################################
2684
2685	def validProtein(Phase,old):
2686
2687	def sumintact(intact):
2688	return {'CC':intact['CC'],'NN':intact['NN'],'OO':intact['OO'],
2689	'CN':(intact['CN']+intact['NC']),'CO':(intact['CO']+intact['OC']),
2690	'NO':(intact['NO']+intact['ON'])}
2691
2692	resNames = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
2693	'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
2694	# data from errat.f
2695	b1_old = np.array([
2696	[1154.343, 600.213, 1051.018, 1132.885, 960.738],
2697	[600.213, 1286.818, 1282.042, 957.156, 612.789],
2698	[1051.018, 1282.042, 3519.471, 991.974, 1226.491],
2699	[1132.885, 957.156, 991.974, 1798.672, 820.355],
2700	[960.738, 612.789, 1226.491, 820.355, 2428.966]
2701	])
2702	avg_old = np.array([ 0.225, 0.281, 0.071, 0.237, 0.044]) #Table 1 3.5A Obsd. Fr. p 1513
2703	# data taken from erratv2.ccp
2704	b1 = np.array([
2705	[5040.279078850848200, 3408.805141583649400, 4152.904423767300600, 4236.200004171890200, 5054.781210204625500],
2706	[3408.805141583648900, 8491.906094010220800, 5958.881777877950300, 1521.387352718486200, 4304.078200827221700],
2707	[4152.904423767301500, 5958.881777877952100, 7637.167089335050100, 6620.715738223072500, 5287.691183798410700],
2708	[4236.200004171890200, 1521.387352718486200, 6620.715738223072500, 18368.343774298410000, 4050.797811118806700],
2709	[5054.781210204625500, 4304.078200827220800, 5287.691183798409800, 4050.797811118806700, 6666.856740479164700]])
2710	avg = np.array([0.192765509919262, 0.195575208778518, 0.275322406824210, 0.059102357035642, 0.233154192767480])
2711	General = Phase['General']
2712	Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
2713	cx,ct,cs,cia = General['AtomPtrs']
2714	Atoms = Phase['Atoms']
2715	cartAtoms = []
2716	xyzmin = 999.*np.ones(3)
2717	xyzmax = -999.*np.ones(3)
2718	#select residue atoms,S,Se --> O make cartesian
2719	for atom in Atoms:
2720	if atom[1] in resNames:
2721	cartAtoms.append(atom[:cx+3])
2722	if atom[4].strip() in ['S','Se']:
2723	if not old:
2724	continue #S,Se skipped for erratv2?
2725	cartAtoms[-1][3] = 'Os'
2726	cartAtoms[-1][4] = 'O'
2727	cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
2728	cartAtoms[-1].append(atom[cia+8])
2729	XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
2730	xyzmin = np.array([np.min(XYZ.T[i]) for i in [0,1,2]])
2731	xyzmax = np.array([np.max(XYZ.T[i]) for i in [0,1,2]])
2732	nbox = list(np.array(np.ceil((xyzmax-xyzmin)/4.),dtype=int))+[15,]
2733	Boxes = np.zeros(nbox,dtype=int)
2734	iBox = np.array([np.trunc((XYZ.T[i]-xyzmin[i])/4.) for i in [0,1,2]],dtype=int).T
2735	for ib,box in enumerate(iBox): #put in a try for too many atoms in box (IndexError)?
2736	try:
2737	Boxes[box[0],box[1],box[2],0] += 1
2738	Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
2739	except IndexError:
2740	G2fil.G2Print('Error: too many atoms in box' )
2741	continue
2742	#Box content checks with errat.f $ erratv2.cpp ibox1 arrays
2743	indices = (-1,0,1)
2744	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2745	dsmax = 3.75**2
2746	if old:
2747	dsmax = 3.5**2
2748	chains = []
2749	resIntAct = []
2750	chainIntAct = []
2751	res = []
2752	resNames = []
2753	resIDs = {}
2754	resname = []
2755	resID = {}
2756	newChain = True
2757	intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2758	for ia,atom in enumerate(cartAtoms):
2759	jntact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2760	if atom[2] not in chains: #get chain id & save residue sequence from last chain
2761	chains.append(atom[2])
2762	if len(resIntAct):
2763	resIntAct.append(sumintact(intact))
2764	chainIntAct.append(resIntAct)
2765	resNames += resname
2766	resIDs.update(resID)
2767	res = []
2768	resname = []
2769	resID = {}
2770	resIntAct = []
2771	intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2772	newChain = True
2773	if atom[0] not in res: #new residue, get residue no.
2774	if res and int(res[-1]) != int(atom[0])-1: #a gap in chain - not new chain
2775	intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2776	ires = int(res[-1])
2777	for i in range(int(atom[0])-ires-1):
2778	res.append(str(ires+i+1))
2779	resname.append('')
2780	resIntAct.append(sumintact(intact))
2781	res.append(atom[0])
2782	name = '%s-%s%s'%(atom[2],atom[0],atom[1])
2783	resname.append(name)
2784	resID[name] = atom[-1]
2785	if not newChain:
2786	resIntAct.append(sumintact(intact))
2787	intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2788	newChain = False
2789	ibox = iBox[ia] #box location of atom
2790	tgts = []
2791	for unit in Units: #assemble list of all possible target atoms
2792	jbox = ibox+unit
2793	if np.all(jbox>=0) and np.all((jbox-nbox[:3])<0):
2794	tgts += list(Boxes[jbox[0],jbox[1],jbox[2]])
2795	tgts = list(set(tgts))
2796	tgts = [tgt for tgt in tgts if atom[:3] != cartAtoms[tgt][:3]] #exclude same residue
2797	tgts = [tgt for tgt in tgts if np.sum((XYZ[ia]-XYZ[tgt])**2) < dsmax]
2798	ires = int(atom[0])
2799	if old:
2800	if atom[3].strip() == 'C':
2801	tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2802	elif atom[3].strip() == 'N':
2803	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])]
2804	elif atom[3].strip() == 'CA':
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	else:
2807	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]]
2808	if atom[3].strip() == 'C':
2809	tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) == ires+1)]
2810	elif atom[3].strip() == 'N':
2811	tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'C' and int(cartAtoms[tgt][0]) == ires-1)]
2812	for tgt in tgts:
2813	dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
2814	mult = 1.0
2815	if dsqt > 3.25 and not old:
2816	mult = 2.*(3.75-dsqt)
2817	intype = atom[4].strip()+cartAtoms[tgt][4].strip()
2818	if 'S' not in intype:
2819	intact[intype] += mult
2820	jntact[intype] += mult
2821	# 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']
2822	resNames += resname
2823	resIDs.update(resID)
2824	resIntAct.append(sumintact(intact))
2825	chainIntAct.append(resIntAct)
2826	chainProb = []
2827	for ich,chn in enumerate(chains):
2828	IntAct = chainIntAct[ich]
2829	nRes = len(IntAct)
2830	Probs = [0.,0.,0.,0.] #skip 1st 4 residues in chain
2831	for i in range(4,nRes-4):
2832	if resNames[i]:
2833	mtrx = np.zeros(5)
2834	summ = 0.
2835	for j in range(i-4,i+5):
2836	summ += np.sum(np.array(list(IntAct[j].values())))
2837	if old:
2838	mtrx[0] += IntAct[j]['CC']
2839	mtrx[1] += IntAct[j]['CO']
2840	mtrx[2] += IntAct[j]['NN']
2841	mtrx[3] += IntAct[j]['NO']
2842	mtrx[4] += IntAct[j]['OO']
2843	else:
2844	mtrx[0] += IntAct[j]['CC']
2845	mtrx[1] += IntAct[j]['CN']
2846	mtrx[2] += IntAct[j]['CO']
2847	mtrx[3] += IntAct[j]['NN']
2848	mtrx[4] += IntAct[j]['NO']
2849	mtrx /= summ
2850	# print i+1,mtrx*summ
2851	if old:
2852	mtrx -= avg_old
2853	prob = np.inner(np.inner(mtrx,b1_old),mtrx)
2854	else:
2855	mtrx -= avg
2856	prob = np.inner(np.inner(mtrx,b1),mtrx)
2857	else: #skip the gaps
2858	prob = 0.0
2859	Probs.append(prob)
2860	Probs += 4*[0.,] #skip last 4 residues in chain
2861	chainProb += Probs
2862	return resNames,chainProb,resIDs
2863
2864	################################################################################
2865	##### Texture fitting stuff
2866	################################################################################
2867
2868	def FitTexture(General,Gangls,refData,keyList,pgbar):
2869	import pytexture as ptx
2870	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
2871
2872	def printSpHarm(textureData,SHtextureSig):
2873	print ('\n Spherical harmonics texture: Order:' + str(textureData['Order']))
2874	names = ['omega','chi','phi']
2875	namstr = ' names :'
2876	ptstr = ' values:'
2877	sigstr = ' esds :'
2878	for name in names:
2879	namstr += '%12s'%('Sample '+name)
2880	ptstr += '%12.3f'%(textureData['Sample '+name][1])
2881	if 'Sample '+name in SHtextureSig:
2882	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2883	else:
2884	sigstr += 12*' '
2885	print (namstr)
2886	print (ptstr)
2887	print (sigstr)
2888	print ('\n Texture coefficients:')
2889	SHcoeff = textureData['SH Coeff'][1]
2890	SHkeys = list(SHcoeff.keys())
2891	nCoeff = len(SHcoeff)
2892	nBlock = nCoeff//10+1
2893	iBeg = 0
2894	iFin = min(iBeg+10,nCoeff)
2895	for block in range(nBlock):
2896	namstr = ' names :'
2897	ptstr = ' values:'
2898	sigstr = ' esds :'
2899	for name in SHkeys[iBeg:iFin]:
2900	if 'C' in name:
2901	namstr += '%12s'%(name)
2902	ptstr += '%12.3f'%(SHcoeff[name])
2903	if name in SHtextureSig:
2904	sigstr += '%12.3f'%(SHtextureSig[name])
2905	else:
2906	sigstr += 12*' '
2907	print (namstr)
2908	print (ptstr)
2909	print (sigstr)
2910	iBeg += 10
2911	iFin = min(iBeg+10,nCoeff)
2912
2913	def Dict2Values(parmdict, varylist):
2914	'''Use before call to leastsq to setup list of values for the parameters
2915	in parmdict, as selected by key in varylist'''
2916	return [parmdict[key] for key in varylist]
2917
2918	def Values2Dict(parmdict, varylist, values):
2919	''' Use after call to leastsq to update the parameter dictionary with
2920	values corresponding to keys in varylist'''
2921	parmdict.update(list(zip(varylist,values)))
2922
2923	def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2924	parmDict.update(list(zip(varyList,values)))
2925	Mat = np.empty(0)
2926	sumObs = 0
2927	Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2928	for hist in Gangls.keys():
2929	Refs = refData[hist]
2930	Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2931	wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2932	# wt = 1./np.max(Refs[:,4],.25)
2933	sumObs += np.sum(wt*Refs[:,5])
2934	Refs[:,6] = 1.
2935	H = Refs[:,:3]
2936	phi,beta = G2lat.CrsAng(H,cell,SGData)
2937	psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2938	for item in parmDict:
2939	if 'C' in item:
2940	L,M,N = eval(item.strip('C'))
2941	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2942	Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2943	Lnorm = G2lat.Lnorm(L)
2944	Refs[:,6] += parmDict[item]LnormKcl*Ksl
2945	mat = wt*(Refs[:,5]-Refs[:,6])
2946	Mat = np.concatenate((Mat,mat))
2947	sumD = np.sum(np.abs(Mat))
2948	R = min(100.,100.*sumD/sumObs)
2949	pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2950	print (' Residual: %.3f%%'%(R))
2951	return Mat
2952
2953	def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2954	Mat = np.empty(0)
2955	Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2956	for hist in Gangls.keys():
2957	mat = np.zeros((len(varyList),len(refData[hist])))
2958	Refs = refData[hist]
2959	H = Refs[:,:3]
2960	wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2961	# wt = 1./np.max(Refs[:,4],.25)
2962	phi,beta = G2lat.CrsAng(H,cell,SGData)
2963	psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2964	for j,item in enumerate(varyList):
2965	if 'C' in item:
2966	L,M,N = eval(item.strip('C'))
2967	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2968	Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2969	Lnorm = G2lat.Lnorm(L)
2970	mat[j] = -wtLnormKcl*Ksl
2971	for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2972	try:
2973	l = varyList.index(itema)
2974	mat[l] -= parmDict[item]wtLnormKcl(dKdpdPdA[k]+dKdgdGdA[k])
2975	except ValueError:
2976	pass
2977	if len(Mat):
2978	Mat = np.concatenate((Mat,mat.T))
2979	else:
2980	Mat = mat.T
2981	print ('deriv')
2982	return Mat
2983
2984	print (' Fit texture for '+General['Name'])
2985	SGData = General['SGData']
2986	cell = General['Cell'][1:7]
2987	Texture = General['SH Texture']
2988	if not Texture['Order']:
2989	return 'No spherical harmonics coefficients'
2990	varyList = []
2991	parmDict = copy.copy(Texture['SH Coeff'][1])
2992	for item in ['Sample omega','Sample chi','Sample phi']:
2993	parmDict[item] = Texture[item][1]
2994	if Texture[item][0]:
2995	varyList.append(item)
2996	if Texture['SH Coeff'][0]:
2997	varyList += list(Texture['SH Coeff'][1].keys())
2998	while True:
2999	begin = time.time()
3000	values = np.array(Dict2Values(parmDict, varyList))
3001	result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
3002	args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
3003	ncyc = int(result[2]['nfev']//2)
3004	if ncyc:
3005	runtime = time.time()-begin
3006	chisq = np.sum(result[2]['fvec']**2)
3007	Values2Dict(parmDict, varyList, result[0])
3008	GOF = chisq/(len(result[2]['fvec'])-len(varyList)) #reduced chi^2
3009	G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(result[2]['fvec']),len(varyList)))
3010	G2fil.G2Print ('refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
3011	try:
3012	sig = np.sqrt(np.diag(result[1])*GOF)
3013	if np.any(np.isnan(sig)):
3014	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *', mode='error')
3015	break #refinement succeeded - finish up!
3016	except ValueError: #result[1] is None on singular matrix
3017	G2fil.G2Print ('** Refinement failed - singular matrix **', mode='error')
3018	return None
3019	else:
3020	break
3021
3022	if ncyc:
3023	for parm in parmDict:
3024	if 'C' in parm:
3025	Texture['SH Coeff'][1][parm] = parmDict[parm]
3026	else:
3027	Texture[parm][1] = parmDict[parm]
3028	sigDict = dict(zip(varyList,sig))
3029	printSpHarm(Texture,sigDict)
3030
3031	return None
3032
3033	################################################################################
3034	##### Fourier & charge flip stuff
3035	################################################################################
3036
3037	def adjHKLmax(SGData,Hmax):
3038	'''default doc string
3039
3040	:param type name: description
3041
3042	:returns: type name: description
3043
3044	'''
3045	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3046	Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
3047	Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
3048	Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
3049	else:
3050	Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
3051	Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
3052	Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
3053
3054	def OmitMap(data,reflDict,pgbar=None):
3055	'''default doc string
3056
3057	:param type name: description
3058
3059	:returns: type name: description
3060
3061	'''
3062	generalData = data['General']
3063	if not generalData['Map']['MapType']:
3064	G2fil.G2Print ('**** ERROR - Fourier map not defined')
3065	return
3066	mapData = generalData['Map']
3067	dmin = mapData['GridStep']*2.
3068	SGData = generalData['SGData']
3069	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3070	SGT = np.array([ops[1] for ops in SGData['SGOps']])
3071	cell = generalData['Cell'][1:8]
3072	A = G2lat.cell2A(cell[:6])
3073	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3074	adjHKLmax(SGData,Hmax)
3075	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3076	time0 = time.time()
3077	for iref,ref in enumerate(reflDict['RefList']):
3078	if ref[4] >= dmin:
3079	Fosq,Fcsq,ph = ref[8:11]
3080	Uniq = np.inner(ref[:3],SGMT)
3081	Phi = np.inner(ref[:3],SGT)
3082	for i,hkl in enumerate(Uniq): #uses uniq
3083	hkl = np.asarray(hkl,dtype='i')
3084	dp = 360.*Phi[i] #and phi
3085	a = cosd(ph+dp)
3086	b = sind(ph+dp)
3087	phasep = complex(a,b)
3088	phasem = complex(a,-b)
3089	if '2Fo-Fc' in mapData['MapType']:
3090	F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
3091	else:
3092	F = np.sqrt(Fosq)
3093	h,k,l = hkl+Hmax
3094	Fhkl[h,k,l] = F*phasep
3095	h,k,l = -hkl+Hmax
3096	Fhkl[h,k,l] = F*phasem
3097	rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3098	M = np.mgrid[0:4,0:4,0:4]
3099	blkIds = np.array(list(zip(M[0].flatten(),M[1].flatten(),M[2].flatten())))
3100	iBeg = blkIds*rho0.shape//4
3101	iFin = (blkIds+1)*rho0.shape//4
3102	rho_omit = np.zeros_like(rho0)
3103	nBlk = 0
3104	for iB,iF in zip(iBeg,iFin):
3105	rho1 = np.copy(rho0)
3106	rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
3107	Fnew = fft.ifftshift(fft.ifftn(rho1))
3108	Fnew = np.where(Fnew,Fnew,1.0) #avoid divide by zero
3109	phase = Fnew/np.absolute(Fnew)
3110	OFhkl = np.absolute(Fhkl)*phase
3111	rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
3112	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]])
3113	nBlk += 1
3114	pgbar.Update(nBlk)
3115	mapData['rho'] = np.real(rho_omit)/cell[6]
3116	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3117	mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3118	G2fil.G2Print ('Omit map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))
3119	return mapData
3120
3121	def FourierMap(data,reflDict):
3122	'''default doc string
3123
3124	:param type name: description
3125
3126	:returns: type name: description
3127
3128	'''
3129	generalData = data['General']
3130	mapData = generalData['Map']
3131	dmin = mapData['GridStep']*2.
3132	SGData = generalData['SGData']
3133	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3134	SGT = np.array([ops[1] for ops in SGData['SGOps']])
3135	cell = generalData['Cell'][1:8]
3136	A = G2lat.cell2A(cell[:6])
3137	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3138	adjHKLmax(SGData,Hmax)
3139	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3140	# Fhkl[0,0,0] = generalData['F000X']
3141	time0 = time.time()
3142	for iref,ref in enumerate(reflDict['RefList']):
3143	if ref[4] > dmin:
3144	Fosq,Fcsq,ph = ref[8:11]
3145	Uniq = np.inner(ref[:3],SGMT)
3146	Phi = np.inner(ref[:3],SGT)
3147	for i,hkl in enumerate(Uniq): #uses uniq
3148	hkl = np.asarray(hkl,dtype='i')
3149	dp = 360.*Phi[i] #and phi
3150	a = cosd(ph+dp)
3151	b = sind(ph+dp)
3152	phasep = complex(a,b)
3153	phasem = complex(a,-b)
3154	if 'Fobs' in mapData['MapType']:
3155	F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
3156	h,k,l = hkl+Hmax
3157	Fhkl[h,k,l] = F*phasep
3158	h,k,l = -hkl+Hmax
3159	Fhkl[h,k,l] = F*phasem
3160	elif 'Fcalc' in mapData['MapType']:
3161	F = np.sqrt(Fcsq)
3162	h,k,l = hkl+Hmax
3163	Fhkl[h,k,l] = F*phasep
3164	h,k,l = -hkl+Hmax
3165	Fhkl[h,k,l] = F*phasem
3166	elif 'delt-F' in mapData['MapType']:
3167	dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3168	h,k,l = hkl+Hmax
3169	Fhkl[h,k,l] = dF*phasep
3170	h,k,l = -hkl+Hmax
3171	Fhkl[h,k,l] = dF*phasem
3172	elif '2*Fo-Fc' in mapData['MapType']:
3173	F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3174	h,k,l = hkl+Hmax
3175	Fhkl[h,k,l] = F*phasep
3176	h,k,l = -hkl+Hmax
3177	Fhkl[h,k,l] = F*phasem
3178	elif 'Patterson' in mapData['MapType']:
3179	h,k,l = hkl+Hmax
3180	Fhkl[h,k,l] = complex(Fosq,0.)
3181	h,k,l = -hkl+Hmax
3182	Fhkl[h,k,l] = complex(Fosq,0.)
3183	rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3184	G2fil.G2Print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))
3185	mapData['Type'] = reflDict['Type']
3186	mapData['rho'] = np.real(rho)
3187	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3188	mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3189
3190	def Fourier4DMap(data,reflDict):
3191	'''default doc string
3192
3193	:param type name: description
3194
3195	:returns: type name: description
3196
3197	'''
3198	generalData = data['General']
3199	map4DData = generalData['4DmapData']
3200	mapData = generalData['Map']
3201	dmin = mapData['GridStep']*2.
3202	SGData = generalData['SGData']
3203	SSGData = generalData['SSGData']
3204	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3205	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3206	cell = generalData['Cell'][1:8]
3207	A = G2lat.cell2A(cell[:6])
3208	maxM = 4
3209	Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
3210	adjHKLmax(SGData,Hmax)
3211	Hmax = np.asarray(Hmax,dtype='i')+1
3212	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3213	time0 = time.time()
3214	for iref,ref in enumerate(reflDict['RefList']):
3215	if ref[5] > dmin:
3216	Fosq,Fcsq,ph = ref[9:12]
3217	Fosq = np.where(Fosq>0.,Fosq,0.) #can't use Fo^2 < 0
3218	Uniq = np.inner(ref[:4],SSGMT)
3219	Phi = np.inner(ref[:4],SSGT)
3220	for i,hkl in enumerate(Uniq): #uses uniq
3221	hkl = np.asarray(hkl,dtype='i')
3222	dp = 360.*Phi[i] #and phi
3223	a = cosd(ph+dp)
3224	b = sind(ph+dp)
3225	phasep = complex(a,b)
3226	phasem = complex(a,-b)
3227	if 'Fobs' in mapData['MapType']:
3228	F = np.sqrt(Fosq)
3229	h,k,l,m = hkl+Hmax
3230	Fhkl[h,k,l,m] = F*phasep
3231	h,k,l,m = -hkl+Hmax
3232	Fhkl[h,k,l,m] = F*phasem
3233	elif 'Fcalc' in mapData['MapType']:
3234	F = np.sqrt(Fcsq)
3235	h,k,l,m = hkl+Hmax
3236	Fhkl[h,k,l,m] = F*phasep
3237	h,k,l,m = -hkl+Hmax
3238	Fhkl[h,k,l,m] = F*phasem
3239	elif 'delt-F' in mapData['MapType']:
3240	dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
3241	h,k,l,m = hkl+Hmax
3242	Fhkl[h,k,l,m] = dF*phasep
3243	h,k,l,m = -hkl+Hmax
3244	Fhkl[h,k,l,m] = dF*phasem
3245	SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6] #4D map
3246	rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6] #3D map
3247	map4DData['rho'] = np.real(SSrho)
3248	map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3249	map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3250	map4DData['Type'] = reflDict['Type']
3251	mapData['Type'] = reflDict['Type']
3252	mapData['rho'] = np.real(rho)
3253	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3254	mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3255	G2fil.G2Print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))
3256
3257	# map printing for testing purposes
3258	def printRho(SGLaue,rho,rhoMax):
3259	'''default doc string
3260
3261	:param type name: description
3262
3263	:returns: type name: description
3264
3265	'''
3266	dim = len(rho.shape)
3267	if dim == 2:
3268	ix,jy = rho.shape
3269	for j in range(jy):
3270	line = ''
3271	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3272	line += (jy-j)*' '
3273	for i in range(ix):
3274	r = int(100*rho[i,j]/rhoMax)
3275	line += '%4d'%(r)
3276	print (line+'\n')
3277	else:
3278	ix,jy,kz = rho.shape
3279	for k in range(kz):
3280	print ('k = %d'%k)
3281	for j in range(jy):
3282	line = ''
3283	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3284	line += (jy-j)*' '
3285	for i in range(ix):
3286	r = int(100*rho[i,j,k]/rhoMax)
3287	line += '%4d'%(r)
3288	print (line+'\n')
3289	## keep this
3290
3291	def findOffset(SGData,A,Fhkl):
3292	'''default doc string
3293
3294	:param type name: description
3295
3296	:returns: type name: description
3297
3298	'''
3299	if SGData['SpGrp'] == 'P 1':
3300	return [0,0,0]
3301	hklShape = Fhkl.shape
3302	hklHalf = np.array(hklShape)/2
3303	sortHKL = np.argsort(Fhkl.flatten())
3304	Fdict = {}
3305	for hkl in sortHKL:
3306	HKL = np.unravel_index(hkl,hklShape)
3307	F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
3308	if F == 0.:
3309	break
3310	Fdict['%.6f'%(np.absolute(F))] = hkl
3311	Flist = np.flipud(np.sort(list(Fdict.keys())))
3312	F = str(1.e6)
3313	i = 0
3314	DH = []
3315	Dphi = []
3316	Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3317	for F in Flist:
3318	hkl = np.unravel_index(Fdict[F],hklShape)
3319	if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
3320	continue
3321	iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
3322	Uniq = np.array(Uniq,dtype='i')
3323	Phi = np.array(Phi)
3324	Uniq = np.concatenate((Uniq,-Uniq))+hklHalf # put in Friedel pairs & make as index to Farray
3325	Phi = np.concatenate((Phi,-Phi)) # and their phase shifts
3326	Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
3327	ang0 = np.angle(Fh0,deg=True)/360.
3328	for H,phi in list(zip(Uniq,Phi))[1:]:
3329	ang = (np.angle(Fhkl[int(H[0]),int(H[1]),int(H[2])],deg=True)/360.-phi)
3330	dH = H-hkl
3331	dang = ang-ang0
3332	DH.append(dH)
3333	Dphi.append((dang+.5) % 1.0)
3334	if i > 20 or len(DH) > 30:
3335	break
3336	i += 1
3337	DH = np.array(DH)
3338	G2fil.G2Print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))
3339	Dphi = np.array(Dphi)
3340	steps = np.array(hklShape)
3341	X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
3342	XYZ = np.array(list(zip(X.flatten(),Y.flatten(),Z.flatten())))
3343	Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
3344	Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3345	hist,bins = np.histogram(Mmap,bins=1000)
3346	# for i,item in enumerate(hist[:10]):
3347	# print item,bins[i]
3348	chisq = np.min(Mmap)
3349	DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3350	G2fil.G2Print (' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2]))
3351	# print (np.dot(DX,DH.T)+.5)%1.-Dphi
3352	return DX
3353
3354	def ChargeFlip(data,reflDict,pgbar):
3355	'''default doc string
3356
3357	:param type name: description
3358
3359	:returns: type name: description
3360
3361	'''
3362	generalData = data['General']
3363	mapData = generalData['Map']
3364	flipData = generalData['Flip']
3365	FFtable = {}
3366	if 'None' not in flipData['Norm element']:
3367	normElem = flipData['Norm element'].upper()
3368	FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3369	for ff in FFs:
3370	if ff['Symbol'] == normElem:
3371	FFtable.update(ff)
3372	dmin = flipData['GridStep']*2.
3373	SGData = generalData['SGData']
3374	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3375	SGT = np.array([ops[1] for ops in SGData['SGOps']])
3376	cell = generalData['Cell'][1:8]
3377	A = G2lat.cell2A(cell[:6])
3378	Vol = cell[6]
3379	im = 0
3380	if generalData['Modulated'] == True:
3381	im = 1
3382	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3383	adjHKLmax(SGData,Hmax)
3384	Ehkl = np.zeros(shape=2*Hmax,dtype='c16') #2X64bits per complex no.
3385	time0 = time.time()
3386	for iref,ref in enumerate(reflDict['RefList']):
3387	dsp = ref[4+im]
3388	if im and ref[3]: #skip super lattice reflections - result is 3D projection
3389	continue
3390	if dsp > dmin:
3391	ff = 0.1Vol #est. no. atoms for ~10A*3/atom
3392	if FFtable:
3393	SQ = 0.25/dsp**2
3394	ff *= G2el.ScatFac(FFtable,SQ)[0]
3395	if ref[8+im] > 0.: #use only +ve Fobs**2
3396	E = np.sqrt(ref[8+im])/ff
3397	else:
3398	E = 0.
3399	ph = ref[10]
3400	ph = rn.uniform(0.,360.)
3401	Uniq = np.inner(ref[:3],SGMT)
3402	Phi = np.inner(ref[:3],SGT)
3403	for i,hkl in enumerate(Uniq): #uses uniq
3404	hkl = np.asarray(hkl,dtype='i')
3405	dp = 360.*Phi[i] #and phi
3406	a = cosd(ph+dp)
3407	b = sind(ph+dp)
3408	phasep = complex(a,b)
3409	phasem = complex(a,-b)
3410	h,k,l = hkl+Hmax
3411	Ehkl[h,k,l] = E*phasep
3412	h,k,l = -hkl+Hmax
3413	Ehkl[h,k,l] = E*phasem
3414	# Ehkl[Hmax] = 0.00001 #this to preserve F[0,0,0]
3415	testHKL = np.array(flipData['testHKL'])+Hmax
3416	CEhkl = copy.copy(Ehkl)
3417	MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3418	Emask = ma.getmask(MEhkl)
3419	sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3420	Ncyc = 0
3421	old = np.seterr(all='raise')
3422	twophases = []
3423	while True:
3424	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3425	CEsig = np.std(CErho)
3426	CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3427	CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho) #solves U atom problem!
3428	CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3429	CFhkl = np.where(CFhkl,CFhkl,1.0) #avoid divide by zero
3430	phase = CFhkl/np.absolute(CFhkl)
3431	twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
3432	CEhkl = np.absolute(Ehkl)*phase
3433	Ncyc += 1
3434	sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3435	DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3436	Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3437	if Rcf < 5.:
3438	break
3439	GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3440	if not GoOn or Ncyc > 10000:
3441	break
3442	np.seterr(**old)
3443	G2fil.G2Print (' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size))
3444	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10. #? to get on same scale as e-map
3445	G2fil.G2Print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))
3446	roll = findOffset(SGData,A,CEhkl) #CEhkl needs to be just the observed set, not the full set!
3447
3448	mapData['Rcf'] = Rcf
3449	mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3450	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3451	mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3452	mapData['Type'] = reflDict['Type']
3453	return mapData,twophases
3454
3455	def findSSOffset(SGData,SSGData,A,Fhklm):
3456	'''default doc string
3457
3458	:param type name: description
3459
3460	:returns: type name: description
3461
3462	'''
3463	if SGData['SpGrp'] == 'P 1':
3464	return [0,0,0,0]
3465	hklmShape = Fhklm.shape
3466	hklmHalf = np.array(hklmShape)/2
3467	sortHKLM = np.argsort(Fhklm.flatten())
3468	Fdict = {}
3469	for hklm in sortHKLM:
3470	HKLM = np.unravel_index(hklm,hklmShape)
3471	F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
3472	if F == 0.:
3473	break
3474	Fdict['%.6f'%(np.absolute(F))] = hklm
3475	Flist = np.flipud(np.sort(list(Fdict.keys())))
3476	F = str(1.e6)
3477	i = 0
3478	DH = []
3479	Dphi = []
3480	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3481	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3482	Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3483	for F in Flist:
3484	hklm = np.unravel_index(Fdict[F],hklmShape)
3485	if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
3486	continue
3487	Uniq = np.inner(hklm-hklmHalf,SSGMT)
3488	Phi = np.inner(hklm-hklmHalf,SSGT)
3489	Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf # put in Friedel pairs & make as index to Farray
3490	Phi = np.concatenate((Phi,-Phi)) # and their phase shifts
3491	Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
3492	ang0 = np.angle(Fh0,deg=True)/360.
3493	for H,phi in list(zip(Uniq,Phi))[1:]:
3494	H = np.array(H,dtype=int)
3495	ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
3496	dH = H-hklm
3497	dang = ang-ang0
3498	DH.append(dH)
3499	Dphi.append((dang+.5) % 1.0)
3500	if i > 20 or len(DH) > 30:
3501	break
3502	i += 1
3503	DH = np.array(DH)
3504	G2fil.G2Print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))
3505	Dphi = np.array(Dphi)
3506	steps = np.array(hklmShape)
3507	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]]
3508	XYZT = np.array(list(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten())))
3509	Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
3510	Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3511	hist,bins = np.histogram(Mmap,bins=1000)
3512	# for i,item in enumerate(hist[:10]):
3513	# print item,bins[i]
3514	chisq = np.min(Mmap)
3515	DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3516	G2fil.G2Print (' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3]))
3517	# print (np.dot(DX,DH.T)+.5)%1.-Dphi
3518	return DX
3519
3520	def SSChargeFlip(data,reflDict,pgbar):
3521	'''default doc string
3522
3523	:param type name: description
3524
3525	:returns: type name: description
3526
3527	'''
3528	generalData = data['General']
3529	mapData = generalData['Map']
3530	map4DData = {}
3531	flipData = generalData['Flip']
3532	FFtable = {}
3533	if 'None' not in flipData['Norm element']:
3534	normElem = flipData['Norm element'].upper()
3535	FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3536	for ff in FFs:
3537	if ff['Symbol'] == normElem:
3538	FFtable.update(ff)
3539	dmin = flipData['GridStep']*2.
3540	SGData = generalData['SGData']
3541	SSGData = generalData['SSGData']
3542	SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3543	SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3544	cell = generalData['Cell'][1:8]
3545	A = G2lat.cell2A(cell[:6])
3546	Vol = cell[6]
3547	maxM = 4
3548	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3549	adjHKLmax(SGData,Hmax)
3550	Ehkl = np.zeros(shape=2*Hmax,dtype='c16') #2X64bits per complex no.
3551	time0 = time.time()
3552	for iref,ref in enumerate(reflDict['RefList']):
3553	dsp = ref[5]
3554	if dsp > dmin:
3555	ff = 0.1Vol #est. no. atoms for ~10A*3/atom
3556	if FFtable:
3557	SQ = 0.25/dsp**2
3558	ff *= G2el.ScatFac(FFtable,SQ)[0]
3559	if ref[9] > 0.: #use only +ve Fobs**2
3560	E = np.sqrt(ref[9])/ff
3561	else:
3562	E = 0.
3563	ph = ref[11]
3564	ph = rn.uniform(0.,360.)
3565	Uniq = np.inner(ref[:4],SSGMT)
3566	Phi = np.inner(ref[:4],SSGT)
3567	for i,hklm in enumerate(Uniq): #uses uniq
3568	hklm = np.asarray(hklm,dtype='i')
3569	dp = 360.*Phi[i] #and phi
3570	a = cosd(ph+dp)
3571	b = sind(ph+dp)
3572	phasep = complex(a,b)
3573	phasem = complex(a,-b)
3574	h,k,l,m = hklm+Hmax
3575	Ehkl[h,k,l,m] = E*phasep
3576	h,k,l,m = -hklm+Hmax #Friedel pair refl.
3577	Ehkl[h,k,l,m] = E*phasem
3578	# Ehkl[Hmax] = 0.00001 #this to preserve F[0,0,0]
3579	CEhkl = copy.copy(Ehkl)
3580	MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3581	Emask = ma.getmask(MEhkl)
3582	sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3583	Ncyc = 0
3584	old = np.seterr(all='raise')
3585	while True:
3586	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3587	CEsig = np.std(CErho)
3588	CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3589	CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho) #solves U atom problem!
3590	CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3591	CFhkl = np.where(CFhkl,CFhkl,1.0) #avoid divide by zero
3592	phase = CFhkl/np.absolute(CFhkl)
3593	CEhkl = np.absolute(Ehkl)*phase
3594	Ncyc += 1
3595	sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3596	DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3597	Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3598	if Rcf < 5.:
3599	break
3600	GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3601	if not GoOn or Ncyc > 10000:
3602	break
3603	np.seterr(**old)
3604	G2fil.G2Print (' Charge flip time: %.4f no. elements: %d'%(time.time()-time0,Ehkl.size))
3605	CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10. #? to get on same scale as e-map
3606	SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10. #? ditto
3607	G2fil.G2Print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))
3608	roll = findSSOffset(SGData,SSGData,A,CEhkl) #CEhkl needs to be just the observed set, not the full set!
3609
3610	mapData['Rcf'] = Rcf
3611	mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3612	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3613	mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3614	mapData['Type'] = reflDict['Type']
3615
3616	map4DData['Rcf'] = Rcf
3617	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))
3618	map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3619	map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3620	map4DData['Type'] = reflDict['Type']
3621	return mapData,map4DData
3622
3623	def getRho(xyz,mapData):
3624	''' get scattering density at a point by 8-point interpolation
3625	param xyz: coordinate to be probed
3626	param: mapData: dict of map data
3627
3628	:returns: density at xyz
3629	'''
3630	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3631	if not len(mapData):
3632	return 0.0
3633	rho = copy.copy(mapData['rho']) #don't mess up original
3634	if not len(rho):
3635	return 0.0
3636	mapShape = np.array(rho.shape)
3637	mapStep = 1./mapShape
3638	X = np.array(xyz)%1. #get into unit cell
3639	I = np.array(X*mapShape,dtype='int')
3640	D = X-I*mapStep #position inside map cell
3641	D12 = D[0]*D[1]
3642	D13 = D[0]*D[2]
3643	D23 = D[1]*D[2]
3644	D123 = np.prod(D)
3645	Rho = rollMap(rho,-I) #shifts map so point is in corner
3646	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]+ \
3647	Rho[1,1,1]D123+Rho[0,1,1](D23-D123)+Rho[1,0,1](D13-D123)+Rho[1,1,0](D12-D123)+ \
3648	Rho[0,0,0](D12+D13+D23-D123)-Rho[0,0,1](D13+D23-D123)- \
3649	Rho[0,1,0](D23+D12-D123)-Rho[1,0,0](D13+D12-D123)
3650	return R
3651
3652	def getRhos(XYZ,rho):
3653	''' get scattering density at an array of point by 8-point interpolation
3654	this is faster than gerRho which is only used for single points. However, getRhos is
3655	replaced by scipy.ndimage.interpolation.map_coordinates which does a better job & is just as fast.
3656	Thus, getRhos is unused in GSAS-II at this time.
3657	param xyz: array coordinates to be probed Nx3
3658	param: rho: array copy of map (NB: don't use original!)
3659
3660	:returns: density at xyz
3661	'''
3662	def getBoxes(rho,I):
3663	Rhos = np.zeros((2,2,2))
3664	Mx,My,Mz = rho.shape
3665	Ix,Iy,Iz = I
3666	Rhos = np.array([[[rho[Ix%Mx,Iy%My,Iz%Mz],rho[Ix%Mx,Iy%My,(Iz+1)%Mz]],
3667	[rho[Ix%Mx,(Iy+1)%My,Iz%Mz],rho[Ix%Mx,(Iy+1)%My,(Iz+1)%Mz]]],
3668	[[rho[(Ix+1)%Mx,Iy%My,Iz%Mz],rho[(Ix+1)%Mx,Iy%My,(Iz+1)%Mz]],
3669	[rho[(Ix+1)%Mx,(Iy+1)%My,Iz%Mz],rho[(Ix+1)%Mx,(Iy+1)%My,(Iz+1)%Mz]]]])
3670	return Rhos
3671
3672	Blk = 400 #400 doesn't seem to matter
3673	nBlk = len(XYZ)//Blk #select Blk so this is an exact divide
3674	mapShape = np.array(rho.shape)
3675	mapStep = 1./mapShape
3676	X = XYZ%1. #get into unit cell
3677	iBeg = 0
3678	R = np.zeros(len(XYZ))
3679	#actually a lot faster!
3680	for iblk in range(nBlk):
3681	iFin = iBeg+Blk
3682	Xs = X[iBeg:iFin]
3683	I = np.array(np.rint(Xs*mapShape),dtype='int')
3684	Rhos = np.array([getBoxes(rho,i) for i in I])
3685	Ds = Xs-I*mapStep
3686	RIJs = Rhos[:,0,:2,:2]*(1.-Ds[:,0][:,nxs,nxs])
3687	RIs = RIJs[:,0](1.-Ds[:,1][:,nxs])+RIJs[:,1]Ds[:,1][:,nxs]
3688	R[iBeg:iFin] = RIs[:,0](1.-Ds[:,2])+RIs[:,1]Ds[:,2]
3689	iBeg += Blk
3690	return R
3691
3692	def SearchMap(generalData,drawingData,Neg=False):
3693	'''Does a search of a density map for peaks meeting the criterion of peak
3694	height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3695	mapData is data['General']['mapData']; the map is also in mapData.
3696
3697	:param generalData: the phase data structure; includes the map
3698	:param drawingData: the drawing data structure
3699	:param Neg: if True then search for negative peaks (i.e. H-atoms & neutron data)
3700
3701	:returns: (peaks,mags,dzeros) where
3702
3703	* peaks : ndarray
3704	x,y,z positions of the peaks found in the map
3705	* mags : ndarray
3706	the magnitudes of the peaks
3707	* dzeros : ndarray
3708	the distance of the peaks from the unit cell origin
3709	* dcent : ndarray
3710	the distance of the peaks from the unit cell center
3711
3712	'''
3713	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3714
3715	norm = 1./(np.sqrt(3.)np.sqrt(2.np.pi)**3)
3716
3717	def fixSpecialPos(xyz,SGData,Amat):
3718	equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3719	X = []
3720	xyzs = [equiv[0] for equiv in equivs]
3721	for x in xyzs:
3722	if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3723	X.append(x)
3724	if len(X) > 1:
3725	return np.average(X,axis=0)
3726	else:
3727	return xyz
3728
3729	def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3730	Mag,x0,y0,z0,sig = parms
3731	z = -((x0-rX)2+(y0-rY)2+(z0-rZ)*2)/(2.sig**2)
3732	# return normMagnp.exp(z)/(sigres*3) #not slower but some faults in LS
3733	return normMag(1.+z+z*2/2.)/(sigres**3)
3734
3735	def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3736	Mag,x0,y0,z0,sig = parms
3737	M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3738	return M
3739
3740	def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3741	Mag,x0,y0,z0,sig = parms
3742	dMdv = np.zeros(([5,]+list(rX.shape)))
3743	delt = .01
3744	for i in range(5):
3745	parms[i] -= delt
3746	rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3747	parms[i] += 2.*delt
3748	rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3749	parms[i] -= delt
3750	dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3751	rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3752	Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3753	dMdv = np.reshape(dMdv,(5,rX.size))
3754	Hess = np.inner(dMdv,dMdv)
3755
3756	return Vec,Hess
3757
3758	SGData = generalData['SGData']
3759	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3760	peaks = []
3761	mags = []
3762	dzeros = []
3763	dcent = []
3764	try:
3765	mapData = generalData['Map']
3766	contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3767	if Neg:
3768	rho = -copy.copy(mapData['rho']) #flip +/-
3769	else:
3770	rho = copy.copy(mapData['rho']) #don't mess up original
3771	mapHalf = np.array(rho.shape)/2
3772	res = mapData['GridStep']*2.
3773	incre = np.array(rho.shape,dtype=np.float)
3774	step = max(1.0,1./res)+1
3775	steps = np.array((3*[step,]),dtype='int32')
3776	except KeyError:
3777	G2fil.G2Print ('**** ERROR - Fourier map not defined')
3778	return peaks,mags
3779	rhoMask = ma.array(rho,mask=(rho<contLevel))
3780	indices = (-1,0,1)
3781	rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3782	for roll in rolls:
3783	if np.any(roll):
3784	rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3785	indx = np.transpose(rhoMask.nonzero())
3786	peaks = indx/incre
3787	mags = rhoMask[rhoMask.nonzero()]
3788	for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3789	rho = rollMap(rho,ind)
3790	rMM = mapHalf-steps
3791	rMP = mapHalf+steps+1
3792	rhoPeak = rho[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3793	peakInt = np.sum(rhoPeak)res*3
3794	rX,rY,rZ = np.mgrid[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3795	x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0] #magnitude, position & width(sig)
3796	result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3797	args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3798	x1 = result[0]
3799	if not np.any(x1 < 0):
3800	peak = (np.array(x1[1:4])-ind)/incre
3801	peak = fixSpecialPos(peak,SGData,Amat)
3802	rho = rollMap(rho,-ind)
3803	cent = np.ones(3)*.5
3804	dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3805	dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3806	if Neg: #want negative magnitudes for negative peaks
3807	return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3808	else:
3809	return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3810
3811	def sortArray(data,pos,reverse=False):
3812	'''data is a list of items
3813	sort by pos in list; reverse if True
3814	'''
3815	T = []
3816	for i,M in enumerate(data):
3817	try:
3818	T.append((M[pos],i))
3819	except IndexError:
3820	return data
3821	D = dict(zip(T,data))
3822	T.sort()
3823	if reverse:
3824	T.reverse()
3825	X = []
3826	for key in T:
3827	X.append(D[key])
3828	return X
3829
3830	def PeaksEquiv(data,Ind):
3831	'''Find the equivalent map peaks for those selected. Works on the
3832	contents of data['Map Peaks'].
3833
3834	:param data: the phase data structure
3835	:param list Ind: list of selected peak indices
3836	:returns: augmented list of peaks including those related by symmetry to the
3837	ones in Ind
3838
3839	'''
3840	def Duplicate(xyz,peaks,Amat):
3841	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3842	return True
3843	return False
3844
3845	generalData = data['General']
3846	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3847	SGData = generalData['SGData']
3848	mapPeaks = data['Map Peaks']
3849	XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3850	Indx = {}
3851	for ind in Ind:
3852	xyz = np.array(mapPeaks[ind][1:4])
3853	xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3854	for jnd,xyz in enumerate(XYZ):
3855	Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3856	Ind = []
3857	for ind in Indx:
3858	if Indx[ind]:
3859	Ind.append(ind)
3860	return Ind
3861
3862	def PeaksUnique(data,Ind):
3863	'''Finds the symmetry unique set of peaks from those selected. Works on the
3864	contents of data['Map Peaks'].
3865
3866	:param data: the phase data structure
3867	:param list Ind: list of selected peak indices
3868	:returns: the list of symmetry unique peaks from among those given in Ind
3869
3870	'''
3871	# XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3872
3873	def noDuplicate(xyz,peaks,Amat):
3874	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3875	return False
3876	return True
3877
3878	generalData = data['General']
3879	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3880	SGData = generalData['SGData']
3881	mapPeaks = data['Map Peaks']
3882	Indx = {}
3883	XYZ = {}
3884	for ind in Ind:
3885	XYZ[ind] = np.array(mapPeaks[ind][1:4])
3886	Indx[ind] = True
3887	for ind in Ind:
3888	if Indx[ind]:
3889	xyz = XYZ[ind]
3890	for jnd in Ind:
3891	if ind != jnd and Indx[jnd]:
3892	Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3893	xyzs = np.array([equiv[0] for equiv in Equiv])
3894	Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3895	Ind = []
3896	for ind in Indx:
3897	if Indx[ind]:
3898	Ind.append(ind)
3899	return Ind
3900
3901	################################################################################
3902	##### Dysnomia setup & return stuff
3903	################################################################################
3904
3905
3906
3907	################################################################################
3908	##### single peak fitting profile fxn stuff
3909	################################################################################
3910
3911	def getCWsig(ins,pos):
3912	'''get CW peak profile sigma^2
3913
3914	:param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3915	as values only
3916	:param float pos: 2-theta of peak
3917	:returns: float getCWsig: peak sigma^2
3918
3919	'''
3920	tp = tand(pos/2.0)
3921	return ins['U']tp2+ins['V']tp+ins['W']
3922
3923	def getCWsigDeriv(pos):
3924	'''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3925
3926	:param float pos: 2-theta of peak
3927
3928	:returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3929
3930	'''
3931	tp = tand(pos/2.0)
3932	return tp**2,tp,1.0
3933
3934	def getCWgam(ins,pos):
3935	'''get CW peak profile gamma
3936
3937	:param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3938	as values only
3939	:param float pos: 2-theta of peak
3940	:returns: float getCWgam: peak gamma
3941
3942	'''
3943	return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)+ins.get('Z',0.0)
3944
3945	def getCWgamDeriv(pos):
3946	'''get derivatives of CW peak profile gamma wrt X, Y & Z
3947
3948	:param float pos: 2-theta of peak
3949
3950	:returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3951
3952	'''
3953	return 1./cosd(pos/2.0),tand(pos/2.0),1.0
3954
3955	def getTOFsig(ins,dsp):
3956	'''get TOF peak profile sigma^2
3957
3958	:param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3959	as values only
3960	:param float dsp: d-spacing of peak
3961
3962	:returns: float getTOFsig: peak sigma^2
3963
3964	'''
3965	return ins['sig-0']+ins['sig-1']dsp2+ins['sig-2']dsp*4+ins['sig-q']dsp
3966
3967	def getTOFsigDeriv(dsp):
3968	'''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3969
3970	:param float dsp: d-spacing of peak
3971
3972	:returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3973
3974	'''
3975	return 1.0,dsp2,dsp4,dsp
3976
3977	def getTOFgamma(ins,dsp):
3978	'''get TOF peak profile gamma
3979
3980	:param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3981	as values only
3982	:param float dsp: d-spacing of peak
3983
3984	:returns: float getTOFgamma: peak gamma
3985
3986	'''
3987	return ins['Z']+ins['X']dsp+ins['Y']dsp**2
3988
3989	def getTOFgammaDeriv(dsp):
3990	'''get derivatives of TOF peak profile gamma wrt X, Y & Z
3991
3992	:param float dsp: d-spacing of peak
3993
3994	:returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3995
3996	'''
3997	return dsp,dsp**2,1.0
3998
3999	def getTOFbeta(ins,dsp):
4000	'''get TOF peak profile beta
4001
4002	:param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
4003	as values only
4004	:param float dsp: d-spacing of peak
4005
4006	:returns: float getTOFbeta: peak beat
4007
4008	'''
4009	return ins['beta-0']+ins['beta-1']/dsp4+ins['beta-q']/dsp2
4010
4011	def getTOFbetaDeriv(dsp):
4012	'''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
4013
4014	:param float dsp: d-spacing of peak
4015
4016	:returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
4017
4018	'''
4019	return 1.0,1./dsp4,1./dsp2
4020
4021	def getTOFalpha(ins,dsp):
4022	'''get TOF peak profile alpha
4023
4024	:param dict ins: instrument parameters with at least 'alpha'
4025	as values only
4026	:param float dsp: d-spacing of peak
4027
4028	:returns: flaot getTOFalpha: peak alpha
4029
4030	'''
4031	return ins['alpha']/dsp
4032
4033	def getTOFalphaDeriv(dsp):
4034	'''get derivatives of TOF peak profile beta wrt alpha
4035
4036	:param float dsp: d-spacing of peak
4037
4038	:returns: float getTOFalphaDeriv: d(alp)/d(alpha)
4039
4040	'''
4041	return 1./dsp
4042
4043	def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
4044	'''set starting peak parameters for single peak fits from plot selection or auto selection
4045
4046	:param dict Parms: instrument parameters dictionary
4047	:param dict Parms2: table lookup for TOF profile coefficients
4048	:param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
4049	:param float mag: peak top magnitude from pick
4050	:param bool ifQ: True if pos in Q
4051	:param bool useFit: True if use fitted CW Parms values (not defaults)
4052
4053	:returns: list XY: peak list entry:
4054	for CW: [pos,0,mag,1,sig,0,gam,0]
4055	for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4056	NB: mag refinement set by default, all others off
4057
4058	'''
4059	ind = 0
4060	if useFit:
4061	ind = 1
4062	ins = {}
4063	if 'C' in Parms['Type'][0]: #CW data - TOF later in an else
4064	for x in ['U','V','W','X','Y','Z']:
4065	ins[x] = Parms.get(x,[0.0,0.0])[ind]
4066	if ifQ: #qplot - convert back to 2-theta
4067	pos = 2.0asind(posgetWave(Parms)/(4*math.pi))
4068	sig = getCWsig(ins,pos)
4069	gam = getCWgam(ins,pos)
4070	XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st
4071	else:
4072	if ifQ:
4073	dsp = 2.*np.pi/pos
4074	pos = Parms['difC']*dsp
4075	else:
4076	dsp = pos/Parms['difC'][1]
4077	if 'Pdabc' in Parms2:
4078	for x in ['sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4079	ins[x] = Parms.get(x,[0.0,0.0])[ind]
4080	Pdabc = Parms2['Pdabc'].T
4081	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
4082	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
4083	else:
4084	for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4085	ins[x] = Parms.get(x,[0.0,0.0])[ind]
4086	alp = getTOFalpha(ins,dsp)
4087	bet = getTOFbeta(ins,dsp)
4088	sig = getTOFsig(ins,dsp)
4089	gam = getTOFgamma(ins,dsp)
4090	XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4091	return XY
4092
4093	################################################################################
4094	##### MC/SA stuff
4095	################################################################################
4096
4097	#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
4098	# Original Author: Travis Oliphant 2002
4099	# Bug-fixes in 2006 by Tim Leslie
4100
4101
4102	import numpy
4103	from numpy import asarray, exp, squeeze, sign, \
4104	all, shape, array, where
4105	from numpy import random
4106
4107	#__all__ = ['anneal']
4108
4109	_double_min = numpy.finfo(float).min
4110	_double_max = numpy.finfo(float).max
4111	class base_schedule(object):
4112	def __init__(self):
4113	self.dwell = 20
4114	self.lower = 0.
4115	self.upper = 1.
4116	self.Ninit = 50
4117	self.accepted = 0
4118	self.tests = 0
4119	self.feval = 0
4120	self.k = 0
4121	self.T = None
4122
4123	def init(self, **options):
4124	self.__dict__.update(options)
4125	self.lower = asarray(self.lower)
4126	self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
4127	self.upper = asarray(self.upper)
4128	self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
4129	self.k = 0
4130	self.accepted = 0
4131	self.feval = 0
4132	self.tests = 0
4133
4134	def getstart_temp(self, best_state):
4135	""" Find a matching starting temperature and starting parameters vector
4136	i.e. find x0 such that func(x0) = T0.
4137
4138	:param best_state: _state
4139	A _state object to store the function value and x0 found.
4140
4141	:returns: x0 : array
4142	The starting parameters vector.
4143	"""
4144
4145	assert(not self.dims is None)
4146	lrange = self.lower
4147	urange = self.upper
4148	fmax = _double_min
4149	fmin = _double_max
4150	for _ in range(self.Ninit):
4151	x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
4152	fval = self.func(x0, *self.args)
4153	self.feval += 1
4154	if fval > fmax:
4155	fmax = fval
4156	if fval < fmin:
4157	fmin = fval
4158	best_state.cost = fval
4159	best_state.x = array(x0)
4160
4161	self.T0 = (fmax-fmin)*1.5
4162	return best_state.x
4163
4164	def set_range(self,x0,frac):
4165	delrange = frac*(self.upper-self.lower)
4166	self.upper = x0+delrange
4167	self.lower = x0-delrange
4168
4169	def accept_test(self, dE):
4170	T = self.T
4171	self.tests += 1
4172	if dE < 0:
4173	self.accepted += 1
4174	return 1
4175	p = exp(-dE*1.0/T)
4176	if (p > random.uniform(0.0, 1.0)):
4177	self.accepted += 1
4178	return 1
4179	return 0
4180
4181	def update_guess(self, x0):
4182	return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
4183
4184	def update_temp(self, x0):
4185	pass
4186
4187	class fast_sa(base_schedule):
4188	def init(self, **options):
4189	self.__dict__.update(options)
4190
4191	def update_guess(self, x0):
4192	x0 = asarray(x0)
4193	u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4194	T = self.T
4195	xc = (sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)+1.0)/2.0
4196	xnew = xc*(self.upper - self.lower)+self.lower
4197	return xnew
4198
4199	def update_temp(self):
4200	self.T = self.T0exp(-self.c self.k**(self.quench))
4201	self.k += 1
4202	return
4203
4204	class log_sa(base_schedule): #OK
4205
4206	def init(self,**options):
4207	self.__dict__.update(options)
4208
4209	def update_guess(self,x0): #same as default #TODO - is this a reasonable update procedure?
4210	u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4211	T = self.T
4212	xc = (sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)+1.0)/2.0
4213	xnew = xc*(self.upper - self.lower)+self.lower
4214	return xnew
4215
4216	def update_temp(self):
4217	self.k += 1
4218	self.T = self.T0self.slope*self.k
4219
4220	class _state(object):
4221	def __init__(self):
4222	self.x = None
4223	self.cost = None
4224
4225	def makeTsched(data):
4226	if data['Algorithm'] == 'fast':
4227	sched = fast_sa()
4228	sched.quench = data['fast parms'][0]
4229	sched.c = data['fast parms'][1]
4230	elif data['Algorithm'] == 'log':
4231	sched = log_sa()
4232	sched.slope = data['log slope']
4233	sched.T0 = data['Annealing'][0]
4234	if not sched.T0:
4235	sched.T0 = 50.
4236	Tf = data['Annealing'][1]
4237	if not Tf:
4238	Tf = 0.001
4239	Tsched = [sched.T0,]
4240	while Tsched[-1] > Tf:
4241	sched.update_temp()
4242	Tsched.append(sched.T)
4243	return Tsched[1:]
4244
4245	def anneal(func, x0, args=(), schedule='fast',
4246	T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4247	feps=1e-6, quench=1.0, c=1.0,
4248	lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4249	ranRange=0.10,autoRan=False,dlg=None):
4250	"""Minimize a function using simulated annealing.
4251
4252	Schedule is a schedule class implementing the annealing schedule.
4253	Available ones are 'fast', 'cauchy', 'boltzmann'
4254
4255	:param callable func: f(x, \*args)
4256	Function to be optimized.
4257	:param ndarray x0:
4258	Initial guess.
4259	:param tuple args:
4260	Extra parameters to `func`.
4261	:param base_schedule schedule:
4262	Annealing schedule to use (a class).
4263	:param float T0:
4264	Initial Temperature (estimated as 1.2 times the largest
4265	cost-function deviation over random points in the range).
4266	:param float Tf:
4267	Final goal temperature.
4268	:param int maxeval:
4269	Maximum function evaluations.
4270	:param int maxaccept:
4271	Maximum changes to accept.
4272	:param int maxiter:
4273	Maximum cooling iterations.
4274	:param float feps:
4275	Stopping relative error tolerance for the function value in
4276	last four coolings.
4277	:param float quench,c:
4278	Parameters to alter fast_sa schedule.
4279	:param float/ndarray lower,upper:
4280	Lower and upper bounds on `x`.
4281	:param int dwell:
4282	The number of times to search the space at each temperature.
4283	:param float slope:
4284	Parameter for log schedule
4285	:param bool ranStart=False:
4286	True for set 10% of ranges about x
4287
4288	:returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4289
4290	* xmin (ndarray): Point giving smallest value found.
4291	* Jmin (float): Minimum value of function found.
4292	* T (float): Final temperature.
4293	* feval (int): Number of function evaluations.
4294	* iters (int): Number of cooling iterations.
4295	* accept (int): Number of tests accepted.
4296	* retval (int): Flag indicating stopping condition:
4297
4298	* 0: Points no longer changing
4299	* 1: Cooled to final temperature
4300	* 2: Maximum function evaluations
4301	* 3: Maximum cooling iterations reached
4302	* 4: Maximum accepted query locations reached
4303	* 5: Final point not the minimum amongst encountered points
4304
4305	Notes:
4306	Simulated annealing is a random algorithm which uses no derivative
4307	information from the function being optimized. In practice it has
4308	been more useful in discrete optimization than continuous
4309	optimization, as there are usually better algorithms for continuous
4310	optimization problems.
4311
4312	Some experimentation by trying the difference temperature
4313	schedules and altering their parameters is likely required to
4314	obtain good performance.
4315
4316	The randomness in the algorithm comes from random sampling in numpy.
4317	To obtain the same results you can call numpy.random.seed with the
4318	same seed immediately before calling scipy.optimize.anneal.
4319
4320	We give a brief description of how the three temperature schedules
4321	generate new points and vary their temperature. Temperatures are
4322	only updated with iterations in the outer loop. The inner loop is
4323	over range(dwell), and new points are generated for
4324	every iteration in the inner loop. (Though whether the proposed
4325	new points are accepted is probabilistic.)
4326
4327	For readability, let d denote the dimension of the inputs to func.
4328	Also, let x_old denote the previous state, and k denote the
4329	iteration number of the outer loop. All other variables not
4330	defined below are input variables to scipy.optimize.anneal itself.
4331
4332	In the 'fast' schedule the updates are ::
4333
4334	u ~ Uniform(0, 1, size=d)
4335	y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4336	xc = y * (upper - lower)
4337	x_new = x_old + xc
4338
4339	T_new = T0 * exp(-c * k**quench)
4340
4341	"""
4342
4343	''' Scipy license:
4344	Copyright (c) 2001, 2002 Enthought, Inc.
4345	All rights reserved.
4346
4347	Copyright (c) 2003-2016 SciPy Developers.
4348	All rights reserved.
4349
4350	Redistribution and use in source and binary forms, with or without
4351	modification, are permitted provided that the following conditions are met:
4352
4353	a. Redistributions of source code must retain the above copyright notice,
4354	this list of conditions and the following disclaimer.
4355	b. Redistributions in binary form must reproduce the above copyright
4356	notice, this list of conditions and the following disclaimer in the
4357	documentation and/or other materials provided with the distribution.
4358	c. Neither the name of Enthought nor the names of the SciPy Developers
4359	may be used to endorse or promote products derived from this software
4360	without specific prior written permission.
4361	'''
4362	x0 = asarray(x0)
4363	lower = asarray(lower)
4364	upper = asarray(upper)
4365
4366	schedule = eval(schedule+'_sa()')
4367	# initialize the schedule
4368	schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4369	c=c, quench=quench, dwell=dwell, slope=slope)
4370
4371	current_state, last_state, best_state = _state(), _state(), _state()
4372	if ranStart:
4373	schedule.set_range(x0,ranRange)
4374	if T0 is None:
4375	x0 = schedule.getstart_temp(best_state)
4376	else:
4377	x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4378	best_state.x = None
4379	best_state.cost = numpy.Inf
4380
4381	last_state.x = asarray(x0).copy()
4382	fval = func(x0,*args)
4383	schedule.feval += 1
4384	last_state.cost = fval
4385	if last_state.cost < best_state.cost:
4386	best_state.cost = fval
4387	best_state.x = asarray(x0).copy()
4388	schedule.T = schedule.T0
4389	fqueue = [100, 300, 500, 700]
4390	iters = 1
4391	keepGoing = True
4392	bestn = 0
4393	while keepGoing:
4394	retval = 0
4395	for n in range(dwell):
4396	current_state.x = schedule.update_guess(last_state.x)
4397	current_state.cost = func(current_state.x,*args)
4398	schedule.feval += 1
4399
4400	dE = current_state.cost - last_state.cost
4401	if schedule.accept_test(dE):
4402	last_state.x = current_state.x.copy()
4403	last_state.cost = current_state.cost
4404	if last_state.cost < best_state.cost:
4405	best_state.x = last_state.x.copy()
4406	best_state.cost = last_state.cost
4407	bestn = n
4408	if best_state.cost < 1.0 and autoRan:
4409	schedule.set_range(x0,best_state.cost/2.)
4410	if dlg:
4411	GoOn = dlg.Update(min(100.,best_state.cost*100),
4412	newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4413	'Best trial:',bestn, \
4414	'MC/SA Residual:',best_state.cost*100,'%', \
4415	))[0]
4416	if not GoOn:
4417	best_state.x = last_state.x.copy()
4418	best_state.cost = last_state.cost
4419	retval = 5
4420	schedule.update_temp()
4421	iters += 1
4422	# Stopping conditions
4423	# 0) last saved values of f from each cooling step
4424	# are all very similar (effectively cooled)
4425	# 1) Tf is set and we are below it
4426	# 2) maxeval is set and we are past it
4427	# 3) maxiter is set and we are past it
4428	# 4) maxaccept is set and we are past it
4429	# 5) user canceled run via progress bar
4430
4431	fqueue.append(squeeze(last_state.cost))
4432	fqueue.pop(0)
4433	af = asarray(fqueue)*1.0
4434	if retval == 5:
4435	G2fil.G2Print ('Error: User terminated run; incomplete MC/SA')
4436	keepGoing = False
4437	break
4438	if all(abs((af-af[0])/af[0]) < feps):
4439	retval = 0
4440	if abs(af[-1]-best_state.cost) > feps*10:
4441	retval = 5
4442	G2fil.G2Print (" Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1))
4443	break
4444	if (Tf is not None) and (schedule.T < Tf):
4445	# print ' Minimum T reached in %d steps'%(iters-1)
4446	retval = 1
4447	break
4448	if (maxeval is not None) and (schedule.feval > maxeval):
4449	retval = 2
4450	break
4451	if (iters > maxiter):
4452	G2fil.G2Print (" Warning: Maximum number of iterations exceeded.")
4453	retval = 3
4454	break
4455	if (maxaccept is not None) and (schedule.accepted > maxaccept):
4456	retval = 4
4457	break
4458
4459	return best_state.x, best_state.cost, schedule.T, \
4460	schedule.feval, iters, schedule.accepted, retval
4461
4462	def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4463	outlist = []
4464	timelist = []
4465	nsflist = []
4466	random.seed(int(time.time())%100000+nprocess) #make sure each process has a different random start
4467	for n in range(iCyc):
4468	result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False) #mcsa result,time,rcov
4469	outlist.append(result[0])
4470	timelist.append(result[1])
4471	nsflist.append(result[2])
4472	G2fil.G2Print (' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1]))
4473	out_q.put(outlist)
4474	out_t.put(timelist)
4475	out_n.put(nsflist)
4476
4477	def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4478	import multiprocessing as mp
4479
4480	out_q = mp.Queue()
4481	out_t = mp.Queue()
4482	out_n = mp.Queue()
4483	procs = []
4484	totsftime = 0.
4485	totnsf = 0
4486	iCyc = np.zeros(nprocs)
4487	for i in range(nCyc):
4488	iCyc[i%nprocs] += 1
4489	for i in range(nprocs):
4490	p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4491	procs.append(p)
4492	p.start()
4493	resultlist = []
4494	for i in range(nprocs):
4495	resultlist += out_q.get()
4496	totsftime += np.sum(out_t.get())
4497	totnsf += np.sum(out_n.get())
4498	for p in procs:
4499	p.join()
4500	return resultlist,totsftime,totnsf
4501
4502	def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4503	'''default doc string
4504
4505	:param type name: description
4506
4507	:returns: type name: description
4508	'''
4509
4510	class RandomDisplacementBounds(object):
4511	"""random displacement with bounds"""
4512	def __init__(self, xmin, xmax, stepsize=0.5):
4513	self.xmin = xmin
4514	self.xmax = xmax
4515	self.stepsize = stepsize
4516
4517	def __call__(self, x):
4518	"""take a random step but ensure the new position is within the bounds"""
4519	while True:
4520	# this could be done in a much more clever way, but it will work for example purposes
4521	steps = self.xmax-self.xmin
4522	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
4523	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4524	break
4525	return xnew
4526
4527	global tsum,nsum
4528	tsum = 0.
4529	nsum = 0
4530
4531	def getMDparms(item,pfx,parmDict,varyList):
4532	parmDict[pfx+'MDaxis'] = item['axis']
4533	parmDict[pfx+'MDval'] = item['Coef'][0]
4534	if item['Coef'][1]:
4535	varyList += [pfx+'MDval',]
4536	limits = item['Coef'][2]
4537	lower.append(limits[0])
4538	upper.append(limits[1])
4539
4540	def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4541	parmDict[pfx+'Atype'] = item['atType']
4542	aTypes \|= set([item['atType'],])
4543	pstr = ['Ax','Ay','Az']
4544	XYZ = [0,0,0]
4545	for i in range(3):
4546	name = pfx+pstr[i]
4547	parmDict[name] = item['Pos'][0][i]
4548	XYZ[i] = parmDict[name]
4549	if item['Pos'][1][i]:
4550	varyList += [name,]
4551	limits = item['Pos'][2][i]
4552	lower.append(limits[0])
4553	upper.append(limits[1])
4554	parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4555
4556	def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4557	parmDict[mfx+'MolCent'] = item['MolCent']
4558	parmDict[mfx+'RBId'] = item['RBId']
4559	pstr = ['Px','Py','Pz']
4560	ostr = ['Qa','Qi','Qj','Qk'] #angle,vector not quaternion
4561	for i in range(3):
4562	name = mfx+pstr[i]
4563	parmDict[name] = item['Pos'][0][i]
4564	if item['Pos'][1][i]:
4565	varyList += [name,]
4566	limits = item['Pos'][2][i]
4567	lower.append(limits[0])
4568	upper.append(limits[1])
4569	AV = item['Ori'][0]
4570	A = AV[0]
4571	V = AV[1:]
4572	for i in range(4):
4573	name = mfx+ostr[i]
4574	if i:
4575	parmDict[name] = V[i-1]