Context Navigation

source: trunk/GSASIImath.py @ 3893

Last change on this file since 3893 was 3893, checked in by vondreele, 4 years ago
remove commented lines & "improve" getRhos - more numpy oriented - not faster than original despite elegance
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 189.3 KB

Line
1	# -- coding: utf-8 --
2	#GSASIImath - major mathematics routines
3	########### SVN repository information ###################
4	# $Date: 2019-04-12 20:35:08 +0000 (Fri, 12 Apr 2019) $
5	# $Author: vondreele $
6	# $Revision: 3893 $
7	# $URL: trunk/GSASIImath.py $
8	# $Id: GSASIImath.py 3893 2019-04-12 20:35:08Z 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: 3893 $")
27	import GSASIIElem as G2el
28	import GSASIIlattice as G2lat
29	import GSASIIspc as G2spc
30	import GSASIIpwd as G2pwd
31	import numpy.fft as fft
32	import scipy.optimize as so
33	try:
34	import pypowder as pwd
35	except ImportError:
36	print ('pypowder is not available - profile calcs. not allowed')
37
38	sind = lambda x: np.sin(x*np.pi/180.)
39	cosd = lambda x: np.cos(x*np.pi/180.)
40	tand = lambda x: np.tan(x*np.pi/180.)
41	asind = lambda x: 180.*np.arcsin(x)/np.pi
42	acosd = lambda x: 180.*np.arccos(x)/np.pi
43	atand = lambda x: 180.*np.arctan(x)/np.pi
44	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
45	twopi = 2.0*np.pi
46	twopisq = 2.0np.pi*2
47	nxs = np.newaxis
48
49	################################################################################
50	##### Hessian least-squares Levenberg-Marquardt routine
51	################################################################################
52
53	def pinv(a, rcond=1e-15 ):
54	"""
55	Compute the (Moore-Penrose) pseudo-inverse of a matrix.
56	Modified from numpy.linalg.pinv; assumes a is Hessian & returns no. zeros found
57	Calculate the generalized inverse of a matrix using its
58	singular-value decomposition (SVD) and including all
59	large singular values.
60
61	:param array a: (M, M) array_like - here assumed to be LS Hessian
62	Matrix to be pseudo-inverted.
63	:param float rcond: Cutoff for small singular values.
64	Singular values smaller (in modulus) than
65	`rcond` * largest_singular_value (again, in modulus)
66	are set to zero.
67
68	:returns: B : (M, M) ndarray
69	The pseudo-inverse of `a`
70
71	Raises: LinAlgError
72	If the SVD computation does not converge.
73
74	Notes:
75	The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
76	defined as: "the matrix that 'solves' [the least-squares problem]
77	:math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
78	:math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
79
80	It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
81	value decomposition of A, then
82	:math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
83	orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
84	of A's so-called singular values, (followed, typically, by
85	zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
86	consisting of the reciprocals of A's singular values
87	(again, followed by zeros). [1]
88
89	References:
90	.. [1] G. Strang, Linear Algebra and Its Applications, 2nd Ed., Orlando, FL, Academic Press, Inc., 1980, pp. 139-142.
91
92	"""
93	u, s, vt = nl.svd(a)
94	cutoff = rcond*np.maximum.reduce(s)
95	s = np.where(s>cutoff,1./s,0.)
96	nzero = s.shape[0]-np.count_nonzero(s)
97	# res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))
98	res = np.dot(vt.T,s[:,nxs]*u.T)
99	return res,nzero
100
101	def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
102
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	print (' 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	print ('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	print ('ouch #1 bad SVD inversion; change parameterization for ',psing)
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	print ('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	print ('ouch #3 chisq1 %g.4 stuck > chisq0 %g.4'%(chisq1,chisq0))
204	break
205	chitol *= 2
206	lamMax = max(lamMax,lam)
207	deltaChi2 = (chisq0-chisq1)/chisq0
208	if Print:
209	print (' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g'%(
210	icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
211	if deltaChi2 < ftol:
212	ifConverged = True
213	if Print: print ("converged")
214	break
215	icycle += 1
216	M = func(x0,*args)
217	nfev += 1
218	Yvec,Amat = Hess(x0,*args)
219	Adiag = np.sqrt(np.diag(Amat))
220	Anorm = np.outer(Adiag,Adiag)
221	Lam = np.eye(Amat.shape[0])*lam
222	Amatlam = Amat/Anorm
223	try:
224	Bmat,Nzero = pinv(Amatlam,xtol) #Moore-Penrose inversion (via SVD) & count of zeros
225	if Print: print ('Found %d SVD zeros'%(Nzero))
226	Bmat = Bmat/Anorm
227	return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
228	except nl.LinAlgError:
229	print ('ouch #2 linear algebra error in making v-cov matrix')
230	psing = []
231	if maxcyc:
232	psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
233	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]
234
235	def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
236
237	"""
238	Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
239	values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
240	where :math:`x = arg min(\sum_{y=0}^{N_{obs}} (func(y)^2,axis=0))`
241
242	:param function func: callable method or function
243	should take at least one (possibly length N vector) argument and
244	returns M floating point numbers.
245	:param np.ndarray x0: The starting estimate for the minimization of length N
246	:param function Hess: callable method or function
247	A required function or method to compute the weighted vector and Hessian for func.
248	It must be a symmetric NxN array
249	:param tuple args: Any extra arguments to func are placed in this tuple.
250	:param float ftol: Relative error desired in the sum of squares.
251	:param float xtol: Relative tolerance of zeros in the SVD solution in nl.pinv.
252	:param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
253	until other limits are met (ftol, xtol)
254	:param bool Print: True for printing results (residuals & times) by cycle
255
256	:returns: (x,cov_x,infodict) where
257
258	* x : ndarray
259	The solution (or the result of the last iteration for an unsuccessful
260	call).
261	* cov_x : ndarray
262	Uses the fjac and ipvt optional outputs to construct an
263	estimate of the jacobian around the solution. ``None`` if a
264	singular matrix encountered (indicates very flat curvature in
265	some direction). This matrix must be multiplied by the
266	residual standard deviation to get the covariance of the
267	parameter estimates -- see curve_fit.
268	* infodict : dict
269	a dictionary of optional outputs with the keys:
270
271	* 'fvec' : the function evaluated at the output
272	* 'num cyc':
273	* 'nfev':
274	* 'lamMax':0.
275	* 'psing':
276	* 'SVD0':
277
278	"""
279
280	ifConverged = False
281	deltaChi2 = -10.
282	x0 = np.array(x0, ndmin=1) #might be redundant?
283	n = len(x0)
284	if type(args) != type(()):
285	args = (args,)
286
287	icycle = 0
288	nfev = 0
289	if Print:
290	print (' Hessian SVD refinement on %d variables:'%(n))
291	while icycle < maxcyc:
292	time0 = time.time()
293	M = func(x0,*args)
294	nfev += 1
295	chisq0 = np.sum(M**2)
296	Yvec,Amat = Hess(x0,*args)
297	Adiag = np.sqrt(np.diag(Amat))
298	psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
299	if np.any(psing): #hard singularity in matrix
300	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
301	Anorm = np.outer(Adiag,Adiag)
302	Yvec /= Adiag
303	Amat /= Anorm
304	if Print:
305	print ('initial chi^2 %.5g'%(chisq0))
306	try:
307	Ainv,Nzeros = pinv(Amat,xtol) #do Moore-Penrose inversion (via SVD)
308	except nl.LinAlgError:
309	print ('ouch #1 bad SVD inversion; change parameterization')
310	psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
311	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
312	Xvec = np.inner(Ainv,Yvec) #solve
313	Xvec /= Adiag
314	M2 = func(x0+Xvec,*args)
315	nfev += 1
316	chisq1 = np.sum(M2**2)
317	deltaChi2 = (chisq0-chisq1)/chisq0
318	if Print:
319	print (' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Delta: %.3g'%(
320	icycle,time.time()-time0,chisq1,deltaChi2))
321	if deltaChi2 < ftol:
322	ifConverged = True
323	if Print: print ("converged")
324	break
325	icycle += 1
326	M = func(x0,*args)
327	nfev += 1
328	Yvec,Amat = Hess(x0,*args)
329	Adiag = np.sqrt(np.diag(Amat))
330	Anorm = np.outer(Adiag,Adiag)
331	Amat = Amat/Anorm
332	try:
333	Bmat,Nzero = pinv(Amat,xtol) #Moore-Penrose inversion (via SVD) & count of zeros
334	print ('Found %d SVD zeros'%(Nzero))
335	# Bmat = nl.inv(Amatlam); Nzeros = 0
336	Bmat = Bmat/Anorm
337	return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':[],
338	'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
339	except nl.LinAlgError:
340	print ('ouch #2 linear algebra error in making v-cov matrix')
341	psing = []
342	if maxcyc:
343	psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
344	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]
345
346	def getVCov(varyNames,varyList,covMatrix):
347	'''obtain variance-covariance terms for a set of variables. NB: the varyList
348	and covMatrix were saved by the last least squares refinement so they must match.
349
350	:param list varyNames: variable names to find v-cov matric for
351	:param list varyList: full list of all variables in v-cov matrix
352	:param nparray covMatrix: full variance-covariance matrix from the last
353	least squares refinement
354
355	:returns: nparray vcov: variance-covariance matrix for the variables given
356	in varyNames
357
358	'''
359	vcov = np.zeros((len(varyNames),len(varyNames)))
360	for i1,name1 in enumerate(varyNames):
361	for i2,name2 in enumerate(varyNames):
362	try:
363	vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
364	except ValueError:
365	vcov[i1][i2] = 0.0
366	# if i1 == i2:
367	# vcov[i1][i2] = 1e-20
368	# else:
369	# vcov[i1][i2] = 0.0
370	return vcov
371
372	################################################################################
373	##### Atom manipulations
374	################################################################################
375
376	def FindMolecule(ind,generalData,atomData): #uses numpy & masks - very fast even for proteins!
377
378	def getNeighbors(atom,radius):
379	Dx = UAtoms-np.array(atom[cx:cx+3])
380	dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of disorder "bonds" < 0.5A
381	sumR = Radii+radius
382	return set(ma.nonzero(ma.masked_greater(dist-factor*sumR,0.))[0]) #get indices of bonded atoms
383
384	import numpy.ma as ma
385	indices = (-1,0,1)
386	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices],dtype='f')
387	cx,ct,cs,ci = generalData['AtomPtrs']
388	DisAglCtls = generalData['DisAglCtls']
389	SGData = generalData['SGData']
390	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
391	radii = DisAglCtls['BondRadii']
392	atomTypes = DisAglCtls['AtomTypes']
393	factor = DisAglCtls['Factors'][0]
394	unit = np.zeros(3)
395	try:
396	indH = atomTypes.index('H')
397	radii[indH] = 0.5
398	except:
399	pass
400	nAtom = len(atomData)
401	Indx = range(nAtom)
402	UAtoms = []
403	Radii = []
404	for atom in atomData:
405	UAtoms.append(np.array(atom[cx:cx+3]))
406	Radii.append(radii[atomTypes.index(atom[ct])])
407	UAtoms = np.array(UAtoms)
408	Radii = np.array(Radii)
409	for nOp,Op in enumerate(SGData['SGOps'][1:]):
410	UAtoms = np.concatenate((UAtoms,(np.inner(Op[0],UAtoms[:nAtom]).T+Op[1])))
411	Radii = np.concatenate((Radii,Radii[:nAtom]))
412	Indx += Indx[:nAtom]
413	for icen,cen in enumerate(SGData['SGCen'][1:]):
414	UAtoms = np.concatenate((UAtoms,(UAtoms+cen)))
415	Radii = np.concatenate((Radii,Radii))
416	Indx += Indx[:nAtom]
417	if SGData['SGInv']:
418	UAtoms = np.concatenate((UAtoms,-UAtoms))
419	Radii = np.concatenate((Radii,Radii))
420	Indx += Indx
421	UAtoms %= 1.
422	mAtoms = len(UAtoms)
423	for unit in Units:
424	if np.any(unit): #skip origin cell
425	UAtoms = np.concatenate((UAtoms,UAtoms[:mAtoms]+unit))
426	Radii = np.concatenate((Radii,Radii[:mAtoms]))
427	Indx += Indx[:mAtoms]
428	UAtoms = np.array(UAtoms)
429	Radii = np.array(Radii)
430	newAtoms = [atomData[ind],]
431	atomData[ind] = None
432	radius = Radii[ind]
433	IndB = getNeighbors(newAtoms[-1],radius)
434	while True:
435	if not len(IndB):
436	break
437	indb = IndB.pop()
438	if atomData[Indx[indb]] == None:
439	continue
440	while True:
441	try:
442	jndb = IndB.index(indb)
443	IndB.remove(jndb)
444	except:
445	break
446	newAtom = copy.copy(atomData[Indx[indb]])
447	newAtom[cx:cx+3] = UAtoms[indb] #NB: thermal Uij, etc. not transformed!
448	newAtoms.append(newAtom)
449	atomData[Indx[indb]] = None
450	IndB = set(list(IndB)+list(getNeighbors(newAtoms[-1],radius)))
451	if len(IndB) > nAtom:
452	return 'Assemble molecule cannot be used on extended structures'
453	for atom in atomData:
454	if atom != None:
455	newAtoms.append(atom)
456	return newAtoms
457
458	def FindAtomIndexByIDs(atomData,loc,IDs,Draw=True):
459	'''finds the set of atom array indices for a list of atom IDs. Will search
460	either the Atom table or the drawAtom table.
461
462	:param list atomData: Atom or drawAtom table containting coordinates, etc.
463	:param int loc: location of atom id in atomData record
464	:param list IDs: atom IDs to be found
465	:param bool Draw: True if drawAtom table to be searched; False if Atom table
466	is searched
467
468	:returns: list indx: atom (or drawAtom) indices
469
470	'''
471	indx = []
472	for i,atom in enumerate(atomData):
473	if Draw and atom[loc] in IDs:
474	indx.append(i)
475	elif atom[loc] in IDs:
476	indx.append(i)
477	return indx
478
479	def FillAtomLookUp(atomData,indx):
480	'''create a dictionary of atom indexes with atom IDs as keys
481
482	:param list atomData: Atom table to be used
483
484	:returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
485
486	'''
487	atomLookUp = {}
488	for iatm,atom in enumerate(atomData):
489	atomLookUp[atom[indx]] = iatm
490	return atomLookUp
491
492	def GetAtomsById(atomData,atomLookUp,IdList):
493	'''gets a list of atoms from Atom table that match a set of atom IDs
494
495	:param list atomData: Atom table to be used
496	:param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
497	:param list IdList: atom IDs to be found
498
499	:returns: list atoms: list of atoms found
500
501	'''
502	atoms = []
503	for Id in IdList:
504	atoms.append(atomData[atomLookUp[Id]])
505	return atoms
506
507	def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
508	'''gets atom parameters for atoms using atom IDs
509
510	:param list atomData: Atom table to be used
511	:param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
512	:param list IdList: atom IDs to be found
513	:param int itemLoc: pointer to desired 1st item in an atom table entry
514	:param int numItems: number of items to be retrieved
515
516	:returns: type name: description
517
518	'''
519	Items = []
520	if not isinstance(IdList,list):
521	IdList = [IdList,]
522	for Id in IdList:
523	if numItems == 1:
524	Items.append(atomData[atomLookUp[Id]][itemLoc])
525	else:
526	Items.append(atomData[atomLookUp[Id]][itemLoc:itemLoc+numItems])
527	return Items
528
529	def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
530	'''default doc string
531
532	:param type name: description
533
534	:returns: type name: description
535
536	'''
537	pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
538	dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
539	XYZ = []
540	for ind in indx:
541	names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
542	dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
543	XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
544	return XYZ
545
546	def GetAtomFracByID(pId,parmDict,AtLookup,indx):
547	'''default doc string
548
549	:param type name: description
550
551	:returns: type name: description
552
553	'''
554	pfx = str(pId)+'::Afrac:'
555	Frac = []
556	for ind in indx:
557	name = pfx+str(AtLookup[ind])
558	Frac.append(parmDict[name])
559	return Frac
560
561	# for Atom in Atoms:
562	# XYZ = Atom[cx:cx+3]
563	# if 'A' in Atom[cia]:
564	# U6 = Atom[cia+2:cia+8]
565
566
567	def ApplySeqData(data,seqData):
568	'''Applies result from seq. refinement to drawing atom positions & Uijs
569	'''
570	generalData = data['General']
571	SGData = generalData['SGData']
572	cx,ct,cs,cia = generalData['AtomPtrs']
573	drawingData = data['Drawing']
574	dcx,dct,dcs,dci = drawingData['atomPtrs']
575	atoms = data['Atoms']
576	drawAtoms = drawingData['Atoms']
577	pId = data['pId']
578	pfx = '%d::'%(pId)
579	parmDict = seqData['parmDict']
580	for ia,atom in enumerate(atoms):
581	dxyz = np.array([parmDict[pfx+'dAx:'+str(ia)],parmDict[pfx+'dAy:'+str(ia)],parmDict[pfx+'dAz:'+str(ia)]])
582	if atom[cia] == 'A':
583	atuij = np.array([parmDict[pfx+'AU11:'+str(ia)],parmDict[pfx+'AU22:'+str(ia)],parmDict[pfx+'AU33:'+str(ia)],
584	parmDict[pfx+'AU12:'+str(ia)],parmDict[pfx+'AU13:'+str(ia)],parmDict[pfx+'AU23:'+str(ia)]])
585	else:
586	atuiso = parmDict[pfx+'AUiso:'+str(ia)]
587	atxyz = G2spc.MoveToUnitCell(np.array(atom[cx:cx+3])+dxyz)[0]
588	indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
589	for ind in indx:
590	drawatom = drawAtoms[ind]
591	opr = drawatom[dcs-1]
592	#how do I handle Sfrac? - fade the atoms?
593	if atom[cia] == 'A':
594	X,U = G2spc.ApplyStringOps(opr,SGData,atxyz,atuij)
595	drawatom[dcx:dcx+3] = X
596	drawatom[dci-6:dci] = U
597	else:
598	X = G2spc.ApplyStringOps(opr,SGData,atxyz)
599	drawatom[dcx:dcx+3] = X
600	drawatom[dci-7] = atuiso
601	return drawAtoms
602
603	def FindNeighbors(phase,FrstName,AtNames,notName=''):
604	General = phase['General']
605	cx,ct,cs,cia = General['AtomPtrs']
606	Atoms = phase['Atoms']
607	atNames = [atom[ct-1] for atom in Atoms]
608	Cell = General['Cell'][1:7]
609	Amat,Bmat = G2lat.cell2AB(Cell)
610	atTypes = General['AtomTypes']
611	Radii = np.array(General['BondRadii'])
612	DisAglCtls = General['DisAglCtls']
613	radiusFactor = DisAglCtls['Factors'][0]
614	AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
615	Orig = atNames.index(FrstName)
616	OId = Atoms[Orig][cia+8]
617	OType = Atoms[Orig][ct]
618	XYZ = getAtomXYZ(Atoms,cx)
619	Neigh = []
620	Ids = []
621	Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
622	dist = np.sqrt(np.sum(Dx**2,axis=1))
623	sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
624	IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))
625	for j in IndB[0]:
626	if j != Orig:
627	if AtNames[j] != notName:
628	Neigh.append([AtNames[j],dist[j],True])
629	Ids.append(Atoms[j][cia+8])
630	return Neigh,[OId,Ids]
631
632	def FindAllNeighbors(phase,FrstName,AtNames,notName=''):
633	General = phase['General']
634	cx,ct,cs,cia = General['AtomPtrs']
635	Atoms = phase['Atoms']
636	atNames = [atom[ct-1] for atom in Atoms]
637	Cell = General['Cell'][1:7]
638	Amat,Bmat = G2lat.cell2AB(Cell)
639	SGData = General['SGData']
640	indices = (-1,0,1)
641	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
642	atTypes = General['AtomTypes']
643	Radii = np.array(General['BondRadii'])
644	DisAglCtls = General['DisAglCtls']
645	radiusFactor = DisAglCtls['Factors'][0]
646	AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
647	Orig = atNames.index(FrstName)
648	OId = Atoms[Orig][cia+8]
649	OType = Atoms[Orig][ct]
650	XYZ = getAtomXYZ(Atoms,cx)
651	Oxyz = XYZ[Orig]
652	Neigh = []
653	Ids = []
654	sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
655	sumR = np.reshape(np.tile(sumR,27),(27,-1))
656	results = []
657	for xyz in XYZ:
658	results.append(G2spc.GenAtom(xyz,SGData,False,Move=False))
659	for iA,result in enumerate(results):
660	if iA != Orig:
661	for [Txyz,Top,Tunit,Spn] in result:
662	Dx = np.array([Txyz-Oxyz+unit for unit in Units])
663	dx = np.inner(Dx,Amat)
664	dist = np.sqrt(np.sum(dx**2,axis=1))
665	IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR[:,iA],0.))
666	# GSASIIpath.IPyBreak()
667	for iU in IndB[0]:
668	if AtNames[iA] != notName:
669	unit = Units[iU]
670	if np.any(unit):
671	Topstr = ' +(%4d)[%2d,%2d,%2d]'%(Top,unit[0],unit[1],unit[2])
672	else:
673	Topstr = ' +(%4d)'%(Top)
674	Neigh.append([AtNames[iA]+Topstr,dist[iU]])
675	Ids.append(Atoms[iA][cia+8])
676	return Neigh,[OId,Ids]
677
678	def calcBond(A,Ax,Bx,MTCU):
679	cell = G2lat.A2cell(A)
680	Amat,Bmat = G2lat.cell2AB(cell)
681	M,T,C,U = MTCU
682	Btx = np.inner(M,Bx)+T+C+U
683	Dx = Btx-Ax
684	dist = np.sqrt(np.inner(Amat,Dx))
685	return dist
686
687	def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
688
689	def getTransMat(RXYZ,OXYZ,TXYZ,Amat):
690	Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T
691	Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,nxs]
692	Mat2 = np.cross(Vec[0],Vec[1]) #UxV
693	Mat2 /= np.sqrt(np.sum(Mat2**2))
694	Mat3 = np.cross(Mat2,Vec[0]) #(UxV)xU
695	return nl.inv(np.array([Vec[0],Mat2,Mat3]))
696
697	cx,ct,cs,cia = General['AtomPtrs']
698	Cell = General['Cell'][1:7]
699	Amat,Bmat = G2lat.cell2AB(Cell)
700	nBonds = AddHydId[-1]+len(AddHydId[1])
701	Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0]
702	OXYZ = np.array(Oatom[cx:cx+3])
703	if 'I' in Oatom[cia]:
704	Uiso = Oatom[cia+1]
705	else:
706	Uiso = (Oatom[cia+2]+Oatom[cia+3]+Oatom[cia+4])/3.0 #simple average
707	Uiso = max(Uiso,0.005) #set floor!
708	Tatoms = GetAtomsById(Atoms,AtLookUp,AddHydId[1])
709	TXYZ = np.array([tatom[cx:cx+3] for tatom in Tatoms]) #3 x xyz
710	if nBonds == 4:
711	if AddHydId[-1] == 1:
712	Vec = TXYZ-OXYZ
713	Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0))
714	Vec = np.sum(Vec/Len,axis=0)
715	Len = np.sqrt(np.sum(Vec**2))
716	Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len
717	HU = 1.1*Uiso
718	return [Hpos,],[HU,]
719	elif AddHydId[-1] == 2:
720	Vec = np.inner(Amat,TXYZ-OXYZ).T
721	Vec[0] += Vec[1] #U - along bisector
722	Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,nxs]
723	Mat2 = np.cross(Vec[0],Vec[1]) #UxV
724	Mat2 /= np.sqrt(np.sum(Mat2**2))
725	Mat3 = np.cross(Mat2,Vec[0]) #(UxV)xU
726	iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
727	Hpos = np.array([[-0.97cosd(54.75),0.97sind(54.75),0.],
728	[-0.97cosd(54.75),-0.97sind(54.75),0.]])
729	HU = 1.2Uisonp.ones(2)
730	Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
731	return Hpos,HU
732	else:
733	Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
734	RXYZ = np.array(Ratom[cx:cx+3])
735	iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
736	a = 0.96*cosd(70.5)
737	b = 0.96*sind(70.5)
738	Hpos = np.array([[a,0.,-b],[a,-bcosd(30.),0.5b],[a,bcosd(30.),0.5b]])
739	Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
740	HU = 1.5Uisonp.ones(3)
741	return Hpos,HU
742	elif nBonds == 3:
743	if AddHydId[-1] == 1:
744	Vec = np.sum(TXYZ-OXYZ,axis=0)
745	Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
746	Vec = -0.93*Vec/Len
747	Hpos = OXYZ+Vec
748	HU = 1.1*Uiso
749	return [Hpos,],[HU,]
750	elif AddHydId[-1] == 2:
751	Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
752	RXYZ = np.array(Ratom[cx:cx+3])
753	iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
754	a = 0.93*cosd(60.)
755	b = 0.93*sind(60.)
756	Hpos = [[a,b,0],[a,-b,0]]
757	Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
758	HU = 1.2Uisonp.ones(2)
759	return Hpos,HU
760	else: #2 bonds
761	if 'C' in Oatom[ct]:
762	Vec = TXYZ[0]-OXYZ
763	Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
764	Vec = -0.93*Vec/Len
765	Hpos = OXYZ+Vec
766	HU = 1.1*Uiso
767	return [Hpos,],[HU,]
768	elif 'O' in Oatom[ct]:
769	mapData = General['Map']
770	Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
771	RXYZ = np.array(Ratom[cx:cx+3])
772	iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
773	a = 0.82*cosd(70.5)
774	b = 0.82*sind(70.5)
775	azm = np.arange(0.,360.,5.)
776	Hpos = np.array([[a,bcosd(x),bsind(x)] for x in azm])
777	Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
778	Rhos = np.array([getRho(pos,mapData) for pos in Hpos])
779	imax = np.argmax(Rhos)
780	HU = 1.5*Uiso
781	return [Hpos[imax],],[HU,]
782	return [],[]
783
784	#def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj): #unfinished & not used
785	# '''default doc string
786	#
787	# :param type name: description
788	#
789	# :returns: type name: description
790	#
791	# '''
792	# for atom in atomData:
793	# XYZ = np.inner(Amat,atom[cx:cx+3])
794	# if atom[cia] == 'A':
795	# UIJ = atom[cia+2:cia+8]
796
797	def TLS2Uij(xyz,g,Amat,rbObj): #not used anywhere, but could be?
798	'''default doc string
799
800	:param type name: description
801
802	:returns: type name: description
803
804	'''
805	TLStype,TLS = rbObj['ThermalMotion'][:2]
806	Tmat = np.zeros((3,3))
807	Lmat = np.zeros((3,3))
808	Smat = np.zeros((3,3))
809	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
810	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
811	if 'T' in TLStype:
812	Tmat = G2lat.U6toUij(TLS[:6])
813	if 'L' in TLStype:
814	Lmat = G2lat.U6toUij(TLS[6:12])
815	if 'S' in TLStype:
816	Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
817	XYZ = np.inner(Amat,xyz)
818	Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
819	Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
820	beta = np.inner(np.inner(g,Umat),g)
821	return G2lat.UijtoU6(beta)*gvec
822
823	def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj): #not used anywhere, but could be?
824	'''default doc string
825
826	:param type name: description
827
828	:returns: type name: description
829
830	'''
831	cx,ct,cs,cia = atPtrs
832	TLStype,TLS = rbObj['ThermalMotion'][:2]
833	Tmat = np.zeros((3,3))
834	Lmat = np.zeros((3,3))
835	Smat = np.zeros((3,3))
836	G,g = G2lat.A2Gmat(Amat)
837	gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
838	if 'T' in TLStype:
839	Tmat = G2lat.U6toUij(TLS[:6])
840	if 'L' in TLStype:
841	Lmat = G2lat.U6toUij(TLS[6:12])
842	if 'S' in TLStype:
843	Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
844	for atom in atomData:
845	XYZ = np.inner(Amat,atom[cx:cx+3])
846	Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
847	if 'U' in TLStype:
848	atom[cia+1] = TLS[0]
849	atom[cia] = 'I'
850	else:
851	atom[cia] = 'A'
852	Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
853	beta = np.inner(np.inner(g,Umat),g)
854	atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
855
856	def GetXYZDist(xyz,XYZ,Amat):
857	'''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
858	and are in crystal coordinates; Amat is crystal to Cart matrix
859
860	:param type name: description
861
862	:returns: type name: description
863
864	'''
865	return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
866
867	def getAtomXYZ(atoms,cx):
868	'''default doc string
869
870	:param type name: description
871
872	:returns: type name: description
873
874	'''
875	XYZ = []
876	for atom in atoms:
877	XYZ.append(atom[cx:cx+3])
878	return np.array(XYZ)
879
880	def getRBTransMat(X,Y):
881	'''Get transformation for Cartesian axes given 2 vectors
882	X will be parallel to new X-axis; X cross Y will be new Z-axis &
883	(X cross Y) cross Y will be new Y-axis
884	Useful for rigid body axes definintion
885
886	:param array X: normalized vector
887	:param array Y: normalized vector
888
889	:returns: array M: transformation matrix
890
891	use as XYZ' = np.inner(M,XYZ) where XYZ are Cartesian
892
893	'''
894	Mat2 = np.cross(X,Y) #UxV-->Z
895	Mat2 /= np.sqrt(np.sum(Mat2**2))
896	Mat3 = np.cross(Mat2,X) #(UxV)xU-->Y
897	Mat3 /= np.sqrt(np.sum(Mat3**2))
898	return np.array([X,Mat3,Mat2])
899
900	def RotateRBXYZ(Bmat,Cart,oriQ):
901	'''rotate & transform cartesian coordinates to crystallographic ones
902	no translation applied. To be used for numerical derivatives
903
904	:param type name: description
905
906	:returns: type name: description
907
908	'''
909	''' returns crystal coordinates for atoms described by RBObj
910	'''
911	XYZ = np.zeros_like(Cart)
912	for i,xyz in enumerate(Cart):
913	XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
914	return XYZ
915
916	def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
917	'''default doc string
918
919	:param type name: description
920
921	:returns: type name: description
922
923	'''
924	''' returns crystal coordinates for atoms described by RBObj
925	'''
926	RBRes = RBData[RBType][RBObj['RBId']]
927	if RBType == 'Vector':
928	vecs = RBRes['rbVect']
929	mags = RBRes['VectMag']
930	Cart = np.zeros_like(vecs[0])
931	for vec,mag in zip(vecs,mags):
932	Cart += vec*mag
933	elif RBType == 'Residue':
934	Cart = np.array(RBRes['rbXYZ'])
935	for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
936	QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
937	Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
938	XYZ = np.zeros_like(Cart)
939	for i,xyz in enumerate(Cart):
940	XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
941	return XYZ,Cart
942
943	def UpdateMCSAxyz(Bmat,MCSA):
944	'''default doc string
945
946	:param type name: description
947
948	:returns: type name: description
949
950	'''
951	xyz = []
952	atTypes = []
953	iatm = 0
954	for model in MCSA['Models'][1:]: #skip the MD model
955	if model['Type'] == 'Atom':
956	xyz.append(model['Pos'][0])
957	atTypes.append(model['atType'])
958	iatm += 1
959	else:
960	RBRes = MCSA['rbData'][model['Type']][model['RBId']]
961	Pos = np.array(model['Pos'][0])
962	Ori = np.array(model['Ori'][0])
963	Qori = AVdeg2Q(Ori[0],Ori[1:])
964	if model['Type'] == 'Vector':
965	vecs = RBRes['rbVect']
966	mags = RBRes['VectMag']
967	Cart = np.zeros_like(vecs[0])
968	for vec,mag in zip(vecs,mags):
969	Cart += vec*mag
970	elif model['Type'] == 'Residue':
971	Cart = np.array(RBRes['rbXYZ'])
972	for itor,seq in enumerate(RBRes['rbSeq']):
973	QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
974	Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
975	if model['MolCent'][1]:
976	Cart -= model['MolCent'][0]
977	for i,x in enumerate(Cart):
978	xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
979	atType = RBRes['rbTypes'][i]
980	atTypes.append(atType)
981	iatm += 1
982	return np.array(xyz),atTypes
983
984	def SetMolCent(model,RBData):
985	'''default doc string
986
987	:param type name: description
988
989	:returns: type name: description
990
991	'''
992	rideList = []
993	RBRes = RBData[model['Type']][model['RBId']]
994	if model['Type'] == 'Vector':
995	vecs = RBRes['rbVect']
996	mags = RBRes['VectMag']
997	Cart = np.zeros_like(vecs[0])
998	for vec,mag in zip(vecs,mags):
999	Cart += vec*mag
1000	elif model['Type'] == 'Residue':
1001	Cart = np.array(RBRes['rbXYZ'])
1002	for seq in RBRes['rbSeq']:
1003	rideList += seq[3]
1004	centList = set(range(len(Cart)))-set(rideList)
1005	cent = np.zeros(3)
1006	for i in centList:
1007	cent += Cart[i]
1008	model['MolCent'][0] = cent/len(centList)
1009
1010	def UpdateRBUIJ(Bmat,Cart,RBObj):
1011	'''default doc string
1012
1013	:param type name: description
1014
1015	:returns: type name: description
1016
1017	'''
1018	''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
1019	'''
1020	TLStype,TLS = RBObj['ThermalMotion'][:2]
1021	T = np.zeros(6)
1022	L = np.zeros(6)
1023	S = np.zeros(8)
1024	if 'T' in TLStype:
1025	T = TLS[:6]
1026	if 'L' in TLStype:
1027	L = np.array(TLS[6:12])(np.pi/180.)*2
1028	if 'S' in TLStype:
1029	S = np.array(TLS[12:])*(np.pi/180.)
1030	g = nl.inv(np.inner(Bmat,Bmat))
1031	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
1032	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
1033	Uout = []
1034	Q = RBObj['Orient'][0]
1035	for X in Cart:
1036	X = prodQVQ(Q,X)
1037	if 'U' in TLStype:
1038	Uout.append(['I',TLS[0],0,0,0,0,0,0])
1039	elif not 'N' in TLStype:
1040	U = [0,0,0,0,0,0]
1041	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])
1042	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])
1043	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])
1044	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]+ \
1045	S[4]X[0]-S[5]X[1]-(S[6]+S[7])*X[2]
1046	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]+ \
1047	S[3]X[2]-S[2]X[0]+S[6]*X[1]
1048	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]+ \
1049	S[0]X[1]-S[1]X[2]+S[7]*X[0]
1050	Umat = G2lat.U6toUij(U)
1051	beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
1052	Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
1053	else:
1054	Uout.append(['N',])
1055	return Uout
1056
1057	def GetSHCoeff(pId,parmDict,SHkeys):
1058	'''default doc string
1059
1060	:param type name: description
1061
1062	:returns: type name: description
1063
1064	'''
1065	SHCoeff = {}
1066	for shkey in SHkeys:
1067	shname = str(pId)+'::'+shkey
1068	SHCoeff[shkey] = parmDict[shname]
1069	return SHCoeff
1070
1071	def getMass(generalData):
1072	'''Computes mass of unit cell contents
1073
1074	:param dict generalData: The General dictionary in Phase
1075
1076	:returns: float mass: Crystal unit cell mass in AMU.
1077
1078	'''
1079	mass = 0.
1080	for i,elem in enumerate(generalData['AtomTypes']):
1081	mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
1082	return max(mass,1.0)
1083
1084	def getDensity(generalData):
1085	'''calculate crystal structure density
1086
1087	:param dict generalData: The General dictionary in Phase
1088
1089	:returns: float density: crystal density in gm/cm^3
1090
1091	'''
1092	mass = getMass(generalData)
1093	Volume = generalData['Cell'][7]
1094	density = mass/(0.6022137*Volume)
1095	return density,Volume/mass
1096
1097	def getWave(Parms):
1098	'''returns wavelength from Instrument parameters dictionary
1099
1100	:param dict Parms: Instrument parameters;
1101	must contain:
1102	Lam: single wavelength
1103	or
1104	Lam1: Ka1 radiation wavelength
1105
1106	:returns: float wave: wavelength
1107
1108	'''
1109	try:
1110	return Parms['Lam'][1]
1111	except KeyError:
1112	return Parms['Lam1'][1]
1113
1114	def getMeanWave(Parms):
1115	'''returns mean wavelength from Instrument parameters dictionary
1116
1117	:param dict Parms: Instrument parameters;
1118	must contain:
1119	Lam: single wavelength
1120	or
1121	Lam1,Lam2: Ka1,Ka2 radiation wavelength
1122	I(L2)/I(L1): Ka2/Ka1 ratio
1123
1124	:returns: float wave: mean wavelength
1125
1126	'''
1127	try:
1128	return Parms['Lam'][1]
1129	except KeyError:
1130	meanLam = (Parms['Lam1'][1]+Parms['I(L2)/I(L1)'][1]*Parms['Lam2'][1])/ \
1131	(1.+Parms['I(L2)/I(L1)'][1])
1132	return meanLam
1133
1134
1135	def El2Mass(Elements):
1136	'''compute molecular weight from Elements
1137
1138	:param dict Elements: elements in molecular formula;
1139	each must contain
1140	Num: number of atoms in formula
1141	Mass: at. wt.
1142
1143	:returns: float mass: molecular weight.
1144
1145	'''
1146	mass = 0
1147	for El in Elements:
1148	mass += Elements[El]['Num']*Elements[El]['Mass']
1149	return mass
1150
1151	def Den2Vol(Elements,density):
1152	'''converts density to molecular volume
1153
1154	:param dict Elements: elements in molecular formula;
1155	each must contain
1156	Num: number of atoms in formula
1157	Mass: at. wt.
1158	:param float density: material density in gm/cm^3
1159
1160	:returns: float volume: molecular volume in A^3
1161
1162	'''
1163	return El2Mass(Elements)/(density*0.6022137)
1164
1165	def Vol2Den(Elements,volume):
1166	'''converts volume to density
1167
1168	:param dict Elements: elements in molecular formula;
1169	each must contain
1170	Num: number of atoms in formula
1171	Mass: at. wt.
1172	:param float volume: molecular volume in A^3
1173
1174	:returns: float density: material density in gm/cm^3
1175
1176	'''
1177	return El2Mass(Elements)/(volume*0.6022137)
1178
1179	def El2EstVol(Elements):
1180	'''Estimate volume from molecular formula; assumes atom volume = 10A^3
1181
1182	:param dict Elements: elements in molecular formula;
1183	each must contain
1184	Num: number of atoms in formula
1185
1186	:returns: float volume: estimate of molecular volume in A^3
1187
1188	'''
1189	vol = 0
1190	for El in Elements:
1191	vol += 10.*Elements[El]['Num']
1192	return vol
1193
1194	def XScattDen(Elements,vol,wave=0.):
1195	'''Estimate X-ray scattering density from molecular formula & volume;
1196	ignores valence, but includes anomalous effects
1197
1198	:param dict Elements: elements in molecular formula;
1199	each element must contain
1200	Num: number of atoms in formula
1201	Z: atomic number
1202	:param float vol: molecular volume in A^3
1203	:param float wave: optional wavelength in A
1204
1205	:returns: float rho: scattering density in 10^10cm^-2;
1206	if wave > 0 the includes f' contribution
1207	:returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
1208	:returns: float fpp: if wave>0 f" in 10^10cm^-2; otherwise 0
1209
1210	'''
1211	rho = 0
1212	mu = 0
1213	fpp = 0
1214	if wave:
1215	Xanom = XAnomAbs(Elements,wave)
1216	for El in Elements:
1217	f0 = Elements[El]['Z']
1218	if wave:
1219	f0 += Xanom[El][0]
1220	fpp += Xanom[El][1]*Elements[El]['Num']
1221	mu += Xanom[El][2]*Elements[El]['Num']
1222	rho += Elements[El]['Num']*f0
1223	return 28.179rho/vol,0.1mu/vol,28.179*fpp/vol
1224
1225	def NCScattDen(Elements,vol,wave=0.):
1226	'''Estimate neutron scattering density from molecular formula & volume;
1227	ignores valence, but includes anomalous effects
1228
1229	:param dict Elements: elements in molecular formula;
1230	each element must contain
1231	Num: number of atoms in formula
1232	Z: atomic number
1233	:param float vol: molecular volume in A^3
1234	:param float wave: optional wavelength in A
1235
1236	:returns: float rho: scattering density in 10^10cm^-2;
1237	if wave > 0 the includes f' contribution
1238	:returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
1239	:returns: float fpp: if wave>0 f" in 10^10cm^-2; otherwise 0
1240
1241	'''
1242	rho = 0
1243	mu = 0
1244	bpp = 0
1245	for El in Elements:
1246	isotope = Elements[El]['Isotope']
1247	b0 = Elements[El]['Isotopes'][isotope]['SL'][0]
1248	mu += Elements[El]['Isotopes'][isotope].get('SA',0.)*Elements[El]['Num']
1249	if wave and 'BW-LS' in Elements[El]['Isotopes'][isotope]:
1250	Re,Im,E0,gam,A,E1,B,E2 = Elements[El]['Isotopes'][isotope]['BW-LS'][1:]
1251	Emev = 81.80703/wave**2
1252	T0 = Emev-E0
1253	T1 = Emev-E1
1254	T2 = Emev-E2
1255	D0 = T02+gam2
1256	D1 = T12+gam2
1257	D2 = T22+gam2
1258	b0 += Re(T0/D0+AT1/D1+B*T2/D2)
1259	bpp += Im*(1/D0+A/D1+B/D2)
1260	else:
1261	bpp += Elements[El]['Isotopes'][isotope]['SL'][1]
1262	rho += Elements[El]['Num']*b0
1263	if wave: mu *= wave
1264	return 100.rho/vol,mu/vol,100.bpp/vol
1265
1266	def wavekE(wavekE):
1267	'''Convert wavelength to energy & vise versa
1268
1269	:param float waveKe:wavelength in A or energy in kE
1270
1271	:returns float waveKe:the other one
1272
1273	'''
1274	return 12.397639/wavekE
1275
1276	def XAnomAbs(Elements,wave):
1277	kE = wavekE(wave)
1278	Xanom = {}
1279	for El in Elements:
1280	Orbs = G2el.GetXsectionCoeff(El)
1281	Xanom[El] = G2el.FPcalc(Orbs, kE)
1282	return Xanom #f',f", mu
1283
1284	################################################################################
1285	#### Modulation math
1286	################################################################################
1287
1288	def makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast):
1289	'''
1290	waveTypes: array nAtoms: 'Fourier','ZigZag' or 'Block'
1291	FSSdata: array 2 x atoms x waves (sin,cos terms)
1292	XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1293	USSdata: array 2x6 x atoms X waves (sin,cos terms)
1294	MSSdata: array 2x3 x atoms X waves (sin,cos terms)
1295
1296	Mast: array orthogonalization matrix for Uij
1297	'''
1298	ngl = 36 #selected for integer steps for 1/6,1/4,1/3...
1299	glTau,glWt = pwd.pygauleg(0.,1.,ngl) #get Gauss-Legendre intervals & weights
1300	mglTau = np.arange(0.,1.,1./ngl)
1301	Ax = np.array(XSSdata[:3]).T #atoms x waves x sin pos mods
1302	Bx = np.array(XSSdata[3:]).T #...cos pos mods
1303	Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms
1304	Bf = np.array(FSSdata[1]).T #cos frac mods...
1305	Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods as betaij
1306	Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods as betaij
1307	Am = np.array(MSSdata[:3]).T #atoms x waves x sin pos mods
1308	Bm = np.array(MSSdata[3:]).T #...cos pos mods
1309	nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1],Am.shape[1]]
1310	if nWaves[0]:
1311	tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau #Fwaves x ngl
1312	FmodA = Af[:,:,nxs]np.sin(twopitauF[nxs,:,:]) #atoms X Fwaves X ngl
1313	FmodB = Bf[:,:,nxs]np.cos(twopitauF[nxs,:,:])
1314	Fmod = np.sum(1.0+FmodA+FmodB,axis=1) #atoms X ngl; sum waves
1315	else:
1316	Fmod = 1.0
1317	XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1318	XmodA = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1319	XmodB = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1320	for iatm in range(Ax.shape[0]):
1321	nx = 0
1322	if 'ZigZag' in waveTypes[iatm]:
1323	nx = 1
1324	Tmm = Ax[iatm][0][:2]
1325	XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
1326	XmodZ[iatm][0] += posZigZag(glTau,Tmm,XYZmax).T
1327	elif 'Block' in waveTypes[iatm]:
1328	nx = 1
1329	Tmm = Ax[iatm][0][:2]
1330	XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])
1331	XmodZ[iatm][0] += posBlock(glTau,Tmm,XYZmax).T
1332	tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau #Xwaves x ngl
1333	if nx:
1334	XmodA[iatm][:-nx] = Ax[iatm,nx:,:,nxs]np.sin(twopitauX[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1335	XmodB[iatm][:-nx] = Bx[iatm,nx:,:,nxs]np.cos(twopitauX[nxs,:,nxs,:]) #ditto
1336	else:
1337	XmodA[iatm] = Ax[iatm,:,:,nxs]np.sin(twopitauX[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1338	XmodB[iatm] = Bx[iatm,:,:,nxs]np.cos(twopitauX[nxs,:,nxs,:]) #ditto
1339	Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X 3 X ngl; sum waves
1340	Xmod = np.swapaxes(Xmod,1,2)
1341	if nWaves[2]:
1342	tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau #Uwaves x ngl
1343	UmodA = Au[:,:,:,:,nxs]np.sin(twopitauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x ngl
1344	UmodB = Bu[:,:,:,:,nxs]np.cos(twopitauU[nxs,:,nxs,nxs,:]) #ditto
1345	Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3) #atoms x 3x3 x ngl; sum waves
1346	else:
1347	Umod = 1.0
1348	if nWaves[3]:
1349	tauM = np.arange(1.,nWaves[3]+1-nx)[:,nxs]*mglTau #Mwaves x ngl
1350	MmodA = Am[:,:,:,nxs]np.sin(twopitauM[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
1351	MmodB = Bm[:,:,:,nxs]np.cos(twopitauM[nxs,:,nxs,:]) #ditto
1352	Mmod = np.sum(MmodA+MmodB,axis=1)
1353	Mmod = np.swapaxes(Mmod,1,2) #Mxyz,Ntau,Natm
1354	else:
1355	Mmod = 1.0
1356	return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt
1357
1358	def MagMod(ngl,XYZ,modQ,MSSdata,SGData,SSGData):
1359	'''
1360	this needs to make magnetic moment modulations & magnitudes as
1361	fxn of ngl tau points
1362	'''
1363	Bm = np.array(MSSdata[:3]).T #atoms x waves x sin pos mods
1364	Am = np.array(MSSdata[3:]).T #...cos pos mods
1365	nWaves = Am.shape[1]
1366	tau = np.arange(ngl)/ngl
1367	if not nWaves:
1368	return 0.0,0.0
1369	SGMT = np.array([ops[0] for ops in SGData['SGOps']]) #not .T!!
1370	Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
1371	SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1372	if SGData['SGInv']:
1373	SGMT = np.vstack((SGMT,-SGMT))
1374	Sinv =np.vstack((Sinv,-Sinv))
1375	SGT = np.vstack((SGT,-SGT))
1376	SGMT = np.vstack([SGMT for cen in SGData['SGCen']])
1377	Sinv = np.vstack([Sinv for cen in SGData['SGCen']])
1378	SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
1379	detSM = nl.det(SGMT)
1380	mst = Sinv[:,3,:3]
1381	epsinv = Sinv[:,3,3]
1382	phi = np.inner(XYZ,modQ).T
1383	TA = np.sum(mst[nxs,:,:]*(XYZ-SGT[:,:3][nxs,:,:]),axis=-1).T
1384	tauT = TA[nxs,:,:] + epsinv[nxs,:,nxs]*(tau[:,nxs,nxs]-SGT[:,3][nxs,:,nxs]+phi[nxs,:,:])
1385	modind = np.arange(nWaves)+1.
1386	phase = (modind[:,nxs,nxs]*tauT) #Nops,Natm,Nwave
1387	psin = np.sin(twopi*phase)
1388	pcos = np.cos(twopi*phase)
1389	MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)
1390	MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)
1391	if SGData['SGGray']:
1392	MmodA = -np.sum(SGMT[nxs,:,nxs,:,:]MmodA[:,:,:,nxs,:],axis=-1)detSM[nxs,:,nxs,nxs]
1393	MmodB = -np.sum(SGMT[nxs,:,nxs,:,:]MmodB[:,:,:,nxs,:],axis=-1)detSM[nxs,:,nxs,nxs]
1394	else:
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	Mmod = MmodA+MmodB
1398	return Mmod,MmodA,MmodB #Ntau,Nops,Natm,,Mxyz; sum,Re & Im parts
1399
1400	def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1401	'''
1402	H: array nRefBlk x ops X hklt
1403	HP: array nRefBlk x ops X hklt proj to hkl
1404	nWaves: list number of waves for frac, pos, uij & mag
1405	Fmod: array 2 x atoms x waves (sin,cos terms)
1406	Xmod: array atoms X 3 X ngl
1407	Umod: array atoms x 3x3 x ngl
1408	glTau,glWt: arrays Gauss-Lorentzian pos & wts
1409	'''
1410
1411	if nWaves[2]: #uij (adp) waves
1412	if len(HP.shape) > 2:
1413	HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1414	else:
1415	HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1416	else:
1417	HbH = 1.0
1418	HdotX = np.inner(HP,Xmod) #refBlk x ops x atoms X ngl
1419	if len(H.shape) > 2:
1420	D = H[:,:,3:]glTau[nxs,nxs,:] #me*tau; refBlk x ops X ngl
1421	HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1422	else:
1423	D = H[:,3:]glTau[nxs,:] #me*tau; refBlk x ops X ngl
1424	HdotXD = twopi*(HdotX+D[:,nxs,:])
1425	cosHA = np.sum(FmodHbHnp.cos(HdotXD)*glWt,axis=-1) #real part; refBlk X ops x atoms; sum for G-L integration
1426	sinHA = np.sum(FmodHbHnp.sin(HdotXD)*glWt,axis=-1) #imag part; ditto
1427	return np.array([cosHA,sinHA]) # 2 x refBlk x SGops x atoms
1428
1429	def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
1430	'''
1431	H: array nRefBlk x ops X hklt
1432	HP: array nRefBlk x ops X hklt proj to hkl
1433	nWaves: list number of waves for frac, pos, uij & mag
1434	Fmod: array 2 x atoms x waves (sin,cos terms)
1435	Xmod: array atoms X 3 X ngl
1436	Umod: 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	print('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	print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(result[2]['fvec']),len(varyList)))
3010	print ('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	print ('* Least squares aborted - some invalid esds possible *')
3015	break #refinement succeeded - finish up!
3016	except ValueError: #result[1] is None on singular matrix
3017	print ('** Refinement failed - singular matrix **')
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	print ('**** ERROR - Fourier map not defined')
3065	return
3066	mapData = generalData['Map']
3067	dmin = mapData['Resolution']
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	print ('Omit map time: %.4f no. elements: %d'%(time.time()-time0,Fhkl.size))
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['Resolution']
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	print ('Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size))
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['Resolution']
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	print ('Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size))
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	print (' 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	print (' 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['Resolution']
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	print (' 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	print (' 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	print (' 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	print (' 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['Resolution']
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	print (' 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	print (' 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	param xyz: array coordinates to be probed Nx3
3655	param: rho: array copy of map (NB: don't use original!)
3656
3657	:returns: density at xyz
3658	'''
3659	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)[:2,:2,:2]
3660	mapShape = np.array(rho.shape)
3661	mapStep = 1./mapShape
3662	X = XYZ%1. #get into unit cell
3663	I = np.array(np.rint(X*mapShape),dtype='int')
3664	R = np.zeros(len(XYZ))
3665	for i,x in enumerate(X):
3666	D = x-I[i]*mapStep #position inside map cell
3667	Rho = rollMap(rho,-I[i]) #shifts map so point is in corner
3668	RIJ = Rho[0,:2,:2]*(1.-D[0])
3669	RI = RIJ[0](1.-D[1])+RIJ[1]D[1]
3670	R[i] = RI[0](1.-D[2])+RI[1]D[2]
3671	return R
3672
3673	def SearchMap(generalData,drawingData,Neg=False):
3674	'''Does a search of a density map for peaks meeting the criterion of peak
3675	height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3676	mapData is data['General']['mapData']; the map is also in mapData.
3677
3678	:param generalData: the phase data structure; includes the map
3679	:param drawingData: the drawing data structure
3680	:param Neg: if True then search for negative peaks (i.e. H-atoms & neutron data)
3681
3682	:returns: (peaks,mags,dzeros) where
3683
3684	* peaks : ndarray
3685	x,y,z positions of the peaks found in the map
3686	* mags : ndarray
3687	the magnitudes of the peaks
3688	* dzeros : ndarray
3689	the distance of the peaks from the unit cell origin
3690	* dcent : ndarray
3691	the distance of the peaks from the unit cell center
3692
3693	'''
3694	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3695
3696	norm = 1./(np.sqrt(3.)np.sqrt(2.np.pi)**3)
3697
3698	def fixSpecialPos(xyz,SGData,Amat):
3699	equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3700	X = []
3701	xyzs = [equiv[0] for equiv in equivs]
3702	for x in xyzs:
3703	if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3704	X.append(x)
3705	if len(X) > 1:
3706	return np.average(X,axis=0)
3707	else:
3708	return xyz
3709
3710	def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3711	Mag,x0,y0,z0,sig = parms
3712	z = -((x0-rX)2+(y0-rY)2+(z0-rZ)*2)/(2.sig**2)
3713	# return normMagnp.exp(z)/(sigres*3) #not slower but some faults in LS
3714	return normMag(1.+z+z*2/2.)/(sigres**3)
3715
3716	def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3717	Mag,x0,y0,z0,sig = parms
3718	M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3719	return M
3720
3721	def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3722	Mag,x0,y0,z0,sig = parms
3723	dMdv = np.zeros(([5,]+list(rX.shape)))
3724	delt = .01
3725	for i in range(5):
3726	parms[i] -= delt
3727	rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3728	parms[i] += 2.*delt
3729	rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3730	parms[i] -= delt
3731	dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3732	rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3733	Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3734	dMdv = np.reshape(dMdv,(5,rX.size))
3735	Hess = np.inner(dMdv,dMdv)
3736
3737	return Vec,Hess
3738
3739	SGData = generalData['SGData']
3740	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3741	peaks = []
3742	mags = []
3743	dzeros = []
3744	dcent = []
3745	try:
3746	mapData = generalData['Map']
3747	contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3748	if Neg:
3749	rho = -copy.copy(mapData['rho']) #flip +/-
3750	else:
3751	rho = copy.copy(mapData['rho']) #don't mess up original
3752	mapHalf = np.array(rho.shape)/2
3753	res = mapData['Resolution']
3754	incre = np.array(rho.shape,dtype=np.float)
3755	step = max(1.0,1./res)+1
3756	steps = np.array((3*[step,]),dtype='int32')
3757	except KeyError:
3758	print ('**** ERROR - Fourier map not defined')
3759	return peaks,mags
3760	rhoMask = ma.array(rho,mask=(rho<contLevel))
3761	indices = (-1,0,1)
3762	rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3763	for roll in rolls:
3764	if np.any(roll):
3765	rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3766	indx = np.transpose(rhoMask.nonzero())
3767	peaks = indx/incre
3768	mags = rhoMask[rhoMask.nonzero()]
3769	for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3770	rho = rollMap(rho,ind)
3771	rMM = mapHalf-steps
3772	rMP = mapHalf+steps+1
3773	rhoPeak = rho[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3774	peakInt = np.sum(rhoPeak)res*3
3775	rX,rY,rZ = np.mgrid[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3776	x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0] #magnitude, position & width(sig)
3777	result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3778	args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3779	x1 = result[0]
3780	if not np.any(x1 < 0):
3781	peak = (np.array(x1[1:4])-ind)/incre
3782	peak = fixSpecialPos(peak,SGData,Amat)
3783	rho = rollMap(rho,-ind)
3784	cent = np.ones(3)*.5
3785	dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3786	dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3787	if Neg: #want negative magnitudes for negative peaks
3788	return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3789	else:
3790	return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3791
3792	def sortArray(data,pos,reverse=False):
3793	'''data is a list of items
3794	sort by pos in list; reverse if True
3795	'''
3796	T = []
3797	for i,M in enumerate(data):
3798	try:
3799	T.append((M[pos],i))
3800	except IndexError:
3801	return data
3802	D = dict(zip(T,data))
3803	T.sort()
3804	if reverse:
3805	T.reverse()
3806	X = []
3807	for key in T:
3808	X.append(D[key])
3809	return X
3810
3811	def PeaksEquiv(data,Ind):
3812	'''Find the equivalent map peaks for those selected. Works on the
3813	contents of data['Map Peaks'].
3814
3815	:param data: the phase data structure
3816	:param list Ind: list of selected peak indices
3817	:returns: augmented list of peaks including those related by symmetry to the
3818	ones in Ind
3819
3820	'''
3821	def Duplicate(xyz,peaks,Amat):
3822	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3823	return True
3824	return False
3825
3826	generalData = data['General']
3827	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3828	SGData = generalData['SGData']
3829	mapPeaks = data['Map Peaks']
3830	XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3831	Indx = {}
3832	for ind in Ind:
3833	xyz = np.array(mapPeaks[ind][1:4])
3834	xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3835	for jnd,xyz in enumerate(XYZ):
3836	Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3837	Ind = []
3838	for ind in Indx:
3839	if Indx[ind]:
3840	Ind.append(ind)
3841	return Ind
3842
3843	def PeaksUnique(data,Ind):
3844	'''Finds the symmetry unique set of peaks from those selected. Works on the
3845	contents of data['Map Peaks'].
3846
3847	:param data: the phase data structure
3848	:param list Ind: list of selected peak indices
3849	:returns: the list of symmetry unique peaks from among those given in Ind
3850
3851	'''
3852	# XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3853
3854	def noDuplicate(xyz,peaks,Amat):
3855	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3856	return False
3857	return True
3858
3859	generalData = data['General']
3860	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3861	SGData = generalData['SGData']
3862	mapPeaks = data['Map Peaks']
3863	Indx = {}
3864	XYZ = {}
3865	for ind in Ind:
3866	XYZ[ind] = np.array(mapPeaks[ind][1:4])
3867	Indx[ind] = True
3868	for ind in Ind:
3869	if Indx[ind]:
3870	xyz = XYZ[ind]
3871	for jnd in Ind:
3872	if ind != jnd and Indx[jnd]:
3873	Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3874	xyzs = np.array([equiv[0] for equiv in Equiv])
3875	Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3876	Ind = []
3877	for ind in Indx:
3878	if Indx[ind]:
3879	Ind.append(ind)
3880	return Ind
3881
3882	################################################################################
3883	##### single peak fitting profile fxn stuff
3884	################################################################################
3885
3886	def getCWsig(ins,pos):
3887	'''get CW peak profile sigma^2
3888
3889	:param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3890	as values only
3891	:param float pos: 2-theta of peak
3892	:returns: float getCWsig: peak sigma^2
3893
3894	'''
3895	tp = tand(pos/2.0)
3896	return ins['U']tp2+ins['V']tp+ins['W']
3897
3898	def getCWsigDeriv(pos):
3899	'''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3900
3901	:param float pos: 2-theta of peak
3902
3903	:returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3904
3905	'''
3906	tp = tand(pos/2.0)
3907	return tp**2,tp,1.0
3908
3909	def getCWgam(ins,pos):
3910	'''get CW peak profile gamma
3911
3912	:param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3913	as values only
3914	:param float pos: 2-theta of peak
3915	:returns: float getCWgam: peak gamma
3916
3917	'''
3918	return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)+ins['Z']
3919
3920	def getCWgamDeriv(pos):
3921	'''get derivatives of CW peak profile gamma wrt X, Y & Z
3922
3923	:param float pos: 2-theta of peak
3924
3925	:returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3926
3927	'''
3928	return 1./cosd(pos/2.0),tand(pos/2.0),1.0
3929
3930	def getTOFsig(ins,dsp):
3931	'''get TOF peak profile sigma^2
3932
3933	:param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3934	as values only
3935	:param float dsp: d-spacing of peak
3936
3937	:returns: float getTOFsig: peak sigma^2
3938
3939	'''
3940	return ins['sig-0']+ins['sig-1']dsp2+ins['sig-2']dsp*4+ins['sig-q']dsp
3941
3942	def getTOFsigDeriv(dsp):
3943	'''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3944
3945	:param float dsp: d-spacing of peak
3946
3947	:returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3948
3949	'''
3950	return 1.0,dsp2,dsp4,dsp
3951
3952	def getTOFgamma(ins,dsp):
3953	'''get TOF peak profile gamma
3954
3955	:param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3956	as values only
3957	:param float dsp: d-spacing of peak
3958
3959	:returns: float getTOFgamma: peak gamma
3960
3961	'''
3962	return ins['Z']+ins['X']dsp+ins['Y']dsp**2
3963
3964	def getTOFgammaDeriv(dsp):
3965	'''get derivatives of TOF peak profile gamma wrt X, Y & Z
3966
3967	:param float dsp: d-spacing of peak
3968
3969	:returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3970
3971	'''
3972	return dsp,dsp**2,1.0
3973
3974	def getTOFbeta(ins,dsp):
3975	'''get TOF peak profile beta
3976
3977	:param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3978	as values only
3979	:param float dsp: d-spacing of peak
3980
3981	:returns: float getTOFbeta: peak beat
3982
3983	'''
3984	return ins['beta-0']+ins['beta-1']/dsp4+ins['beta-q']/dsp2
3985
3986	def getTOFbetaDeriv(dsp):
3987	'''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3988
3989	:param float dsp: d-spacing of peak
3990
3991	:returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3992
3993	'''
3994	return 1.0,1./dsp4,1./dsp2
3995
3996	def getTOFalpha(ins,dsp):
3997	'''get TOF peak profile alpha
3998
3999	:param dict ins: instrument parameters with at least 'alpha'
4000	as values only
4001	:param float dsp: d-spacing of peak
4002
4003	:returns: flaot getTOFalpha: peak alpha
4004
4005	'''
4006	return ins['alpha']/dsp
4007
4008	def getTOFalphaDeriv(dsp):
4009	'''get derivatives of TOF peak profile beta wrt alpha
4010
4011	:param float dsp: d-spacing of peak
4012
4013	:returns: float getTOFalphaDeriv: d(alp)/d(alpha)
4014
4015	'''
4016	return 1./dsp
4017
4018	def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
4019	'''set starting peak parameters for single peak fits from plot selection or auto selection
4020
4021	:param dict Parms: instrument parameters dictionary
4022	:param dict Parms2: table lookup for TOF profile coefficients
4023	:param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
4024	:param float mag: peak top magnitude from pick
4025	:param bool ifQ: True if pos in Q
4026	:param bool useFit: True if use fitted CW Parms values (not defaults)
4027
4028	:returns: list XY: peak list entry:
4029	for CW: [pos,0,mag,1,sig,0,gam,0]
4030	for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4031	NB: mag refinement set by default, all others off
4032
4033	'''
4034	ind = 0
4035	if useFit:
4036	ind = 1
4037	ins = {}
4038	if 'C' in Parms['Type'][0]: #CW data - TOF later in an else
4039	for x in ['U','V','W','X','Y','Z']:
4040	ins[x] = Parms[x][ind]
4041	if ifQ: #qplot - convert back to 2-theta
4042	pos = 2.0asind(posgetWave(Parms)/(4*math.pi))
4043	sig = getCWsig(ins,pos)
4044	gam = getCWgam(ins,pos)
4045	XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st
4046	else:
4047	if ifQ:
4048	dsp = 2.*np.pi/pos
4049	pos = Parms['difC']*dsp
4050	else:
4051	dsp = pos/Parms['difC'][1]
4052	if 'Pdabc' in Parms2:
4053	for x in ['sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4054	ins[x] = Parms[x][ind]
4055	Pdabc = Parms2['Pdabc'].T
4056	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
4057	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
4058	else:
4059	for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4060	ins[x] = Parms[x][ind]
4061	alp = getTOFalpha(ins,dsp)
4062	bet = getTOFbeta(ins,dsp)
4063	sig = getTOFsig(ins,dsp)
4064	gam = getTOFgamma(ins,dsp)
4065	XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4066	return XY
4067
4068	################################################################################
4069	##### MC/SA stuff
4070	################################################################################
4071
4072	#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
4073	# Original Author: Travis Oliphant 2002
4074	# Bug-fixes in 2006 by Tim Leslie
4075
4076
4077	import numpy
4078	from numpy import asarray, exp, squeeze, sign, \
4079	all, shape, array, where
4080	from numpy import random
4081
4082	#__all__ = ['anneal']
4083
4084	_double_min = numpy.finfo(float).min
4085	_double_max = numpy.finfo(float).max
4086	class base_schedule(object):
4087	def __init__(self):
4088	self.dwell = 20
4089	self.lower = 0.
4090	self.upper = 1.
4091	self.Ninit = 50
4092	self.accepted = 0
4093	self.tests = 0
4094	self.feval = 0
4095	self.k = 0
4096	self.T = None
4097
4098	def init(self, **options):
4099	self.__dict__.update(options)
4100	self.lower = asarray(self.lower)
4101	self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
4102	self.upper = asarray(self.upper)
4103	self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
4104	self.k = 0
4105	self.accepted = 0
4106	self.feval = 0
4107	self.tests = 0
4108
4109	def getstart_temp(self, best_state):
4110	""" Find a matching starting temperature and starting parameters vector
4111	i.e. find x0 such that func(x0) = T0.
4112
4113	:param best_state: _state
4114	A _state object to store the function value and x0 found.
4115
4116	:returns: x0 : array
4117	The starting parameters vector.
4118	"""
4119
4120	assert(not self.dims is None)
4121	lrange = self.lower
4122	urange = self.upper
4123	fmax = _double_min
4124	fmin = _double_max
4125	for _ in range(self.Ninit):
4126	x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
4127	fval = self.func(x0, *self.args)
4128	self.feval += 1
4129	if fval > fmax:
4130	fmax = fval
4131	if fval < fmin:
4132	fmin = fval
4133	best_state.cost = fval
4134	best_state.x = array(x0)
4135
4136	self.T0 = (fmax-fmin)*1.5
4137	return best_state.x
4138
4139	def set_range(self,x0,frac):
4140	delrange = frac*(self.upper-self.lower)
4141	self.upper = x0+delrange
4142	self.lower = x0-delrange
4143
4144	def accept_test(self, dE):
4145	T = self.T
4146	self.tests += 1
4147	if dE < 0:
4148	self.accepted += 1
4149	return 1
4150	p = exp(-dE*1.0/T)
4151	if (p > random.uniform(0.0, 1.0)):
4152	self.accepted += 1
4153	return 1
4154	return 0
4155
4156	def update_guess(self, x0):
4157	return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
4158
4159	def update_temp(self, x0):
4160	pass
4161
4162	class fast_sa(base_schedule):
4163	def init(self, **options):
4164	self.__dict__.update(options)
4165
4166	def update_guess(self, x0):
4167	x0 = asarray(x0)
4168	u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4169	T = self.T
4170	xc = (sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)+1.0)/2.0
4171	xnew = xc*(self.upper - self.lower)+self.lower
4172	return xnew
4173
4174	def update_temp(self):
4175	self.T = self.T0exp(-self.c self.k**(self.quench))
4176	self.k += 1
4177	return
4178
4179	class log_sa(base_schedule): #OK
4180
4181	def init(self,**options):
4182	self.__dict__.update(options)
4183
4184	def update_guess(self,x0): #same as default #TODO - is this a reasonable update procedure?
4185	u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4186	T = self.T
4187	xc = (sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)+1.0)/2.0
4188	xnew = xc*(self.upper - self.lower)+self.lower
4189	return xnew
4190
4191	def update_temp(self):
4192	self.k += 1
4193	self.T = self.T0self.slope*self.k
4194
4195	class _state(object):
4196	def __init__(self):
4197	self.x = None
4198	self.cost = None
4199
4200	def makeTsched(data):
4201	if data['Algorithm'] == 'fast':
4202	sched = fast_sa()
4203	sched.quench = data['fast parms'][0]
4204	sched.c = data['fast parms'][1]
4205	elif data['Algorithm'] == 'log':
4206	sched = log_sa()
4207	sched.slope = data['log slope']
4208	sched.T0 = data['Annealing'][0]
4209	if not sched.T0:
4210	sched.T0 = 50.
4211	Tf = data['Annealing'][1]
4212	if not Tf:
4213	Tf = 0.001
4214	Tsched = [sched.T0,]
4215	while Tsched[-1] > Tf:
4216	sched.update_temp()
4217	Tsched.append(sched.T)
4218	return Tsched[1:]
4219
4220	def anneal(func, x0, args=(), schedule='fast',
4221	T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4222	feps=1e-6, quench=1.0, c=1.0,
4223	lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4224	ranRange=0.10,autoRan=False,dlg=None):
4225	"""Minimize a function using simulated annealing.
4226
4227	Schedule is a schedule class implementing the annealing schedule.
4228	Available ones are 'fast', 'cauchy', 'boltzmann'
4229
4230	:param callable func: f(x, \*args)
4231	Function to be optimized.
4232	:param ndarray x0:
4233	Initial guess.
4234	:param tuple args:
4235	Extra parameters to `func`.
4236	:param base_schedule schedule:
4237	Annealing schedule to use (a class).
4238	:param float T0:
4239	Initial Temperature (estimated as 1.2 times the largest
4240	cost-function deviation over random points in the range).
4241	:param float Tf:
4242	Final goal temperature.
4243	:param int maxeval:
4244	Maximum function evaluations.
4245	:param int maxaccept:
4246	Maximum changes to accept.
4247	:param int maxiter:
4248	Maximum cooling iterations.
4249	:param float feps:
4250	Stopping relative error tolerance for the function value in
4251	last four coolings.
4252	:param float quench,c:
4253	Parameters to alter fast_sa schedule.
4254	:param float/ndarray lower,upper:
4255	Lower and upper bounds on `x`.
4256	:param int dwell:
4257	The number of times to search the space at each temperature.
4258	:param float slope:
4259	Parameter for log schedule
4260	:param bool ranStart=False:
4261	True for set 10% of ranges about x
4262
4263	:returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4264
4265	* xmin (ndarray): Point giving smallest value found.
4266	* Jmin (float): Minimum value of function found.
4267	* T (float): Final temperature.
4268	* feval (int): Number of function evaluations.
4269	* iters (int): Number of cooling iterations.
4270	* accept (int): Number of tests accepted.
4271	* retval (int): Flag indicating stopping condition:
4272
4273	* 0: Points no longer changing
4274	* 1: Cooled to final temperature
4275	* 2: Maximum function evaluations
4276	* 3: Maximum cooling iterations reached
4277	* 4: Maximum accepted query locations reached
4278	* 5: Final point not the minimum amongst encountered points
4279
4280	Notes:
4281	Simulated annealing is a random algorithm which uses no derivative
4282	information from the function being optimized. In practice it has
4283	been more useful in discrete optimization than continuous
4284	optimization, as there are usually better algorithms for continuous
4285	optimization problems.
4286
4287	Some experimentation by trying the difference temperature
4288	schedules and altering their parameters is likely required to
4289	obtain good performance.
4290
4291	The randomness in the algorithm comes from random sampling in numpy.
4292	To obtain the same results you can call numpy.random.seed with the
4293	same seed immediately before calling scipy.optimize.anneal.
4294
4295	We give a brief description of how the three temperature schedules
4296	generate new points and vary their temperature. Temperatures are
4297	only updated with iterations in the outer loop. The inner loop is
4298	over range(dwell), and new points are generated for
4299	every iteration in the inner loop. (Though whether the proposed
4300	new points are accepted is probabilistic.)
4301
4302	For readability, let d denote the dimension of the inputs to func.
4303	Also, let x_old denote the previous state, and k denote the
4304	iteration number of the outer loop. All other variables not
4305	defined below are input variables to scipy.optimize.anneal itself.
4306
4307	In the 'fast' schedule the updates are ::
4308
4309	u ~ Uniform(0, 1, size=d)
4310	y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4311	xc = y * (upper - lower)
4312	x_new = x_old + xc
4313
4314	T_new = T0 * exp(-c * k**quench)
4315
4316	"""
4317
4318	''' Scipy license:
4319	Copyright (c) 2001, 2002 Enthought, Inc.
4320	All rights reserved.
4321
4322	Copyright (c) 2003-2016 SciPy Developers.
4323	All rights reserved.
4324
4325	Redistribution and use in source and binary forms, with or without
4326	modification, are permitted provided that the following conditions are met:
4327
4328	a. Redistributions of source code must retain the above copyright notice,
4329	this list of conditions and the following disclaimer.
4330	b. Redistributions in binary form must reproduce the above copyright
4331	notice, this list of conditions and the following disclaimer in the
4332	documentation and/or other materials provided with the distribution.
4333	c. Neither the name of Enthought nor the names of the SciPy Developers
4334	may be used to endorse or promote products derived from this software
4335	without specific prior written permission.
4336	'''
4337	x0 = asarray(x0)
4338	lower = asarray(lower)
4339	upper = asarray(upper)
4340
4341	schedule = eval(schedule+'_sa()')
4342	# initialize the schedule
4343	schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4344	c=c, quench=quench, dwell=dwell, slope=slope)
4345
4346	current_state, last_state, best_state = _state(), _state(), _state()
4347	if ranStart:
4348	schedule.set_range(x0,ranRange)
4349	if T0 is None:
4350	x0 = schedule.getstart_temp(best_state)
4351	else:
4352	x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4353	best_state.x = None
4354	best_state.cost = numpy.Inf
4355
4356	last_state.x = asarray(x0).copy()
4357	fval = func(x0,*args)
4358	schedule.feval += 1
4359	last_state.cost = fval
4360	if last_state.cost < best_state.cost:
4361	best_state.cost = fval
4362	best_state.x = asarray(x0).copy()
4363	schedule.T = schedule.T0
4364	fqueue = [100, 300, 500, 700]
4365	iters = 1
4366	keepGoing = True
4367	bestn = 0
4368	while keepGoing:
4369	retval = 0
4370	for n in range(dwell):
4371	current_state.x = schedule.update_guess(last_state.x)
4372	current_state.cost = func(current_state.x,*args)
4373	schedule.feval += 1
4374
4375	dE = current_state.cost - last_state.cost
4376	if schedule.accept_test(dE):
4377	last_state.x = current_state.x.copy()
4378	last_state.cost = current_state.cost
4379	if last_state.cost < best_state.cost:
4380	best_state.x = last_state.x.copy()
4381	best_state.cost = last_state.cost
4382	bestn = n
4383	if best_state.cost < 1.0 and autoRan:
4384	schedule.set_range(x0,best_state.cost/2.)
4385	if dlg:
4386	GoOn = dlg.Update(min(100.,best_state.cost*100),
4387	newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4388	'Best trial:',bestn, \
4389	'MC/SA Residual:',best_state.cost*100,'%', \
4390	))[0]
4391	if not GoOn:
4392	best_state.x = last_state.x.copy()
4393	best_state.cost = last_state.cost
4394	retval = 5
4395	schedule.update_temp()
4396	iters += 1
4397	# Stopping conditions
4398	# 0) last saved values of f from each cooling step
4399	# are all very similar (effectively cooled)
4400	# 1) Tf is set and we are below it
4401	# 2) maxeval is set and we are past it
4402	# 3) maxiter is set and we are past it
4403	# 4) maxaccept is set and we are past it
4404	# 5) user canceled run via progress bar
4405
4406	fqueue.append(squeeze(last_state.cost))
4407	fqueue.pop(0)
4408	af = asarray(fqueue)*1.0
4409	if retval == 5:
4410	print (' User terminated run; incomplete MC/SA')
4411	keepGoing = False
4412	break
4413	if all(abs((af-af[0])/af[0]) < feps):
4414	retval = 0
4415	if abs(af[-1]-best_state.cost) > feps*10:
4416	retval = 5
4417	print (" Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1))
4418	break
4419	if (Tf is not None) and (schedule.T < Tf):
4420	# print ' Minimum T reached in %d steps'%(iters-1)
4421	retval = 1
4422	break
4423	if (maxeval is not None) and (schedule.feval > maxeval):
4424	retval = 2
4425	break
4426	if (iters > maxiter):
4427	print (" Warning: Maximum number of iterations exceeded.")
4428	retval = 3
4429	break
4430	if (maxaccept is not None) and (schedule.accepted > maxaccept):
4431	retval = 4
4432	break
4433
4434	return best_state.x, best_state.cost, schedule.T, \
4435	schedule.feval, iters, schedule.accepted, retval
4436
4437	def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4438	outlist = []
4439	timelist = []
4440	nsflist = []
4441	random.seed(int(time.time())%100000+nprocess) #make sure each process has a different random start
4442	for n in range(iCyc):
4443	result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False) #mcsa result,time,rcov
4444	outlist.append(result[0])
4445	timelist.append(result[1])
4446	nsflist.append(result[2])
4447	print (' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1]))
4448	out_q.put(outlist)
4449	out_t.put(timelist)
4450	out_n.put(nsflist)
4451
4452	def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4453	import multiprocessing as mp
4454
4455	out_q = mp.Queue()
4456	out_t = mp.Queue()
4457	out_n = mp.Queue()
4458	procs = []
4459	totsftime = 0.
4460	totnsf = 0
4461	iCyc = np.zeros(nprocs)
4462	for i in range(nCyc):
4463	iCyc[i%nprocs] += 1
4464	for i in range(nprocs):
4465	p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4466	procs.append(p)
4467	p.start()
4468	resultlist = []
4469	for i in range(nprocs):
4470	resultlist += out_q.get()
4471	totsftime += np.sum(out_t.get())
4472	totnsf += np.sum(out_n.get())
4473	for p in procs:
4474	p.join()
4475	return resultlist,totsftime,totnsf
4476
4477	def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4478	'''default doc string
4479
4480	:param type name: description
4481
4482	:returns: type name: description
4483	'''
4484
4485	class RandomDisplacementBounds(object):
4486	"""random displacement with bounds"""
4487	def __init__(self, xmin, xmax, stepsize=0.5):
4488	self.xmin = xmin
4489	self.xmax = xmax
4490	self.stepsize = stepsize
4491
4492	def __call__(self, x):
4493	"""take a random step but ensure the new position is within the bounds"""
4494	while True:
4495	# this could be done in a much more clever way, but it will work for example purposes
4496	steps = self.xmax-self.xmin
4497	xnew = x + np.random.uniform(-self.stepsizesteps, self.stepsizesteps, np.shape(x))
4498	if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4499	break
4500	return xnew
4501
4502	global tsum,nsum
4503	tsum = 0.
4504	nsum = 0
4505
4506	def getMDparms(item,pfx,parmDict,varyList):
4507	parmDict[pfx+'MDaxis'] = item['axis']
4508	parmDict[pfx+'MDval'] = item['Coef'][0]
4509	if item['Coef'][1]:
4510	varyList += [pfx+'MDval',]
4511	limits = item['Coef'][2]
4512	lower.append(limits[0])
4513	upper.append(limits[1])
4514
4515	def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4516	parmDict[pfx+'Atype'] = item['atType']
4517	aTypes \|= set([item['atType'],])
4518	pstr = ['Ax','Ay','Az']
4519	XYZ = [0,0,0]
4520	for i in range(3):
4521	name = pfx+pstr[i]
4522	parmDict[name] = item['Pos'][0][i]
4523	XYZ[i] = parmDict[name]
4524	if item['Pos'][1][i]:
4525	varyList += [name,]
4526	limits = item['Pos'][2][i]
4527	lower.append(limits[0])
4528	upper.append(limits[1])
4529	parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4530
4531	def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4532	parmDict[mfx+'MolCent'] = item['MolCent']
4533	parmDict[mfx+'RBId'] = item['RBId']
4534	pstr = ['Px','Py','Pz']
4535	ostr = ['Qa','Qi','Qj','Qk'] #angle,vector not quaternion
4536	for i in range(3):
4537	name = mfx+pstr[i]
4538	parmDict[name] = item['Pos'][0][i]
4539	if item['Pos'][1][i]:
4540	varyList += [name,]
4541	limits = item['Pos'][2][i]
4542	lower.append(limits[0])
4543	upper.append(limits[1])
4544	AV = item['Ori'][0]
4545	A = AV[0]
4546	V = AV[1:]
4547	for i in range(4):
4548	name = mfx+ostr[i]
4549	if i:
4550	parmDict[name] = V[i-1]
4551	else:
4552	parmDict[name] = A
4553	if item['Ovar'] == 'AV':
4554	varyList += [name,]
4555	limits = item['Ori'][2][i]
4556	lower.append(limits[0])
4557	upper.append(limits[1])
4558	elif item['Ovar'] == 'A' and not i:
4559	varyList += [name,]
4560	limits = item['Ori'][2][i]
4561	lower.append(limits[0])
4562	upper.append(limits[1])
4563	if 'Tor' in item: #'Tor' not there for 'Vector' RBs
4564	for i in range(len(item['Tor'][0])):
4565	name = mfx+'Tor'+str(i)
4566	parmDict[name] = item['Tor'][0][i]
4567	if item['Tor'][1][i]:
4568	varyList += [name,]
4569	�