Context Navigation

source: trunk/GSASIImath.py @ 2025

Last change on this file since 2025 was 2025, checked in by vondreele, 8 years ago

Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 149.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: