Context Navigation

source: trunk/GSASIImath.py @ 3865

Last change on this file since 3865 was 3865, checked in by vondreele, 4 years ago
a new attempt at incommensurate mag. str. fctrs.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 188.6 KB

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