Context Navigation

source: trunk/GSASIIstruct.py @ 737

Last change on this file since 737 was 735, checked in by vondreele, 13 years ago
back to 701 (modified)
Property svn:keywords set to `Date Author Revision URL Id`
File size: 166.8 KB

Line
1	# -- coding: utf-8 --
2	#GSASIIstructure - structure computation routines
3	########### SVN repository information ###################
4	# $Date: 2012-08-23 22:22:35 +0000 (Thu, 23 Aug 2012) $
5	# $Author: toby $
6	# $Revision: 735 $
7	# $URL: trunk/GSASIIstruct.py $
8	# $Id: GSASIIstruct.py 735 2012-08-23 22:22:35Z toby $
9	########### SVN repository information ###################
10	import sys
11	import os
12	import os.path as ospath
13	import time
14	import math
15	import cPickle
16	import numpy as np
17	import numpy.linalg as nl
18	import scipy.optimize as so
19	import GSASIIpath
20	GSASIIpath.SetVersionNumber("$Revision: 735 $")
21	import GSASIIElem as G2el
22	import GSASIIlattice as G2lat
23	import GSASIIspc as G2spc
24	import GSASIIpwd as G2pwd
25	import GSASIImapvars as G2mv
26	import GSASIImath as G2mth
27
28	sind = lambda x: np.sin(x*np.pi/180.)
29	cosd = lambda x: np.cos(x*np.pi/180.)
30	tand = lambda x: np.tan(x*np.pi/180.)
31	asind = lambda x: 180.*np.arcsin(x)/np.pi
32	acosd = lambda x: 180.*np.arccos(x)/np.pi
33	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
34
35	ateln2 = 8.0*math.log(2.0)
36	DEBUG = False
37
38	def GetControls(GPXfile):
39	''' Returns dictionary of control items found in GSASII gpx file
40	input:
41	GPXfile = .gpx full file name
42	return:
43	Controls = dictionary of control items
44	'''
45	Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
46	'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
47	fl = open(GPXfile,'rb')
48	while True:
49	try:
50	data = cPickle.load(fl)
51	except EOFError:
52	break
53	datum = data[0]
54	if datum[0] == 'Controls':
55	Controls.update(datum[1])
56	fl.close()
57	return Controls
58
59	def GetConstraints(GPXfile):
60	'''Read the constraints from the GPX file and interpret them
61	'''
62	constList = []
63	fl = open(GPXfile,'rb')
64	while True:
65	try:
66	data = cPickle.load(fl)
67	except EOFError:
68	break
69	datum = data[0]
70	if datum[0] == 'Constraints':
71	constDict = datum[1]
72	for item in constDict:
73	constList += constDict[item]
74	fl.close()
75	constDict,fixedList,ignored = ProcessConstraints(constList)
76	if ignored:
77	print ignored,'old-style Constraints were rejected'
78	return constDict,fixedList
79
80	def ProcessConstraints(constList):
81	"interpret constraints"
82	constDict = []
83	fixedList = []
84	ignored = 0
85	for item in constList:
86	if item[-1] == 'h':
87	# process a hold
88	fixedList.append('0')
89	constDict.append({item[0][1]:0.0})
90	elif item[-1] == 'f':
91	# process a new variable
92	fixedList.append(None)
93	constDict.append({})
94	for term in item[:-3]:
95	constDict[-1][term[1]] = term[0]
96	#constFlag[-1] = ['Vary']
97	elif item[-1] == 'c':
98	# process a contraint relationship
99	fixedList.append(str(item[-3]))
100	constDict.append({})
101	for term in item[:-3]:
102	constDict[-1][term[1]] = term[0]
103	#constFlag[-1] = ['VaryFree']
104	elif item[-1] == 'e':
105	# process an equivalence
106	firstmult = None
107	eqlist = []
108	for term in item[:-3]:
109	if term[0] == 0: term[0] = 1.0
110	if firstmult is None:
111	firstmult,firstvar = term
112	else:
113	eqlist.append([term[1],firstmult/term[0]])
114	G2mv.StoreEquivalence(firstvar,eqlist)
115	else:
116	ignored += 1
117	return constDict,fixedList,ignored
118
119	def CheckConstraints(GPXfile):
120	'''Load constraints and related info and return any error or warning messages'''
121	# init constraints
122	G2mv.InitVars()
123	# get variables
124	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
125	if not Phases:
126	return 'Error: No Phases!',''
127	if not Histograms:
128	return 'Error: no diffraction data',''
129	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,Print=False)
130	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
131	histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
132	varyList = phaseVary+hapVary+histVary
133	constrDict,fixedList = GetConstraints(GPXfile)
134	return G2mv.CheckConstraints(varyList,constrDict,fixedList)
135
136	def GetRestraints(GPXfile):
137	'''Read the restraints from the GPX file
138	'''
139	fl = open(GPXfile,'rb')
140	while True:
141	try:
142	data = cPickle.load(fl)
143	except EOFError:
144	break
145	datum = data[0]
146	if datum[0] == 'Restraints':
147	restraintDict = datum[1]
148	fl.close()
149	return restraintDict
150
151	def GetPhaseNames(GPXfile):
152	''' Returns a list of phase names found under 'Phases' in GSASII gpx file
153	input:
154	GPXfile = gpx full file name
155	return:
156	PhaseNames = list of phase names
157	'''
158	fl = open(GPXfile,'rb')
159	PhaseNames = []
160	while True:
161	try:
162	data = cPickle.load(fl)
163	except EOFError:
164	break
165	datum = data[0]
166	if 'Phases' == datum[0]:
167	for datus in data[1:]:
168	PhaseNames.append(datus[0])
169	fl.close()
170	return PhaseNames
171
172	def GetAllPhaseData(GPXfile,PhaseName):
173	''' Returns the entire dictionary for PhaseName from GSASII gpx file
174	input:
175	GPXfile = gpx full file name
176	PhaseName = phase name
177	return:
178	phase dictionary
179	'''
180	fl = open(GPXfile,'rb')
181	General = {}
182	Atoms = []
183	while True:
184	try:
185	data = cPickle.load(fl)
186	except EOFError:
187	break
188	datum = data[0]
189	if 'Phases' == datum[0]:
190	for datus in data[1:]:
191	if datus[0] == PhaseName:
192	break
193	fl.close()
194	return datus[1]
195
196	def GetHistograms(GPXfile,hNames):
197	""" Returns a dictionary of histograms found in GSASII gpx file
198	input:
199	GPXfile = .gpx full file name
200	hNames = list of histogram names
201	return:
202	Histograms = dictionary of histograms (types = PWDR & HKLF)
203	"""
204	fl = open(GPXfile,'rb')
205	Histograms = {}
206	while True:
207	try:
208	data = cPickle.load(fl)
209	except EOFError:
210	break
211	datum = data[0]
212	hist = datum[0]
213	if hist in hNames:
214	if 'PWDR' in hist[:4]:
215	PWDRdata = {}
216	try:
217	PWDRdata.update(datum[1][0]) #weight factor
218	except ValueError:
219	PWDRdata['wtFactor'] = 1.0 #patch
220	PWDRdata['Data'] = datum[1][1] #powder data arrays
221	PWDRdata[data[2][0]] = data[2][1] #Limits
222	PWDRdata[data[3][0]] = data[3][1] #Background
223	PWDRdata[data[4][0]] = data[4][1] #Instrument parameters
224	PWDRdata[data[5][0]] = data[5][1] #Sample parameters
225	try:
226	PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing
227	except IndexError:
228	PWDRdata['Reflection Lists'] = {}
229
230	Histograms[hist] = PWDRdata
231	elif 'HKLF' in hist[:4]:
232	HKLFdata = {}
233	try:
234	HKLFdata.update(datum[1][0]) #weight factor
235	except ValueError:
236	HKLFdata['wtFactor'] = 1.0 #patch
237	HKLFdata['Data'] = datum[1][1]
238	HKLFdata[data[1][0]] = data[1][1] #Instrument parameters
239	HKLFdata['Reflection Lists'] = None
240	Histograms[hist] = HKLFdata
241	fl.close()
242	return Histograms
243
244	def GetHistogramNames(GPXfile,hType):
245	""" Returns a list of histogram names found in GSASII gpx file
246	input:
247	GPXfile = .gpx full file name
248	hType = list ['PWDR','HKLF']
249	return:
250	HistogramNames = list of histogram names (types = PWDR & HKLF)
251	"""
252	fl = open(GPXfile,'rb')
253	HistogramNames = []
254	while True:
255	try:
256	data = cPickle.load(fl)
257	except EOFError:
258	break
259	datum = data[0]
260	if datum[0][:4] in hType:
261	HistogramNames.append(datum[0])
262	fl.close()
263	return HistogramNames
264
265	def GetUsedHistogramsAndPhases(GPXfile):
266	''' Returns all histograms that are found in any phase
267	and any phase that uses a histogram
268	input:
269	GPXfile = .gpx full file name
270	return:
271	Histograms = dictionary of histograms as {name:data,...}
272	Phases = dictionary of phases that use histograms
273	'''
274	phaseNames = GetPhaseNames(GPXfile)
275	histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
276	allHistograms = GetHistograms(GPXfile,histoList)
277	phaseData = {}
278	for name in phaseNames:
279	phaseData[name] = GetAllPhaseData(GPXfile,name)
280	Histograms = {}
281	Phases = {}
282	for phase in phaseData:
283	Phase = phaseData[phase]
284	if Phase['Histograms']:
285	if phase not in Phases:
286	pId = phaseNames.index(phase)
287	Phase['pId'] = pId
288	Phases[phase] = Phase
289	for hist in Phase['Histograms']:
290	if 'Use' not in Phase['Histograms'][hist]: #patch
291	Phase['Histograms'][hist]['Use'] = True
292	if hist not in Histograms and Phase['Histograms'][hist]['Use']:
293	Histograms[hist] = allHistograms[hist]
294	#future restraint, etc. histograms here
295	hId = histoList.index(hist)
296	Histograms[hist]['hId'] = hId
297	return Histograms,Phases
298
299	def getBackupName(GPXfile,makeBack):
300	GPXpath,GPXname = ospath.split(GPXfile)
301	if GPXpath == '': GPXpath = '.'
302	Name = ospath.splitext(GPXname)[0]
303	files = os.listdir(GPXpath)
304	last = 0
305	for name in files:
306	name = name.split('.')
307	if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
308	if makeBack:
309	last = max(last,int(name[1].strip('bak'))+1)
310	else:
311	last = max(last,int(name[1].strip('bak')))
312	GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
313	return GPXback
314
315	def GPXBackup(GPXfile,makeBack=True):
316	import distutils.file_util as dfu
317	GPXback = getBackupName(GPXfile,makeBack)
318	dfu.copy_file(GPXfile,GPXback)
319	return GPXback
320
321	def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
322	''' Updates gpxfile from all histograms that are found in any phase
323	and any phase that used a histogram
324	input:
325	GPXfile = .gpx full file name
326	Histograms = dictionary of histograms as {name:data,...}
327	Phases = dictionary of phases that use histograms
328	CovData = dictionary of refined variables, varyList, & covariance matrix
329	makeBack = True if new backup of .gpx file is to be made; else use the last one made
330	'''
331
332	GPXback = GPXBackup(GPXfile,makeBack)
333	print 'Read from file:',GPXback
334	print 'Save to file :',GPXfile
335	infile = open(GPXback,'rb')
336	outfile = open(GPXfile,'wb')
337	while True:
338	try:
339	data = cPickle.load(infile)
340	except EOFError:
341	break
342	datum = data[0]
343	# print 'read: ',datum[0]
344	if datum[0] == 'Phases':
345	for iphase in range(len(data)):
346	if data[iphase][0] in Phases:
347	phaseName = data[iphase][0]
348	data[iphase][1].update(Phases[phaseName])
349	elif datum[0] == 'Covariance':
350	data[0][1] = CovData
351	try:
352	histogram = Histograms[datum[0]]
353	# print 'found ',datum[0]
354	data[0][1][1] = histogram['Data']
355	for datus in data[1:]:
356	# print ' read: ',datus[0]
357	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
358	datus[1] = histogram[datus[0]]
359	except KeyError:
360	pass
361
362	cPickle.dump(data,outfile,1)
363	infile.close()
364	outfile.close()
365	print 'GPX file save successful'
366
367	def SetSeqResult(GPXfile,Histograms,SeqResult):
368	GPXback = GPXBackup(GPXfile)
369	print 'Read from file:',GPXback
370	print 'Save to file :',GPXfile
371	infile = open(GPXback,'rb')
372	outfile = open(GPXfile,'wb')
373	while True:
374	try:
375	data = cPickle.load(infile)
376	except EOFError:
377	break
378	datum = data[0]
379	if datum[0] == 'Sequental results':
380	data[0][1] = SeqResult
381	try:
382	histogram = Histograms[datum[0]]
383	data[0][1][1] = histogram['Data']
384	for datus in data[1:]:
385	if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
386	datus[1] = histogram[datus[0]]
387	except KeyError:
388	pass
389
390	cPickle.dump(data,outfile,1)
391	infile.close()
392	outfile.close()
393	print 'GPX file save successful'
394
395	def ShowBanner(pFile=None):
396	print >>pFile,80''
397	print >>pFile,' General Structure Analysis System-II Crystal Structure Refinement'
398	print >>pFile,' by Robert B. Von Dreele & Brian H. Toby'
399	print >>pFile,' Argonne National Laboratory(C), 2010'
400	print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,'
401	print >>pFile,' as Operator of Argonne National Laboratory.'
402	print >>pFile,80'','\n'
403
404	def ShowControls(Controls,pFile=None):
405	print >>pFile,' Least squares controls:'
406	print >>pFile,' Refinement type: ',Controls['deriv type']
407	if 'Hessian' in Controls['deriv type']:
408	print >>pFile,' Maximum number of cycles:',Controls['max cyc']
409	else:
410	print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
411	print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
412
413	def GetFFtable(General):
414	''' returns a dictionary of form factor data for atom types found in General
415	input:
416	General = dictionary of phase info.; includes AtomTypes
417	return:
418	FFtable = dictionary of form factor data; key is atom type
419	'''
420	atomTypes = General['AtomTypes']
421	FFtable = {}
422	for El in atomTypes:
423	FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
424	for item in FFs:
425	if item['Symbol'] == El.upper():
426	FFtable[El] = item
427	return FFtable
428
429	def GetBLtable(General):
430	''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
431	input:
432	General = dictionary of phase info.; includes AtomTypes & Isotopes
433	return:
434	BLtable = dictionary of scattering length data; key is atom type
435	'''
436	atomTypes = General['AtomTypes']
437	BLtable = {}
438	isotopes = General['Isotopes']
439	isotope = General['Isotope']
440	for El in atomTypes:
441	BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
442	return BLtable
443
444	def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
445	# if SGLaue in ['-1','2/m','mmm']:
446	# return #no Pawley symmetry required constraints
447	eqvDict = {}
448	for i,varyI in enumerate(pawleyVary):
449	eqvDict[varyI] = []
450	refI = int(varyI.split(':')[-1])
451	ih,ik,il = PawleyRef[refI][:3]
452	dspI = PawleyRef[refI][4]
453	for varyJ in pawleyVary[i+1:]:
454	refJ = int(varyJ.split(':')[-1])
455	jh,jk,jl = PawleyRef[refJ][:3]
456	dspJ = PawleyRef[refJ][4]
457	if SGLaue in ['4/m','4/mmm']:
458	isum = ih2+ik2
459	jsum = jh2+jk2
460	if abs(il) == abs(jl) and isum == jsum:
461	eqvDict[varyI].append(varyJ)
462	elif SGLaue in ['3R','3mR']:
463	isum = ih2+ik2+il**2
464	jsum = jh2+jk2jl*2
465	isum2 = ihik+ihil+ik*il
466	jsum2 = jhjk+jhjl+jk*jl
467	if isum == jsum and isum2 == jsum2:
468	eqvDict[varyI].append(varyJ)
469	elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
470	isum = ih2+ik2+ih*ik
471	jsum = jh2+jk2+jh*jk
472	if abs(il) == abs(jl) and isum == jsum:
473	eqvDict[varyI].append(varyJ)
474	elif SGLaue in ['m3','m3m']:
475	isum = ih2+ik2+il**2
476	jsum = jh2+jk2+jl**2
477	if isum == jsum:
478	eqvDict[varyI].append(varyJ)
479	elif abs(dspI-dspJ)/dspI < 1.e-4:
480	eqvDict[varyI].append(varyJ)
481	for item in pawleyVary:
482	if eqvDict[item]:
483	for item2 in pawleyVary:
484	if item2 in eqvDict[item]:
485	eqvDict[item2] = []
486	G2mv.StoreEquivalence(item,eqvDict[item])
487
488	def cellVary(pfx,SGData):
489	if SGData['SGLaue'] in ['-1',]:
490	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
491	elif SGData['SGLaue'] in ['2/m',]:
492	if SGData['SGUniq'] == 'a':
493	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
494	elif SGData['SGUniq'] == 'b':
495	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
496	else:
497	return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
498	elif SGData['SGLaue'] in ['mmm',]:
499	return [pfx+'A0',pfx+'A1',pfx+'A2']
500	elif SGData['SGLaue'] in ['4/m','4/mmm']:
501	return [pfx+'A0',pfx+'A2']
502	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
503	return [pfx+'A0',pfx+'A2']
504	elif SGData['SGLaue'] in ['3R', '3mR']:
505	return [pfx+'A0',pfx+'A3']
506	elif SGData['SGLaue'] in ['m3m','m3']:
507	return [pfx+'A0',]
508
509	################################################################################
510	##### Phase data
511	################################################################################
512
513	def GetPhaseData(PhaseData,RestraintDict=None,Print=True,pFile=None):
514
515	def PrintFFtable(FFtable):
516	print >>pFile,'\n X-ray scattering factors:'
517	print >>pFile,' Symbol fa fb fc'
518	print >>pFile,99*'-'
519	for Ename in FFtable:
520	ffdata = FFtable[Ename]
521	fa = ffdata['fa']
522	fb = ffdata['fb']
523	print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \
524	(Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
525
526	def PrintBLtable(BLtable):
527	print >>pFile,'\n Neutron scattering factors:'
528	print >>pFile,' Symbol isotope mass b resonant terms'
529	print >>pFile,99*'-'
530	for Ename in BLtable:
531	bldata = BLtable[Ename]
532	isotope = bldata[0]
533	mass = bldata[1][0]
534	blen = bldata[1][1]
535	bres = []
536	if len(bldata[1]) > 2:
537	bres = bldata[1][2:]
538	line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
539	for item in bres:
540	line += '%10.5g'%(item)
541	print >>pFile,line
542
543	def PrintAtoms(General,Atoms):
544	print >>pFile,'\n Atoms:'
545	line = ' name type refine? x y z '+ \
546	' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23'
547	if General['Type'] == 'magnetic':
548	line += ' Mx My Mz'
549	elif General['Type'] == 'macromolecular':
550	line = ' res no residue chain '+line
551	print >>pFile,line
552	if General['Type'] == 'nuclear':
553	print >>pFile,135*'-'
554	for i,at in enumerate(Atoms):
555	line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
556	'%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
557	if at[9] == 'I':
558	line += '%8.4f'%(at[10])+48*' '
559	else:
560	line += 8*' '
561	for j in range(6):
562	line += '%8.4f'%(at[11+j])
563	print >>pFile,line
564
565	def PrintTexture(textureData):
566	topstr = '\n Spherical harmonics texture: Order:' + \
567	str(textureData['Order'])
568	if textureData['Order']:
569	print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
570	else:
571	print >>pFile,topstr
572	return
573	names = ['omega','chi','phi']
574	line = '\n'
575	for name in names:
576	line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
577	print >>pFile,line
578	print >>pFile,'\n Texture coefficients:'
579	ptlbls = ' names :'
580	ptstr = ' values:'
581	SHcoeff = textureData['SH Coeff'][1]
582	for item in SHcoeff:
583	ptlbls += '%12s'%(item)
584	ptstr += '%12.4f'%(SHcoeff[item])
585	print >>pFile,ptlbls
586	print >>pFile,ptstr
587
588	def PrintRestraints(phaseRest):
589	if phaseRest:
590	print >>pFile,'\n Restraints:'
591	names = ['Bonds','Angles','Planes','Volumes']
592	for i,name in enumerate(['Bond','Angle','Plane','Chiral']):
593	itemRest = phaseRest[name]
594	if itemRest[names[i]]:
595	print >>pFile,'\n %30s %10.3f'%(name+' restraint weight factor',itemRest['wtFactor'])
596	print >>pFile,' atoms(symOp), calc, obs, sig: '
597	for item in phaseRest[name][names[i]]:
598	text = ' '
599	for a,at in enumerate(item[0]):
600	text += at+'+('+item[1][a]+') '
601	if (a+1)%5 == 0:
602	text += '\n '
603	print >>pFile,text,' %.3f %.3f %.3f'%(item[3],item[4],item[5])
604
605	if Print:print >>pFile,' Phases:'
606	phaseVary = []
607	phaseDict = {}
608	phaseConstr = {}
609	pawleyLookup = {}
610	FFtables = {} #scattering factors - xrays
611	BLtables = {} # neutrons
612	Natoms = {}
613	AtMults = {}
614	AtIA = {}
615	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
616	SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
617	atomIndx = {}
618	for name in PhaseData:
619	General = PhaseData[name]['General']
620	pId = PhaseData[name]['pId']
621	pfx = str(pId)+'::'
622	FFtable = GetFFtable(General)
623	BLtable = GetBLtable(General)
624	FFtables.update(FFtable)
625	BLtables.update(BLtable)
626	Atoms = PhaseData[name]['Atoms']
627	try:
628	PawleyRef = PhaseData[name]['Pawley ref']
629	except KeyError:
630	PawleyRef = []
631	SGData = General['SGData']
632	SGtext = G2spc.SGPrint(SGData)
633	cell = General['Cell']
634	A = G2lat.cell2A(cell[1:7])
635	phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
636	if cell[0]:
637	phaseVary += cellVary(pfx,SGData)
638	Natoms[pfx] = 0
639	if Atoms and not General.get('doPawley'):
640	if General['Type'] == 'nuclear':
641	Natoms[pfx] = len(Atoms)
642	for i,at in enumerate(Atoms):
643	atomIndx[at[-1]] = [pfx,i] #lookup table for restraints
644	phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
645	pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
646	pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z
647	pfx+'AI/A:'+str(i):at[9],})
648	if at[9] == 'I':
649	phaseDict[pfx+'AUiso:'+str(i)] = at[10]
650	else:
651	phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
652	pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
653	if 'F' in at[2]:
654	phaseVary.append(pfx+'Afrac:'+str(i))
655	if 'X' in at[2]:
656	xId,xCoef = G2spc.GetCSxinel(at[7])
657	names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
658	equivs = [[],[],[]]
659	for j in range(3):
660	if xId[j] > 0:
661	phaseVary.append(names[j])
662	equivs[xId[j]-1].append([names[j],xCoef[j]])
663	for equiv in equivs:
664	if len(equiv) > 1:
665	name = equiv[0][0]
666	for eqv in equiv[1:]:
667	G2mv.StoreEquivalence(name,(eqv,))
668	if 'U' in at[2]:
669	if at[9] == 'I':
670	phaseVary.append(pfx+'AUiso:'+str(i))
671	else:
672	uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
673	names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
674	pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
675	equivs = [[],[],[],[],[],[]]
676	for j in range(6):
677	if uId[j] > 0:
678	phaseVary.append(names[j])
679	equivs[uId[j]-1].append([names[j],uCoef[j]])
680	for equiv in equivs:
681	if len(equiv) > 1:
682	name = equiv[0][0]
683	for eqv in equiv[1:]:
684	G2mv.StoreEquivalence(name,(eqv,))
685	# elif General['Type'] == 'magnetic':
686	# elif General['Type'] == 'macromolecular':
687	textureData = General['SH Texture']
688	if textureData['Order']:
689	phaseDict[pfx+'SHorder'] = textureData['Order']
690	phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
691	for name in ['omega','chi','phi']:
692	phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
693	if textureData['Sample '+name][0]:
694	phaseVary.append(pfx+'SH '+name)
695	for name in textureData['SH Coeff'][1]:
696	phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
697	if textureData['SH Coeff'][0]:
698	phaseVary.append(pfx+name)
699
700	if Print:
701	print >>pFile,'\n Phase name: ',General['Name']
702	print >>pFile,135*'-'
703	PrintFFtable(FFtable)
704	PrintBLtable(BLtable)
705	print >>pFile,''
706	for line in SGtext: print >>pFile,line
707	PrintAtoms(General,Atoms)
708	print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
709	' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
710	'%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
711	PrintTexture(textureData)
712	if name in RestraintDict:
713	PrintRestraints(RestraintDict[name])
714
715	elif PawleyRef:
716	pawleyVary = []
717	for i,refl in enumerate(PawleyRef):
718	phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
719	pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
720	if refl[5]:
721	pawleyVary.append(pfx+'PWLref:'+str(i))
722	GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence
723	phaseVary += pawleyVary
724
725	return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
726
727	def cellFill(pfx,SGData,parmDict,sigDict):
728	if SGData['SGLaue'] in ['-1',]:
729	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
730	parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
731	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
732	sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
733	elif SGData['SGLaue'] in ['2/m',]:
734	if SGData['SGUniq'] == 'a':
735	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
736	parmDict[pfx+'A3'],0,0]
737	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
738	sigDict[pfx+'A3'],0,0]
739	elif SGData['SGUniq'] == 'b':
740	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
741	0,parmDict[pfx+'A4'],0]
742	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
743	0,sigDict[pfx+'A4'],0]
744	else:
745	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
746	0,0,parmDict[pfx+'A5']]
747	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
748	0,0,sigDict[pfx+'A5']]
749	elif SGData['SGLaue'] in ['mmm',]:
750	A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
751	sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
752	elif SGData['SGLaue'] in ['4/m','4/mmm']:
753	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
754	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
755	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
756	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
757	parmDict[pfx+'A0'],0,0]
758	sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
759	elif SGData['SGLaue'] in ['3R', '3mR']:
760	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
761	parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
762	sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
763	elif SGData['SGLaue'] in ['m3m','m3']:
764	A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
765	sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
766	return A,sigA
767
768	def getCellEsd(pfx,SGData,A,covData):
769	dpr = 180./np.pi
770	rVsq = G2lat.calc_rVsq(A)
771	G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors
772	cell = np.array(G2lat.Gmat2cell(g)) #real cell
773	cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
774	scos = cosd(cellst[3:6])
775	ssin = sind(cellst[3:6])
776	scot = scos/ssin
777	rcos = cosd(cell[3:6])
778	rsin = sind(cell[3:6])
779	rcot = rcos/rsin
780	RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
781	varyList = covData['varyList']
782	covMatrix = covData['covMatrix']
783	vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
784	Ax = np.array(A)
785	Ax[3:] /= 2.
786	drVdA = np.array([Ax[1]Ax[2]-Ax[5]2,Ax[0]Ax[2]-Ax[4]*2,Ax[0]Ax[1]-Ax[3]**2,
787	Ax[4]Ax[5]-Ax[2]Ax[3],Ax[3]Ax[5]-Ax[1]Ax[4],Ax[3]Ax[4]-Ax[0]Ax[5]])
788	srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
789	Vol = 1/np.sqrt(rVsq)
790	sigVol = Vol*3np.sqrt(srcvlsq)/2.
791	R123 = Ax[0]Ax[1]Ax[2]
792	dsasdg = np.zeros((3,6))
793	dadg = np.zeros((6,6))
794	for i0 in range(3): #0 1 2
795	i1 = (i0+1)%3 #1 2 0
796	i2 = (i1+1)%3 #2 0 1
797	i3 = 5-i2 #3 5 4
798	i4 = 5-i1 #4 3 5
799	i5 = 5-i0 #5 4 3
800	dsasdg[i0][i1] = 0.5scot[i0]scos[i0]/Ax[i1]
801	dsasdg[i0][i2] = 0.5scot[i0]scos[i0]/Ax[i2]
802	dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
803	denmsq = Ax[i0](R123-Ax[i1]Ax[i4]*2-Ax[i2]Ax[i3]*2+(Ax[i4]Ax[i3])**2)
804	denom = np.sqrt(denmsq)
805	dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq(R123-0.5Ax[i1]Ax[i4]2-0.5Ax[i2]Ax[i3]*2)
806	dadg[i5][i1] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i2]-Ax[i0]Ax[i4]*2)
807	dadg[i5][i2] = -0.5rcos[i0]/denmsq(Ax[i0]*2Ax[i1]-Ax[i0]Ax[i3]*2)
808	dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i2]Ax[i3]-Ax[i3]Ax[i4]**2)
809	dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq(Ax[i0]Ax[i1]Ax[i4]-Ax[i3]2Ax[i4])
810	dadg[i5][i5] = -Ax[i0]/denom
811	for i0 in range(3):
812	i1 = (i0+1)%3
813	i2 = (i1+1)%3
814	i3 = 5-i2
815	for ij in range(6):
816	dadg[i0][ij] = cell[i0](rcot[i2]dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
817	if ij == i0:
818	dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
819	dadg[i3][ij] = -dadg[i3][ij]rsin[2-i0]dpr
820	sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
821	var = np.diag(sigMat)
822	CS = np.where(var>0.,np.sqrt(var),0.)
823	cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order
824	return cellSig
825
826	def SetPhaseData(parmDict,sigDict,Phases,covData,pFile=None):
827
828	def PrintAtomsAndSig(General,Atoms,atomsSig):
829	print >>pFile,'\n Atoms:'
830	line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23'
831	if General['Type'] == 'magnetic':
832	line += ' Mx My Mz'
833	elif General['Type'] == 'macromolecular':
834	line = ' res no residue chain '+line
835	print >>pFile,line
836	if General['Type'] == 'nuclear':
837	print >>pFile,135*'-'
838	fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
839	11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
840	noFXsig = {3:[10' ','%10s'],4:[10' ','%10s'],5:[10' ','%10s'],6:[8' ','%8s']}
841	for atyp in General['AtomTypes']: #zero composition
842	General['NoAtoms'][atyp] = 0.
843	for i,at in enumerate(Atoms):
844	General['NoAtoms'][at[1]] += at[6]*float(at[8]) #new composition
845	name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
846	valstr = ' values:'
847	sigstr = ' sig :'
848	for ind in [3,4,5,6]:
849	sigind = str(i)+':'+str(ind)
850	valstr += fmt[ind]%(at[ind])
851	if sigind in atomsSig:
852	sigstr += fmt[ind]%(atomsSig[sigind])
853	else:
854	sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
855	if at[9] == 'I':
856	valstr += fmt[10]%(at[10])
857	if str(i)+':10' in atomsSig:
858	sigstr += fmt[10]%(atomsSig[str(i)+':10'])
859	else:
860	sigstr += 8*' '
861	else:
862	valstr += 8*' '
863	sigstr += 8*' '
864	for ind in [11,12,13,14,15,16]:
865	sigind = str(i)+':'+str(ind)
866	valstr += fmt[ind]%(at[ind])
867	if sigind in atomsSig:
868	sigstr += fmt[ind]%(atomsSig[sigind])
869	else:
870	sigstr += 8*' '
871	print >>pFile,name
872	print >>pFile,valstr
873	print >>pFile,sigstr
874
875	def PrintSHtextureAndSig(textureData,SHtextureSig):
876	print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
877	names = ['omega','chi','phi']
878	namstr = ' names :'
879	ptstr = ' values:'
880	sigstr = ' esds :'
881	for name in names:
882	namstr += '%12s'%(name)
883	ptstr += '%12.3f'%(textureData['Sample '+name][1])
884	if 'Sample '+name in SHtextureSig:
885	sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
886	else:
887	sigstr += 12*' '
888	print >>pFile,namstr
889	print >>pFile,ptstr
890	print >>pFile,sigstr
891	print >>pFile,'\n Texture coefficients:'
892	namstr = ' names :'
893	ptstr = ' values:'
894	sigstr = ' esds :'
895	SHcoeff = textureData['SH Coeff'][1]
896	for name in SHcoeff:
897	namstr += '%12s'%(name)
898	ptstr += '%12.3f'%(SHcoeff[name])
899	if name in SHtextureSig:
900	sigstr += '%12.3f'%(SHtextureSig[name])
901	else:
902	sigstr += 12*' '
903	print >>pFile,namstr
904	print >>pFile,ptstr
905	print >>pFile,sigstr
906
907
908	print >>pFile,'\n Phases:'
909	for phase in Phases:
910	print >>pFile,' Result for phase: ',phase
911	Phase = Phases[phase]
912	General = Phase['General']
913	SGData = General['SGData']
914	Atoms = Phase['Atoms']
915	cell = General['Cell']
916	pId = Phase['pId']
917	pfx = str(pId)+'::'
918	if cell[0]:
919	A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
920	cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol
921	print >>pFile,' Reciprocal metric tensor: '
922	ptfmt = "%15.9f"
923	names = ['A11','A22','A33','A12','A13','A23']
924	namstr = ' names :'
925	ptstr = ' values:'
926	sigstr = ' esds :'
927	for name,a,siga in zip(names,A,sigA):
928	namstr += '%15s'%(name)
929	ptstr += ptfmt%(a)
930	if siga:
931	sigstr += ptfmt%(siga)
932	else:
933	sigstr += 15*' '
934	print >>pFile,namstr
935	print >>pFile,ptstr
936	print >>pFile,sigstr
937	cell[1:7] = G2lat.A2cell(A)
938	cell[7] = G2lat.calc_V(A)
939	print >>pFile,' New unit cell:'
940	ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
941	names = ['a','b','c','alpha','beta','gamma','Volume']
942	namstr = ' names :'
943	ptstr = ' values:'
944	sigstr = ' esds :'
945	for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
946	namstr += '%12s'%(name)
947	ptstr += fmt%(a)
948	if siga:
949	sigstr += fmt%(siga)
950	else:
951	sigstr += 12*' '
952	print >>pFile,namstr
953	print >>pFile,ptstr
954	print >>pFile,sigstr
955
956	if Phase['General'].get('doPawley'):
957	pawleyRef = Phase['Pawley ref']
958	for i,refl in enumerate(pawleyRef):
959	key = pfx+'PWLref:'+str(i)
960	refl[6] = parmDict[key]
961	if key in sigDict:
962	refl[7] = sigDict[key]
963	else:
964	refl[7] = 0
965	else:
966	atomsSig = {}
967	if General['Type'] == 'nuclear':
968	for i,at in enumerate(Atoms):
969	names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
970	10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
971	14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
972	for ind in [3,4,5,6]:
973	at[ind] = parmDict[names[ind]]
974	if ind in [3,4,5]:
975	name = names[ind].replace('A','dA')
976	else:
977	name = names[ind]
978	if name in sigDict:
979	atomsSig[str(i)+':'+str(ind)] = sigDict[name]
980	if at[9] == 'I':
981	at[10] = parmDict[names[10]]
982	if names[10] in sigDict:
983	atomsSig[str(i)+':10'] = sigDict[names[10]]
984	else:
985	for ind in [11,12,13,14,15,16]:
986	at[ind] = parmDict[names[ind]]
987	if names[ind] in sigDict:
988	atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
989	PrintAtomsAndSig(General,Atoms,atomsSig)
990
991	textureData = General['SH Texture']
992	if textureData['Order']:
993	SHtextureSig = {}
994	for name in ['omega','chi','phi']:
995	aname = pfx+'SH '+name
996	textureData['Sample '+name][1] = parmDict[aname]
997	if aname in sigDict:
998	SHtextureSig['Sample '+name] = sigDict[aname]
999	for name in textureData['SH Coeff'][1]:
1000	aname = pfx+name
1001	textureData['SH Coeff'][1][name] = parmDict[aname]
1002	if aname in sigDict:
1003	SHtextureSig[name] = sigDict[aname]
1004	PrintSHtextureAndSig(textureData,SHtextureSig)
1005
1006	################################################################################
1007	##### Histogram & Phase data
1008	################################################################################
1009
1010	def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None):
1011
1012	def PrintSize(hapData):
1013	if hapData[0] in ['isotropic','uniaxial']:
1014	line = '\n Size model : %9s'%(hapData[0])
1015	line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1016	if hapData[0] == 'uniaxial':
1017	line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1018	line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1019	print >>pFile,line
1020	else:
1021	print >>pFile,'\n Size model : %s'%(hapData[0])+ \
1022	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1023	Snames = ['S11','S22','S33','S12','S13','S23']
1024	ptlbls = ' names :'
1025	ptstr = ' values:'
1026	varstr = ' refine:'
1027	for i,name in enumerate(Snames):
1028	ptlbls += '%12s' % (name)
1029	ptstr += '%12.6f' % (hapData[4][i])
1030	varstr += '%12s' % (str(hapData[5][i]))
1031	print >>pFile,ptlbls
1032	print >>pFile,ptstr
1033	print >>pFile,varstr
1034
1035	def PrintMuStrain(hapData,SGData):
1036	if hapData[0] in ['isotropic','uniaxial']:
1037	line = '\n Mustrain model: %9s'%(hapData[0])
1038	line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1039	if hapData[0] == 'uniaxial':
1040	line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1041	line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1042	print >>pFile,line
1043	else:
1044	print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1045	'\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1046	Snames = G2spc.MustrainNames(SGData)
1047	ptlbls = ' names :'
1048	ptstr = ' values:'
1049	varstr = ' refine:'
1050	for i,name in enumerate(Snames):
1051	ptlbls += '%12s' % (name)
1052	ptstr += '%12.6f' % (hapData[4][i])
1053	varstr += '%12s' % (str(hapData[5][i]))
1054	print >>pFile,ptlbls
1055	print >>pFile,ptstr
1056	print >>pFile,varstr
1057
1058	def PrintHStrain(hapData,SGData):
1059	print >>pFile,'\n Hydrostatic/elastic strain: '
1060	Hsnames = G2spc.HStrainNames(SGData)
1061	ptlbls = ' names :'
1062	ptstr = ' values:'
1063	varstr = ' refine:'
1064	for i,name in enumerate(Hsnames):
1065	ptlbls += '%12s' % (name)
1066	ptstr += '%12.6f' % (hapData[0][i])
1067	varstr += '%12s' % (str(hapData[1][i]))
1068	print >>pFile,ptlbls
1069	print >>pFile,ptstr
1070	print >>pFile,varstr
1071
1072	def PrintSHPO(hapData):
1073	print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1074	str(hapData[4])+' Refine? '+str(hapData[2])
1075	ptlbls = ' names :'
1076	ptstr = ' values:'
1077	for item in hapData[5]:
1078	ptlbls += '%12s'%(item)
1079	ptstr += '%12.3f'%(hapData[5][item])
1080	print >>pFile,ptlbls
1081	print >>pFile,ptstr
1082
1083	hapDict = {}
1084	hapVary = []
1085	controlDict = {}
1086	poType = {}
1087	poAxes = {}
1088	spAxes = {}
1089	spType = {}
1090
1091	for phase in Phases:
1092	HistoPhase = Phases[phase]['Histograms']
1093	SGData = Phases[phase]['General']['SGData']
1094	cell = Phases[phase]['General']['Cell'][1:7]
1095	A = G2lat.cell2A(cell)
1096	pId = Phases[phase]['pId']
1097	histoList = HistoPhase.keys()
1098	histoList.sort()
1099	for histogram in histoList:
1100	try:
1101	Histogram = Histograms[histogram]
1102	except KeyError:
1103	#skip if histogram not included e.g. in a sequential refinement
1104	continue
1105	hapData = HistoPhase[histogram]
1106	hId = Histogram['hId']
1107	if 'PWDR' in histogram:
1108	limits = Histogram['Limits'][1]
1109	inst = Histogram['Instrument Parameters']
1110	inst = dict(zip(inst[3],inst[1]))
1111	Zero = inst['Zero']
1112	if 'C' in inst['Type']:
1113	try:
1114	wave = inst['Lam']
1115	except KeyError:
1116	wave = inst['Lam1']
1117	dmin = wave/(2.0*sind(limits[1]/2.0))
1118	pfx = str(pId)+':'+str(hId)+':'
1119	for item in ['Scale','Extinction']:
1120	hapDict[pfx+item] = hapData[item][0]
1121	if hapData[item][1]:
1122	hapVary.append(pfx+item)
1123	names = G2spc.HStrainNames(SGData)
1124	for i,name in enumerate(names):
1125	hapDict[pfx+name] = hapData['HStrain'][0][i]
1126	if hapData['HStrain'][1][i]:
1127	hapVary.append(pfx+name)
1128	controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1129	if hapData['Pref.Ori.'][0] == 'MD':
1130	hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1131	controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1132	if hapData['Pref.Ori.'][2]:
1133	hapVary.append(pfx+'MD')
1134	else: #'SH' spherical harmonics
1135	controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1136	controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1137	for item in hapData['Pref.Ori.'][5]:
1138	hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1139	if hapData['Pref.Ori.'][2]:
1140	hapVary.append(pfx+item)
1141	for item in ['Mustrain','Size']:
1142	controlDict[pfx+item+'Type'] = hapData[item][0]
1143	hapDict[pfx+item+':mx'] = hapData[item][1][2]
1144	if hapData[item][2][2]:
1145	hapVary.append(pfx+item+':mx')
1146	if hapData[item][0] in ['isotropic','uniaxial']:
1147	hapDict[pfx+item+':i'] = hapData[item][1][0]
1148	if hapData[item][2][0]:
1149	hapVary.append(pfx+item+':i')
1150	if hapData[item][0] == 'uniaxial':
1151	controlDict[pfx+item+'Axis'] = hapData[item][3]
1152	hapDict[pfx+item+':a'] = hapData[item][1][1]
1153	if hapData[item][2][1]:
1154	hapVary.append(pfx+item+':a')
1155	else: #generalized for mustrain or ellipsoidal for size
1156	Nterms = len(hapData[item][4])
1157	if item == 'Mustrain':
1158	names = G2spc.MustrainNames(SGData)
1159	pwrs = []
1160	for name in names:
1161	h,k,l = name[1:]
1162	pwrs.append([int(h),int(k),int(l)])
1163	controlDict[pfx+'MuPwrs'] = pwrs
1164	for i in range(Nterms):
1165	sfx = ':'+str(i)
1166	hapDict[pfx+item+sfx] = hapData[item][4][i]
1167	if hapData[item][5][i]:
1168	hapVary.append(pfx+item+sfx)
1169
1170	if Print:
1171	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1172	print >>pFile,135*'-'
1173	print >>pFile,' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1174	print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1175	if hapData['Pref.Ori.'][0] == 'MD':
1176	Ax = hapData['Pref.Ori.'][3]
1177	print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1178	' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1179	else: #'SH' for spherical harmonics
1180	PrintSHPO(hapData['Pref.Ori.'])
1181	PrintSize(hapData['Size'])
1182	PrintMuStrain(hapData['Mustrain'],SGData)
1183	PrintHStrain(hapData['HStrain'],SGData)
1184	HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1185	refList = []
1186	for h,k,l,d in HKLd:
1187	ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1188	mul *= 2 # for powder overlap of Friedel pairs
1189	if ext:
1190	continue
1191	if 'C' in inst['Type']:
1192	pos = 2.0asind(wave/(2.0d))+Zero
1193	if limits[0] < pos < limits[1]:
1194	refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1195	#last item should contain form factor values by atom type
1196	else:
1197	raise ValueError
1198	Histogram['Reflection Lists'][phase] = refList
1199	elif 'HKLF' in histogram:
1200	inst = Histogram['Instrument Parameters']
1201	inst = dict(zip(inst[2],inst[1]))
1202	hId = Histogram['hId']
1203	hfx = ':%d:'%(hId)
1204	for item in inst:
1205	hapDict[hfx+item] = inst[item]
1206	pfx = str(pId)+':'+str(hId)+':'
1207	hapDict[pfx+'Scale'] = hapData['Scale'][0]
1208	if hapData['Scale'][1]:
1209	hapVary.append(pfx+'Scale')
1210	extApprox,extType,extParms = hapData['Extinction']
1211	controlDict[pfx+'EType'] = extType
1212	controlDict[pfx+'EApprox'] = extApprox
1213	if 'Primary' in extType:
1214	Ekey = ['Ep',]
1215	elif 'Secondary Type II' == extType:
1216	Ekey = ['Es',]
1217	elif 'Secondary Type I' == extType:
1218	Ekey = ['Eg',]
1219	else:
1220	Ekey = ['Eg','Es']
1221	for eKey in Ekey:
1222	hapDict[pfx+eKey] = extParms[eKey][0]
1223	if extParms[eKey][1]:
1224	hapVary.append(pfx+eKey)
1225	if Print:
1226	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1227	print >>pFile,135*'-'
1228	print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1229	print >>pFile,' Extinction approx: %10s'%(extApprox),' Type: %15s'%(extType),' tbar: %6.3f'%(extParms['Tbar'])
1230	text = ' Parameters :'
1231	for eKey in Ekey:
1232	text += ' %4s : %10.3g Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1233	print >>pFile,text
1234	Histogram['Reflection Lists'] = phase
1235
1236	return hapVary,hapDict,controlDict
1237
1238	def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1239
1240	def PrintSizeAndSig(hapData,sizeSig):
1241	line = '\n Size model: %9s'%(hapData[0])
1242	refine = False
1243	if hapData[0] in ['isotropic','uniaxial']:
1244	line += ' equatorial:%12.3f'%(hapData[1][0])
1245	if sizeSig[0][0]:
1246	line += ', sig:%8.3f'%(sizeSig[0][0])
1247	refine = True
1248	if hapData[0] == 'uniaxial':
1249	line += ' axial:%12.3f'%(hapData[1][1])
1250	if sizeSig[0][1]:
1251	refine = True
1252	line += ', sig:%8.3f'%(sizeSig[0][1])
1253	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1254	if sizeSig[0][2]:
1255	refine = True
1256	line += ', sig:%8.3f'%(sizeSig[0][2])
1257	if refine:
1258	print >>pFile,line
1259	else:
1260	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1261	if sizeSig[0][2]:
1262	refine = True
1263	line += ', sig:%8.3f'%(sizeSig[0][2])
1264	Snames = ['S11','S22','S33','S12','S13','S23']
1265	ptlbls = ' name :'
1266	ptstr = ' value :'
1267	sigstr = ' sig :'
1268	for i,name in enumerate(Snames):
1269	ptlbls += '%12s' % (name)
1270	ptstr += '%12.6f' % (hapData[4][i])
1271	if sizeSig[1][i]:
1272	refine = True
1273	sigstr += '%12.6f' % (sizeSig[1][i])
1274	else:
1275	sigstr += 12*' '
1276	if refine:
1277	print >>pFile,line
1278	print >>pFile,ptlbls
1279	print >>pFile,ptstr
1280	print >>pFile,sigstr
1281
1282	def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1283	line = '\n Mustrain model: %9s'%(hapData[0])
1284	refine = False
1285	if hapData[0] in ['isotropic','uniaxial']:
1286	line += ' equatorial:%12.1f'%(hapData[1][0])
1287	if mustrainSig[0][0]:
1288	line += ', sig: %8.1f'%(mustrainSig[0][0])
1289	refine = True
1290	if hapData[0] == 'uniaxial':
1291	line += ' axial:%12.1f'%(hapData[1][1])
1292	if mustrainSig[0][1]:
1293	line += ', sig:%8.1f'%(mustrainSig[0][1])
1294	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1295	if mustrainSig[0][2]:
1296	refine = True
1297	line += ', sig:%8.3f'%(mustrainSig[0][2])
1298	if refine:
1299	print >>pFile,line
1300	else:
1301	line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1302	if mustrainSig[0][2]:
1303	refine = True
1304	line += ', sig:%8.3f'%(mustrainSig[0][2])
1305	Snames = G2spc.MustrainNames(SGData)
1306	ptlbls = ' name :'
1307	ptstr = ' value :'
1308	sigstr = ' sig :'
1309	for i,name in enumerate(Snames):
1310	ptlbls += '%12s' % (name)
1311	ptstr += '%12.6f' % (hapData[4][i])
1312	if mustrainSig[1][i]:
1313	refine = True
1314	sigstr += '%12.6f' % (mustrainSig[1][i])
1315	else:
1316	sigstr += 12*' '
1317	if refine:
1318	print >>pFile,line
1319	print >>pFile,ptlbls
1320	print >>pFile,ptstr
1321	print >>pFile,sigstr
1322
1323	def PrintHStrainAndSig(hapData,strainSig,SGData):
1324	Hsnames = G2spc.HStrainNames(SGData)
1325	ptlbls = ' name :'
1326	ptstr = ' value :'
1327	sigstr = ' sig :'
1328	refine = False
1329	for i,name in enumerate(Hsnames):
1330	ptlbls += '%12s' % (name)
1331	ptstr += '%12.6g' % (hapData[0][i])
1332	if name in strainSig:
1333	refine = True
1334	sigstr += '%12.6g' % (strainSig[name])
1335	else:
1336	sigstr += 12*' '
1337	if refine:
1338	print >>pFile,'\n Hydrostatic/elastic strain: '
1339	print >>pFile,ptlbls
1340	print >>pFile,ptstr
1341	print >>pFile,sigstr
1342
1343	def PrintSHPOAndSig(hapData,POsig):
1344	print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1345	ptlbls = ' names :'
1346	ptstr = ' values:'
1347	sigstr = ' sig :'
1348	for item in hapData[5]:
1349	ptlbls += '%12s'%(item)
1350	ptstr += '%12.3f'%(hapData[5][item])
1351	if item in POsig:
1352	sigstr += '%12.3f'%(POsig[item])
1353	else:
1354	sigstr += 12*' '
1355	print >>pFile,ptlbls
1356	print >>pFile,ptstr
1357	print >>pFile,sigstr
1358
1359	for phase in Phases:
1360	HistoPhase = Phases[phase]['Histograms']
1361	SGData = Phases[phase]['General']['SGData']
1362	pId = Phases[phase]['pId']
1363	histoList = HistoPhase.keys()
1364	histoList.sort()
1365	for histogram in histoList:
1366	try:
1367	Histogram = Histograms[histogram]
1368	except KeyError:
1369	#skip if histogram not included e.g. in a sequential refinement
1370	continue
1371	print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1372	print >>pFile,130*'-'
1373	hapData = HistoPhase[histogram]
1374	hId = Histogram['hId']
1375	if 'PWDR' in histogram:
1376	pfx = str(pId)+':'+str(hId)+':'
1377	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
1378	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1379
1380	PhFrExtPOSig = {}
1381	for item in ['Scale','Extinction']:
1382	hapData[item][0] = parmDict[pfx+item]
1383	if pfx+item in sigDict:
1384	PhFrExtPOSig[item] = sigDict[pfx+item]
1385	if hapData['Pref.Ori.'][0] == 'MD':
1386	hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1387	if pfx+'MD' in sigDict:
1388	PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1389	else: #'SH' spherical harmonics
1390	for item in hapData['Pref.Ori.'][5]:
1391	hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1392	if pfx+item in sigDict:
1393	PhFrExtPOSig[item] = sigDict[pfx+item]
1394	if Print:
1395	if 'Scale' in PhFrExtPOSig:
1396	print >>pFile,' Phase fraction : %10.5f, sig %10.5f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1397	if 'Extinction' in PhFrExtPOSig:
1398	print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1399	if hapData['Pref.Ori.'][0] == 'MD':
1400	if 'MD' in PhFrExtPOSig:
1401	print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1402	else:
1403	PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1404	SizeMuStrSig = {'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1405	'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1406	'HStrain':{}}
1407	for item in ['Mustrain','Size']:
1408	hapData[item][1][2] = parmDict[pfx+item+':mx']
1409	hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1410	if pfx+item+':mx' in sigDict:
1411	SizeMuStrSig[item][0][2] = sigDict[pfx+item+':mx']
1412	if hapData[item][0] in ['isotropic','uniaxial']:
1413	hapData[item][1][0] = parmDict[pfx+item+':i']
1414	if item == 'Size':
1415	hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1416	if pfx+item+':i' in sigDict:
1417	SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
1418	if hapData[item][0] == 'uniaxial':
1419	hapData[item][1][1] = parmDict[pfx+item+':a']
1420	if item == 'Size':
1421	hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))
1422	if pfx+item+':a' in sigDict:
1423	SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
1424	else: #generalized for mustrain or ellipsoidal for size
1425	Nterms = len(hapData[item][4])
1426	for i in range(Nterms):
1427	sfx = ':'+str(i)
1428	hapData[item][4][i] = parmDict[pfx+item+sfx]
1429	if pfx+item+sfx in sigDict:
1430	SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1431	names = G2spc.HStrainNames(SGData)
1432	for i,name in enumerate(names):
1433	hapData['HStrain'][0][i] = parmDict[pfx+name]
1434	if pfx+name in sigDict:
1435	SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1436	if Print:
1437	PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1438	PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1439	PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1440
1441	elif 'HKLF' in histogram:
1442	pfx = str(pId)+':'+str(hId)+':'
1443	print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \
1444	%(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1445	print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
1446	ScalExtSig = {}
1447	for item in ['Scale','Ep','Eg','Es']:
1448	if parmDict.get(pfx+item):
1449	hapData[item][0] = parmDict[pfx+item]
1450	if pfx+item in sigDict:
1451	ScalExtSig[item] = sigDict[pfx+item]
1452	if Print:
1453	if 'Scale' in ScalExtSig:
1454	print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig['Scale'])
1455	# fix after it runs!
1456	# print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1457	# print >>pFile,135*'-'
1458	# print >>pFile,' Scale factor : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1459	# print >>pFile,' Extinction approx: %10s'%(extApprox),' Type: %15s'%(extType),' tbar: %6.3f'%(extParms['Tbar'])
1460	# text = ' Parameters :'
1461	# for eKey in Ekey:
1462	# text += ' %4s : %10.3g Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1463	# print >>pFile,text
1464
1465	################################################################################
1466	##### Histogram data
1467	################################################################################
1468
1469	def GetHistogramData(Histograms,Print=True,pFile=None):
1470
1471	def GetBackgroundParms(hId,Background):
1472	Back = Background[0]
1473	DebyePeaks = Background[1]
1474	bakType,bakFlag = Back[:2]
1475	backVals = Back[3:]
1476	backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1477	backDict = dict(zip(backNames,backVals))
1478	backVary = []
1479	if bakFlag:
1480	backVary = backNames
1481	backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
1482	backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
1483	debyeDict = {}
1484	debyeList = []
1485	for i in range(DebyePeaks['nDebye']):
1486	debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1487	debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
1488	debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
1489	debyeVary = []
1490	for item in debyeList:
1491	if item[1]:
1492	debyeVary.append(item[0])
1493	backDict.update(debyeDict)
1494	backVary += debyeVary
1495	peakDict = {}
1496	peakList = []
1497	for i in range(DebyePeaks['nPeaks']):
1498	peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
1499	':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
1500	peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
1501	peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
1502	peakVary = []
1503	for item in peakList:
1504	if item[1]:
1505	peakVary.append(item[0])
1506	backDict.update(peakDict)
1507	backVary += peakVary
1508	return bakType,backDict,backVary
1509
1510	def GetInstParms(hId,Inst):
1511	insVals,insFlags,insNames = Inst[1:4]
1512	dataType = insVals[0]
1513	instDict = {}
1514	insVary = []
1515	pfx = ':'+str(hId)+':'
1516	for i,flag in enumerate(insFlags):
1517	insName = pfx+insNames[i]
1518	instDict[insName] = insVals[i]
1519	if flag:
1520	insVary.append(insName)
1521	instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1522	instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1523	instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1524	return dataType,instDict,insVary
1525
1526	def GetSampleParms(hId,Sample):
1527	sampVary = []
1528	hfx = ':'+str(hId)+':'
1529	sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1530	hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1531	Type = Sample['Type']
1532	if 'Bragg' in Type: #Bragg-Brentano
1533	for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering?
1534	sampDict[hfx+item] = Sample[item][0]
1535	if Sample[item][1]:
1536	sampVary.append(hfx+item)
1537	elif 'Debye' in Type: #Debye-Scherrer
1538	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1539	sampDict[hfx+item] = Sample[item][0]
1540	if Sample[item][1]:
1541	sampVary.append(hfx+item)
1542	return Type,sampDict,sampVary
1543
1544	def PrintBackground(Background):
1545	Back = Background[0]
1546	DebyePeaks = Background[1]
1547	print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
1548	line = ' Coefficients: '
1549	for i,back in enumerate(Back[3:]):
1550	line += '%10.3f'%(back)
1551	if i and not i%10:
1552	line += '\n'+15*' '
1553	print >>pFile,line
1554	if DebyePeaks['nDebye']:
1555	print >>pFile,'\n Debye diffuse scattering coefficients'
1556	parms = ['DebyeA','DebyeR','DebyeU']
1557	line = ' names : '
1558	for parm in parms:
1559	line += '%8s refine?'%(parm)
1560	print >>pFile,line
1561	for j,term in enumerate(DebyePeaks['debyeTerms']):
1562	line = ' term'+'%2d'%(j)+':'
1563	for i in range(3):
1564	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
1565	print >>pFile,line
1566	if DebyePeaks['nPeaks']:
1567	print >>pFile,'\n Single peak coefficients'
1568	parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1569	line = ' names : '
1570	for parm in parms:
1571	line += '%8s refine?'%(parm)
1572	print >>pFile,line
1573	for j,term in enumerate(DebyePeaks['peaksList']):
1574	line = ' peak'+'%2d'%(j)+':'
1575	for i in range(4):
1576	line += '%10.3f %5s'%(term[2i],bool(term[2i+1]))
1577	print >>pFile,line
1578
1579	def PrintInstParms(Inst):
1580	print >>pFile,'\n Instrument Parameters:'
1581	ptlbls = ' name :'
1582	ptstr = ' value :'
1583	varstr = ' refine:'
1584	instNames = Inst[3][1:]
1585	for i,name in enumerate(instNames):
1586	ptlbls += '%12s' % (name)
1587	ptstr += '%12.6f' % (Inst[1][i+1])
1588	if name in ['Lam1','Lam2','Azimuth']:
1589	varstr += 12*' '
1590	else:
1591	varstr += '%12s' % (str(bool(Inst[2][i+1])))
1592	print >>pFile,ptlbls
1593	print >>pFile,ptstr
1594	print >>pFile,varstr
1595
1596	def PrintSampleParms(Sample):
1597	print >>pFile,'\n Sample Parameters:'
1598	print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1599	(Sample['Omega'],Sample['Chi'],Sample['Phi'])
1600	ptlbls = ' name :'
1601	ptstr = ' value :'
1602	varstr = ' refine:'
1603	if 'Bragg' in Sample['Type']:
1604	for item in ['Scale','Shift','Transparency']:
1605	ptlbls += '%14s'%(item)
1606	ptstr += '%14.4f'%(Sample[item][0])
1607	varstr += '%14s'%(str(bool(Sample[item][1])))
1608
1609	elif 'Debye' in Type: #Debye-Scherrer
1610	for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1611	ptlbls += '%14s'%(item)
1612	ptstr += '%14.4f'%(Sample[item][0])
1613	varstr += '%14s'%(str(bool(Sample[item][1])))
1614
1615	print >>pFile,ptlbls
1616	print >>pFile,ptstr
1617	print >>pFile,varstr
1618
1619
1620	histDict = {}
1621	histVary = []
1622	controlDict = {}
1623	histoList = Histograms.keys()
1624	histoList.sort()
1625	for histogram in histoList:
1626	if 'PWDR' in histogram:
1627	Histogram = Histograms[histogram]
1628	hId = Histogram['hId']
1629	pfx = ':'+str(hId)+':'
1630	controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
1631	controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1632
1633	Background = Histogram['Background']
1634	Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1635	controlDict[pfx+'bakType'] = Type
1636	histDict.update(bakDict)
1637	histVary += bakVary
1638
1639	Inst = Histogram['Instrument Parameters']
1640	Type,instDict,insVary = GetInstParms(hId,Inst)
1641	controlDict[pfx+'histType'] = Type
1642	if pfx+'Lam1' in instDict:
1643	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1644	else:
1645	controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']
1646	histDict.update(instDict)
1647	histVary += insVary
1648
1649	Sample = Histogram['Sample Parameters']
1650	Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1651	controlDict[pfx+'instType'] = Type
1652	histDict.update(sampDict)
1653	histVary += sampVary
1654
1655
1656	if Print:
1657	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
1658	print >>pFile,135*'-'
1659	Units = {'C':' deg','T':' msec'}
1660	units = Units[controlDict[pfx+'histType'][2]]
1661	Limits = controlDict[pfx+'Limits']
1662	print >>pFile,' Instrument type: ',Sample['Type']
1663	print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
1664	PrintSampleParms(Sample)
1665	PrintInstParms(Inst)
1666	PrintBackground(Background)
1667	elif 'HKLF' in histogram:
1668	Histogram = Histograms[histogram]
1669	hId = Histogram['hId']
1670	pfx = ':'+str(hId)+':'
1671	controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
1672	Inst = Histogram['Instrument Parameters']
1673	controlDict[pfx+'histType'] = Inst[1][0]
1674	histDict[pfx+'Lam'] = Inst[1][1]
1675	controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
1676	return histVary,histDict,controlDict
1677
1678	def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
1679
1680	def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1681	Back = Background[0]
1682	DebyePeaks = Background[1]
1683	lenBack = len(Back[3:])
1684	backSig = [0 for i in range(lenBack+3DebyePeaks['nDebye']+4DebyePeaks['nPeaks'])]
1685	for i in range(lenBack):
1686	Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1687	if pfx+'Back:'+str(i) in sigDict:
1688	backSig[i] = sigDict[pfx+'Back:'+str(i)]
1689	if DebyePeaks['nDebye']:
1690	for i in range(DebyePeaks['nDebye']):
1691	names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1692	for j,name in enumerate(names):
1693	DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
1694	if name in sigDict:
1695	backSig[lenBack+3*i+j] = sigDict[name]
1696	if DebyePeaks['nPeaks']:
1697	for i in range(DebyePeaks['nPeaks']):
1698	names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
1699	pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
1700	for j,name in enumerate(names):
1701	DebyePeaks['peaksList'][i][2*j] = parmDict[name]
1702	if name in sigDict:
1703	backSig[lenBack+3DebyePeaks['nDebye']+4i+j] = sigDict[name]
1704	return backSig
1705
1706	def SetInstParms(pfx,Inst,parmDict,sigDict):
1707	insVals,insFlags,insNames = Inst[1:4]
1708	instSig = [0 for i in range(len(insVals))]
1709	for i,flag in enumerate(insFlags):
1710	insName = pfx+insNames[i]
1711	insVals[i] = parmDict[insName]
1712	if insName in sigDict:
1713	instSig[i] = sigDict[insName]
1714	return instSig
1715
1716	def SetSampleParms(pfx,Sample,parmDict,sigDict):
1717	if 'Bragg' in Sample['Type']: #Bragg-Brentano
1718	sampSig = [0 for i in range(3)]
1719	for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering?
1720	Sample[item][0] = parmDict[pfx+item]
1721	if pfx+item in sigDict:
1722	sampSig[i] = sigDict[pfx+item]
1723	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1724	sampSig = [0 for i in range(4)]
1725	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1726	Sample[item][0] = parmDict[pfx+item]
1727	if pfx+item in sigDict:
1728	sampSig[i] = sigDict[pfx+item]
1729	return sampSig
1730
1731	def PrintBackgroundSig(Background,backSig):
1732	Back = Background[0]
1733	DebyePeaks = Background[1]
1734	lenBack = len(Back[3:])
1735	valstr = ' value : '
1736	sigstr = ' sig : '
1737	refine = False
1738	for i,back in enumerate(Back[3:]):
1739	valstr += '%10.4g'%(back)
1740	if Back[1]:
1741	refine = True
1742	sigstr += '%10.4g'%(backSig[i])
1743	else:
1744	sigstr += 10*' '
1745	if refine:
1746	print >>pFile,'\n Background function: ',Back[0]
1747	print >>pFile,valstr
1748	print >>pFile,sigstr
1749	if DebyePeaks['nDebye']:
1750	ifAny = False
1751	ptfmt = "%12.3f"
1752	names = ' names :'
1753	ptstr = ' values:'
1754	sigstr = ' esds :'
1755	for item in sigDict:
1756	if 'Debye' in item:
1757	ifAny = True
1758	names += '%12s'%(item)
1759	ptstr += ptfmt%(parmDict[item])
1760	sigstr += ptfmt%(sigDict[item])
1761	if ifAny:
1762	print >>pFile,'\n Debye diffuse scattering coefficients'
1763	print >>pFile,names
1764	print >>pFile,ptstr
1765	print >>pFile,sigstr
1766	if DebyePeaks['nPeaks']:
1767	ifAny = False
1768	ptfmt = "%14.3f"
1769	names = ' names :'
1770	ptstr = ' values:'
1771	sigstr = ' esds :'
1772	for item in sigDict:
1773	if 'BkPk' in item:
1774	ifAny = True
1775	names += '%14s'%(item)
1776	ptstr += ptfmt%(parmDict[item])
1777	sigstr += ptfmt%(sigDict[item])
1778	if ifAny:
1779	print >>pFile,'\n Single peak coefficients'
1780	print >>pFile,names
1781	print >>pFile,ptstr
1782	print >>pFile,sigstr
1783
1784	def PrintInstParmsSig(Inst,instSig):
1785	ptlbls = ' names :'
1786	ptstr = ' value :'
1787	sigstr = ' sig :'
1788	instNames = Inst[3][1:]
1789	refine = False
1790	for i,name in enumerate(instNames):
1791	ptlbls += '%12s' % (name)
1792	ptstr += '%12.6f' % (Inst[1][i+1])
1793	if instSig[i+1]:
1794	refine = True
1795	sigstr += '%12.6f' % (instSig[i+1])
1796	else:
1797	sigstr += 12*' '
1798	if refine:
1799	print >>pFile,'\n Instrument Parameters:'
1800	print >>pFile,ptlbls
1801	print >>pFile,ptstr
1802	print >>pFile,sigstr
1803
1804	def PrintSampleParmsSig(Sample,sampleSig):
1805	ptlbls = ' names :'
1806	ptstr = ' values:'
1807	sigstr = ' sig :'
1808	refine = False
1809	if 'Bragg' in Sample['Type']:
1810	for i,item in enumerate(['Scale','Shift','Transparency']):
1811	ptlbls += '%14s'%(item)
1812	ptstr += '%14.4f'%(Sample[item][0])
1813	if sampleSig[i]:
1814	refine = True
1815	sigstr += '%14.4f'%(sampleSig[i])
1816	else:
1817	sigstr += 14*' '
1818
1819	elif 'Debye' in Sample['Type']: #Debye-Scherrer
1820	for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1821	ptlbls += '%14s'%(item)
1822	ptstr += '%14.4f'%(Sample[item][0])
1823	if sampleSig[i]:
1824	refine = True
1825	sigstr += '%14.4f'%(sampleSig[i])
1826	else:
1827	sigstr += 14*' '
1828
1829	if refine:
1830	print >>pFile,'\n Sample Parameters:'
1831	print >>pFile,ptlbls
1832	print >>pFile,ptstr
1833	print >>pFile,sigstr
1834
1835	histoList = Histograms.keys()
1836	histoList.sort()
1837	for histogram in histoList:
1838	if 'PWDR' in histogram:
1839	Histogram = Histograms[histogram]
1840	hId = Histogram['hId']
1841	pfx = ':'+str(hId)+':'
1842	Background = Histogram['Background']
1843	backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1844
1845	Inst = Histogram['Instrument Parameters']
1846	instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1847
1848	Sample = Histogram['Sample Parameters']
1849	sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1850
1851	print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
1852	print >>pFile,135*'-'
1853	print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
1854	print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
1855	if Print:
1856	print >>pFile,' Instrument type: ',Sample['Type']
1857	PrintSampleParmsSig(Sample,sampSig)
1858	PrintInstParmsSig(Inst,instSig)
1859	PrintBackgroundSig(Background,backSig)
1860
1861	################################################################################
1862	##### Penalty & restraint functions
1863	################################################################################
1864
1865	def penaltyFxn(parmDict,varyList):
1866	pFxn = np.zeros(len(varyList))
1867	for i,item in enumerate(varyList):
1868	if 'PWLref' in item and parmDict[item] < 0.:
1869	pFxn[i] = -parmDict[item]**2 #checked OK
1870	return pFxn/1000.
1871
1872	def penaltyDeriv(parmDict,varyList):
1873	pDerv = np.zeros(len(varyList))
1874	for i,item in enumerate(varyList):
1875	if 'PWLref' in item and parmDict[item] < 0.:
1876	pDerv[i] += 2.*parmDict[item]
1877	return pDerv/1000.
1878
1879	################################################################################
1880	##### Function & derivative calculations
1881	################################################################################
1882
1883	def GetAtomFXU(pfx,calcControls,parmDict):
1884	Natoms = calcControls['Natoms'][pfx]
1885	Tdata = Natoms*[' ',]
1886	Mdata = np.zeros(Natoms)
1887	IAdata = Natoms*[' ',]
1888	Fdata = np.zeros(Natoms)
1889	FFdata = []
1890	BLdata = []
1891	Xdata = np.zeros((3,Natoms))
1892	dXdata = np.zeros((3,Natoms))
1893	Uisodata = np.zeros(Natoms)
1894	Uijdata = np.zeros((6,Natoms))
1895	keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1896	'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1897	'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1898	'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1899	'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1900	for iatm in range(Natoms):
1901	for key in keys:
1902	parm = pfx+key+str(iatm)
1903	if parm in parmDict:
1904	keys[key][iatm] = parmDict[parm]
1905	return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1906
1907	def getFFvalues(FFtables,SQ):
1908	FFvals = {}
1909	for El in FFtables:
1910	FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
1911	return FFvals
1912
1913	def getBLvalues(BLtables):
1914	BLvals = {}
1915	for El in BLtables:
1916	BLvals[El] = BLtables[El][1][1]
1917	return BLvals
1918
1919	def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1920	''' Compute structure factors for all h,k,l for phase
1921	input:
1922	refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1923	G: reciprocal metric tensor
1924	pfx: phase id string
1925	SGData: space group info. dictionary output from SpcGroup
1926	calcControls:
1927	ParmDict:
1928	puts result F^2 in each ref[8] in refList
1929	'''
1930	twopi = 2.0*np.pi
1931	twopisq = 2.0np.pi*2
1932	ast = np.sqrt(np.diag(G))
1933	Mast = twopisq*np.multiply.outer(ast,ast)
1934	FFtables = calcControls['FFtables']
1935	BLtables = calcControls['BLtables']
1936	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1937	FF = np.zeros(len(Tdata))
1938	if 'N' in parmDict[hfx+'Type']:
1939	FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
1940	else:
1941	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1942	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1943	maxPos = len(SGData['SGOps'])
1944	Uij = np.array(G2lat.U6toUij(Uijdata))
1945	bij = Mast*Uij.T
1946	for refl in refList:
1947	fbs = np.array([0,0])
1948	H = refl[:3]
1949	SQ = 1./(2.refl[4])*2
1950	if not len(refl[-1]): #no form factors
1951	if 'N' in parmDict[hfx+'Type']:
1952	refl[-1] = getBLvalues(BLtables)
1953	else: #'X'
1954	refl[-1] = getFFvalues(FFtables,SQ)
1955	for i,El in enumerate(Tdata):
1956	FF[i] = refl[-1][El]
1957	SQfactor = 4.0SQtwopisq
1958	Uniq = refl[11]
1959	phi = refl[12]
1960	phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1961	sinp = np.sin(phase)
1962	cosp = np.cos(phase)
1963	occ = Mdata*Fdata/len(Uniq)
1964	biso = -SQfactor*Uisodata
1965	Tiso = np.where(biso<1.,np.exp(biso),1.0)
1966	HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1967	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1968	Tcorr = Tiso*Tuij
1969	fa = np.array([(FF+FP)occcospTcorr,-FPPoccsinpTcorr])
1970	fas = np.sum(np.sum(fa,axis=1),axis=1) #real
1971	if not SGData['SGInv']:
1972	fb = np.array([(FF+FP)occsinpTcorr,FPPocccospTcorr])
1973	fbs = np.sum(np.sum(fb,axis=1),axis=1)
1974	fasq = fas**2
1975	fbsq = fbs**2 #imaginary
1976	refl[9] = np.sum(fasq)+np.sum(fbsq)
1977	refl[10] = atan2d(fbs[0],fas[0])
1978	return refList
1979
1980	def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1981	twopi = 2.0*np.pi
1982	twopisq = 2.0np.pi*2
1983	ast = np.sqrt(np.diag(G))
1984	Mast = twopisq*np.multiply.outer(ast,ast)
1985	FFtables = calcControls['FFtables']
1986	BLtables = calcControls['BLtables']
1987	Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1988	FF = np.zeros(len(Tdata))
1989	if 'N' in parmDict[hfx+'Type']:
1990	FP = 0.
1991	FPP = 0.
1992	else:
1993	FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1994	FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1995	maxPos = len(SGData['SGOps'])
1996	Uij = np.array(G2lat.U6toUij(Uijdata))
1997	bij = Mast*Uij.T
1998	dFdvDict = {}
1999	dFdfr = np.zeros((len(refList),len(Mdata)))
2000	dFdx = np.zeros((len(refList),len(Mdata),3))
2001	dFdui = np.zeros((len(refList),len(Mdata)))
2002	dFdua = np.zeros((len(refList),len(Mdata),6))
2003	for iref,refl in enumerate(refList):
2004	H = np.array(refl[:3])
2005	SQ = 1./(2.refl[4])2 # or (sin(theta)/lambda)*2
2006	for i,El in enumerate(Tdata):
2007	FF[i] = refl[-1][El]
2008	SQfactor = 8.0SQnp.pi**2
2009	Uniq = refl[11]
2010	phi = refl[12]
2011	phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
2012	sinp = np.sin(phase)
2013	cosp = np.cos(phase)
2014	occ = Mdata*Fdata/len(Uniq)
2015	biso = -SQfactor*Uisodata
2016	Tiso = np.where(biso<1.,np.exp(biso),1.0)
2017	HbH = -np.inner(H,np.inner(bij,H))
2018	Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
2019	Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
2020	Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2021	Tcorr = Tiso*Tuij
2022	fot = (FF+FP)occTcorr
2023	fotp = FPPoccTcorr
2024	fa = np.array([fot[:,np.newaxis]cosp,fotp[:,np.newaxis]cosp]) #non positions
2025	fb = np.array([fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp])
2026
2027	fas = np.sum(np.sum(fa,axis=1),axis=1)
2028	fbs = np.sum(np.sum(fb,axis=1),axis=1)
2029	fax = np.array([-fot[:,np.newaxis]sinp,-fotp[:,np.newaxis]sinp]) #positions
2030	fbx = np.array([fot[:,np.newaxis]cosp,-fot[:,np.newaxis]cosp])
2031	#sum below is over Uniq
2032	dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
2033	dfadx = np.sum(twopiUniqfax[:,:,:,np.newaxis],axis=2)
2034	dfadui = np.sum(-SQfactor*fa,axis=2)
2035	dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
2036	#NB: the above have been checked against PA(1:10,1:2) in strfctr.for
2037	dFdfr[iref] = 2.(fas[0]dfadfr[0]+fas[1]dfadfr[1])Mdata/len(Uniq)
2038	dFdx[iref] = 2.(fas[0]dfadx[0]+fas[1]*dfadx[1])
2039	dFdui[iref] = 2.(fas[0]dfadui[0]+fas[1]*dfadui[1])
2040	dFdua[iref] = 2.(fas[0]dfadua[0]+fas[1]*dfadua[1])
2041	if not SGData['SGInv']:
2042	dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom
2043	dfbdx = np.sum(twopiUniqfbx[:,:,:,np.newaxis],axis=2)
2044	dfbdui = np.sum(-SQfactor*fb,axis=2)
2045	dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
2046	dFdfr[iref] += 2.(fbs[0]dfbdfr[0]-fbs[1]dfbdfr[1])Mdata/len(Uniq)
2047	dFdx[iref] += 2.(fbs[0]dfbdx[0]+fbs[1]*dfbdx[1])
2048	dFdui[iref] += 2.(fbs[0]dfbdui[0]-fbs[1]*dfbdui[1])
2049	dFdua[iref] += 2.(fbs[0]dfbdua[0]+fbs[1]*dfbdua[1])
2050	#loop over atoms - each dict entry is list of derivatives for all the reflections
2051	for i in range(len(Mdata)):
2052	dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
2053	dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
2054	dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
2055	dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
2056	dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
2057	dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
2058	dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
2059	dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
2060	dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
2061	dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
2062	dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
2063	return dFdvDict
2064
2065	def Dict2Values(parmdict, varylist):
2066	'''Use before call to leastsq to setup list of values for the parameters
2067	in parmdict, as selected by key in varylist'''
2068	return [parmdict[key] for key in varylist]
2069
2070	def Values2Dict(parmdict, varylist, values):
2071	''' Use after call to leastsq to update the parameter dictionary with
2072	values corresponding to keys in varylist'''
2073	parmdict.update(zip(varylist,values))
2074
2075	def GetNewCellParms(parmDict,varyList):
2076	newCellDict = {}
2077	Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
2078	for item in varyList:
2079	keys = item.split(':')
2080	if keys[2] in Ddict:
2081	key = keys[0]+'::'+Ddict[keys[2]]
2082	parm = keys[0]+'::'+keys[2]
2083	newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
2084	return newCellDict
2085
2086	def ApplyXYZshifts(parmDict,varyList):
2087	''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
2088	input:
2089	parmDict - parameter dictionary
2090	varyList - list of variables
2091	returns:
2092	newAtomDict - dictitemionary of new atomic coordinate names & values;
2093	key is parameter shift name
2094	'''
2095	newAtomDict = {}
2096	for item in parmDict:
2097	if 'dA' in item:
2098	parm = ''.join(item.split('d'))
2099	parmDict[parm] += parmDict[item]
2100	newAtomDict[item] = [parm,parmDict[parm]]
2101	return newAtomDict
2102
2103	def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2104	#Spherical harmonics texture
2105	IFCoup = 'Bragg' in calcControls[hfx+'instType']
2106	odfCor = 1.0
2107	H = refl[:3]
2108	cell = G2lat.Gmat2cell(g)
2109	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2110	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2111	phi,beta = G2lat.CrsAng(H,cell,SGData)
2112	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
2113	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2114	for item in SHnames:
2115	L,M,N = eval(item.strip('C'))
2116	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2117	Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2118	Lnorm = G2lat.Lnorm(L)
2119	odfCor += parmDict[pfx+item]LnormKcl*Ksl
2120	return odfCor
2121
2122	def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2123	#Spherical harmonics texture derivatives
2124	FORPI = 12.5663706143592
2125	IFCoup = 'Bragg' in calcControls[hfx+'instType']
2126	odfCor = 1.0
2127	dFdODF = {}
2128	dFdSA = [0,0,0]
2129	H = refl[:3]
2130	cell = G2lat.Gmat2cell(g)
2131	Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2132	Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2133	phi,beta = G2lat.CrsAng(H,cell,SGData)
2134	psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
2135	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2136	for item in SHnames:
2137	L,M,N = eval(item.strip('C'))
2138	Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2139	Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2140	Lnorm = G2lat.Lnorm(L)
2141	odfCor += parmDict[pfx+item]LnormKcl*Ksl
2142	dFdODF[pfx+item] = LnormKclKsl
2143	for i in range(3):
2144	dFdSA[i] += parmDict[pfx+item]LnormKcl(dKsdpdPSdA[i]+dKsdg*dGMdA[i])
2145	return odfCor,dFdODF,dFdSA
2146
2147	def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2148	#sphericaql harmonics preferred orientation (cylindrical symmetry only)
2149	odfCor = 1.0
2150	H = refl[:3]
2151	cell = G2lat.Gmat2cell(g)
2152	Sangl = [0.,0.,0.]
2153	if 'Bragg' in calcControls[hfx+'instType']:
2154	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2155	IFCoup = True
2156	else:
2157	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2158	IFCoup = False
2159	phi,beta = G2lat.CrsAng(H,cell,SGData)
2160	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2161	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2162	for item in SHnames:
2163	L,N = eval(item.strip('C'))
2164	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
2165	odfCor += parmDict[phfx+item]LnormKcsl
2166	return odfCor
2167
2168	def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2169	#sphericaql harmonics preferred orientation derivatives (cylindrical symmetry only)
2170	FORPI = 12.5663706143592
2171	odfCor = 1.0
2172	dFdODF = {}
2173	H = refl[:3]
2174	cell = G2lat.Gmat2cell(g)
2175	Sangl = [0.,0.,0.]
2176	if 'Bragg' in calcControls[hfx+'instType']:
2177	Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2178	IFCoup = True
2179	else:
2180	Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2181	IFCoup = False
2182	phi,beta = G2lat.CrsAng(H,cell,SGData)
2183	psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2184	SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2185	for item in SHnames:
2186	L,N = eval(item.strip('C'))
2187	Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
2188	odfCor += parmDict[phfx+item]LnormKcsl
2189	dFdODF[phfx+item] = Kcsl*Lnorm
2190	return odfCor,dFdODF
2191
2192	def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2193	POcorr = 1.0
2194	if calcControls[phfx+'poType'] == 'MD':
2195	MD = parmDict[phfx+'MD']
2196	if MD != 1.0:
2197	MDAxis = calcControls[phfx+'MDAxis']
2198	sumMD = 0
2199	for H in refl[11]:
2200	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2201	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
2202	sumMD += A**3
2203	POcorr = sumMD/len(refl[11])
2204	else: #spherical harmonics
2205	if calcControls[pfx+'SHord']:
2206	POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2207	return POcorr
2208
2209	def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2210	POcorr = 1.0
2211	POderv = {}
2212	if calcControls[phfx+'poType'] == 'MD':
2213	MD = parmDict[phfx+'MD']
2214	MDAxis = calcControls[phfx+'MDAxis']
2215	sumMD = 0
2216	sumdMD = 0
2217	for H in refl[11]:
2218	cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2219	A = 1.0/np.sqrt((MDcosP)2+sinP*2/MD)
2220	sumMD += A**3
2221	sumdMD -= (1.5A5)(2.0MDcosP2-(sinP/MD)2)
2222	POcorr = sumMD/len(refl[11])
2223	POderv[phfx+'MD'] = sumdMD/len(refl[11])
2224	else: #spherical harmonics
2225	if calcControls[pfx+'SHord']:
2226	POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2227	return POcorr,POderv
2228
2229	def GetAbsorb(refl,hfx,calcControls,parmDict):
2230	if 'Debye' in calcControls[hfx+'instType']:
2231	return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
2232	else:
2233	return 1.0
2234
2235	def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
2236	if 'Debye' in calcControls[hfx+'instType']:
2237	return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
2238	else:
2239	return 0.0
2240
2241	def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2242	Icorr = parmDict[phfx+'Scale']parmDict[hfx+'Scale']refl[3] #scale*multiplicity
2243	if 'X' in parmDict[hfx+'Type']:
2244	Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
2245	Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2246	if pfx+'SHorder' in parmDict:
2247	Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2248	Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
2249	refl[13] = Icorr
2250
2251	def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2252	dIdsh = 1./parmDict[hfx+'Scale']
2253	dIdsp = 1./parmDict[phfx+'Scale']
2254	if 'X' in parmDict[hfx+'Type']:
2255	pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2256	dIdPola /= pola
2257	else: #'N'
2258	dIdPola = 0.0
2259	POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2260	for iPO in dIdPO:
2261	dIdPO[iPO] /= POcorr
2262	dFdODF = {}
2263	dFdSA = [0,0,0]
2264	if pfx+'SHorder' in parmDict:
2265	odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2266	for iSH in dFdODF:
2267	dFdODF[iSH] /= odfCor
2268	for i in range(3):
2269	dFdSA[i] /= odfCor
2270	dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
2271	return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
2272
2273	def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2274	costh = cosd(refl[5]/2.)
2275	#crystallite size
2276	if calcControls[phfx+'SizeType'] == 'isotropic':
2277	Sgam = 1.8wave/(np.piparmDict[phfx+'Size:i']*costh)
2278	elif calcControls[phfx+'SizeType'] == 'uniaxial':
2279	H = np.array(refl[:3])
2280	P = np.array(calcControls[phfx+'SizeAxis'])
2281	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2282	Sgam = (1.8wave/np.pi)/(parmDict[phfx+'Size:i']parmDict[phfx+'Size:a']*costh)
2283	Sgam = np.sqrt((sinPparmDict[phfx+'Size:a'])*2+(cosPparmDict[phfx+'Size:i'])**2)
2284	else: #ellipsoidal crystallites
2285	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2286	H = np.array(refl[:3])
2287	lenR = G2pwd.ellipseSize(H,Sij,GB)
2288	Sgam = 1.8wave/(np.picosth*lenR)
2289	#microstrain
2290	if calcControls[phfx+'MustrainType'] == 'isotropic':
2291	Mgam = 0.018parmDict[phfx+'Mustrain:i']tand(refl[5]/2.)/np.pi
2292	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2293	H = np.array(refl[:3])
2294	P = np.array(calcControls[phfx+'MustrainAxis'])
2295	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2296	Si = parmDict[phfx+'Mustrain:i']
2297	Sa = parmDict[phfx+'Mustrain:a']
2298	Mgam = 0.018SiSatand(refl[5]/2.)/(np.pinp.sqrt((SicosP)2+(SasinP)**2))
2299	else: #generalized - P.W. Stephens model
2300	pwrs = calcControls[phfx+'MuPwrs']
2301	sum = 0
2302	for i,pwr in enumerate(pwrs):
2303	sum += parmDict[phfx+'Mustrain:'+str(i)]refl[0]pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
2304	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
2305	gam = SgamparmDict[phfx+'Size:mx']+MgamparmDict[phfx+'Mustrain:mx']
2306	sig = (Sgam(1.-parmDict[phfx+'Size:mx']))2+(Mgam(1.-parmDict[phfx+'Mustrain:mx']))**2
2307	sig /= ateln2
2308	return sig,gam
2309
2310	def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2311	gamDict = {}
2312	sigDict = {}
2313	costh = cosd(refl[5]/2.)
2314	tanth = tand(refl[5]/2.)
2315	#crystallite size derivatives
2316	if calcControls[phfx+'SizeType'] == 'isotropic':
2317	Sgam = 1.8wave/(np.piparmDict[phfx+'Size:i']*costh)
2318	gamDict[phfx+'Size:i'] = -900.waveparmDict[phfx+'Size:mx']/(np.pi*costh)
2319	sigDict[phfx+'Size:i'] = -1800.Sgamwave(1.-parmDict[phfx+'Size:mx'])2/(np.picosth*ateln2)
2320	elif calcControls[phfx+'SizeType'] == 'uniaxial':
2321	H = np.array(refl[:3])
2322	P = np.array(calcControls[phfx+'SizeAxis'])
2323	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2324	Si = parmDict[phfx+'Size:i']
2325	Sa = parmDict[phfx+'Size:a']
2326	gami = (1.8wave/np.pi)/(SiSa)
2327	sqtrm = np.sqrt((sinPSa)2+(cosPSi)**2)
2328	Sgam = gami*sqtrm
2329	gam = Sgam/costh
2330	dsi = (gamiSicosP*2/(sqtrmcosth)-gam/Si)
2331	dsa = (gamiSasinP*2/(sqtrmcosth)-gam/Sa)
2332	gamDict[phfx+'Size:i'] = dsi*parmDict[phfx+'Size:mx']
2333	gamDict[phfx+'Size:a'] = dsa*parmDict[phfx+'Size:mx']
2334	sigDict[phfx+'Size:i'] = 2.dsiSgam(1.-parmDict[phfx+'Size:mx'])*2/ateln2
2335	sigDict[phfx+'Size:a'] = 2.dsaSgam(1.-parmDict[phfx+'Size:mx'])*2/ateln2
2336	else: #ellipsoidal crystallites
2337	const = 1.8wave/(np.picosth)
2338	Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2339	H = np.array(refl[:3])
2340	lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2341	Sgam = 1.8wave/(np.picosth*lenR)
2342	for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2343	gamDict[item] = -(const/lenR*2)dRdS[i]*parmDict[phfx+'Size:mx']
2344	sigDict[item] = -2.Sgam(const/lenR*2)dRdS[i](1.-parmDict[phfx+'Size:mx'])*2/ateln2
2345	gamDict[phfx+'Size:mx'] = Sgam
2346	sigDict[phfx+'Size:mx'] = -2.Sgam2(1.-parmDict[phfx+'Size:mx'])/ateln2
2347
2348	#microstrain derivatives
2349	if calcControls[phfx+'MustrainType'] == 'isotropic':
2350	Mgam = 0.018parmDict[phfx+'Mustrain:i']tand(refl[5]/2.)/np.pi
2351	gamDict[phfx+'Mustrain:i'] = 0.018tanthparmDict[phfx+'Mustrain:mx']/np.pi
2352	sigDict[phfx+'Mustrain:i'] = 0.036Mgamtanth(1.-parmDict[phfx+'Mustrain:mx'])2/(np.piateln2)
2353	elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2354	H = np.array(refl[:3])
2355	P = np.array(calcControls[phfx+'MustrainAxis'])
2356	cosP,sinP = G2lat.CosSinAngle(H,P,G)
2357	Si = parmDict[phfx+'Mustrain:i']
2358	Sa = parmDict[phfx+'Mustrain:a']
2359	gami = 0.018SiSa*tanth/np.pi
2360	sqtrm = np.sqrt((SicosP)2+(SasinP)**2)
2361	Mgam = gami/sqtrm
2362	dsi = -gamiSicosP2/sqtrm3
2363	dsa = -gamiSasinP2/sqtrm3
2364	gamDict[phfx+'Mustrain:i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain:mx']
2365	gamDict[phfx+'Mustrain:a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain:mx']
2366	sigDict[phfx+'Mustrain:i'] = 2(Mgam/Si+dsi)Mgam(1.-parmDict[phfx+'Mustrain:mx'])*2/ateln2
2367	sigDict[phfx+'Mustrain:a'] = 2(Mgam/Sa+dsa)Mgam(1.-parmDict[phfx+'Mustrain:mx'])*2/ateln2
2368	else: #generalized - P.W. Stephens model
2369	pwrs = calcControls[phfx+'MuPwrs']
2370	const = 0.018refl[4]2tanth
2371	sum = 0
2372	for i,pwr in enumerate(pwrs):
2373	term = refl[0]*pwr[0]refl[1]*pwr[1]refl[2]**pwr[2]
2374	sum += parmDict[phfx+'Mustrain:'+str(i)]*term
2375	gamDict[phfx+'Mustrain:'+str(i)] = consttermparmDict[phfx+'Mustrain:mx']
2376	sigDict[phfx+'Mustrain:'+str(i)] = \
2377	2.constterm(1.-parmDict[phfx+'Mustrain:mx'])*2/ateln2
2378	Mgam = 0.018refl[4]2tand(refl[5]/2.)*sum
2379	for i in range(len(pwrs)):
2380	sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
2381	gamDict[phfx+'Mustrain:mx'] = Mgam
2382	sigDict[phfx+'Mustrain:mx'] = -2.Mgam2(1.-parmDict[phfx+'Mustrain:mx'])/ateln2
2383	return sigDict,gamDict
2384
2385	def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2386	h,k,l = refl[:3]
2387	dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2388	d = np.sqrt(dsq)
2389
2390	refl[4] = d
2391	pos = 2.0asind(wave/(2.0d))+parmDict[hfx+'Zero']
2392	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2393	if 'Bragg' in calcControls[hfx+'instType']:
2394	pos -= const(4.parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2395	parmDict[hfx+'Transparency']sind(pos)100.0) #trans(=1/mueff) in cm
2396	else: #Debye-Scherrer - simple but maybe not right
2397	pos -= const(parmDict[hfx+'DisplaceX']cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2398	return pos
2399
2400	def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2401	dpr = 180./np.pi
2402	h,k,l = refl[:3]
2403	dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2404	dst = np.sqrt(dstsq)
2405	pos = refl[5]-parmDict[hfx+'Zero']
2406	const = dpr/np.sqrt(1.0-wave*2dstsq/4.0)
2407	dpdw = const*dst
2408	dpdA = np.array([h2,k2,l*2,hk,hl,kl])
2409	dpdA = constwave/(2.0*dst)
2410	dpdZ = 1.0
2411	const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns
2412	if 'Bragg' in calcControls[hfx+'instType']:
2413	dpdSh = -4.constcosd(pos/2.0)
2414	dpdTr = -constsind(pos)100.0
2415	return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2416	else: #Debye-Scherrer - simple but maybe not right
2417	dpdXd = -const*cosd(pos)
2418	dpdYd = -const*sind(pos)
2419	return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2420
2421	def GetHStrainShift(refl,SGData,phfx,parmDict):
2422	laue = SGData['SGLaue']
2423	uniq = SGData['SGUniq']
2424	h,k,l = refl[:3]
2425	if laue in ['m3','m3m']:
2426	Dij = parmDict[phfx+'D11'](h2+k2+l*2)+ \
2427	refl[4]*2parmDict[phfx+'eA']((hk)*2+(hl)*2+(kl)2)/(h2+k2+l2)**2
2428	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2429	Dij = parmDict[phfx+'D11'](h2+k2+hk)+parmDict[phfx+'D33']l*2
2430	elif laue in ['3R','3mR']:
2431	Dij = parmDict[phfx+'D11'](h2+k2+l2)+parmDict[phfx+'D12'](hk+hl+k*l)
2432	elif laue in ['4/m','4/mmm']:
2433	Dij = parmDict[phfx+'D11'](h2+k2)+parmDict[phfx+'D33']l**2
2434	elif laue in ['mmm']:
2435	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2436	elif laue in ['2/m']:
2437	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2
2438	if uniq == 'a':
2439	Dij += parmDict[phfx+'D23']kl
2440	elif uniq == 'b':
2441	Dij += parmDict[phfx+'D13']hl
2442	elif uniq == 'c':
2443	Dij += parmDict[phfx+'D12']hk
2444	else:
2445	Dij = parmDict[phfx+'D11']h2+parmDict[phfx+'D22']k*2+parmDict[phfx+'D33']l**2+ \
2446	parmDict[phfx+'D12']hk+parmDict[phfx+'D13']hl+parmDict[phfx+'D23']kl
2447	return -Dijrefl[4]2tand(refl[5]/2.0)
2448
2449	def GetHStrainShiftDerv(refl,SGData,phfx):
2450	laue = SGData['SGLaue']
2451	uniq = SGData['SGUniq']
2452	h,k,l = refl[:3]
2453	if laue in ['m3','m3m']:
2454	dDijDict = {phfx+'D11':h2+k2+l**2,
2455	phfx+'eA':refl[4]*2((hk)2+(hl)*2+(kl)2)/(h2+k2+l2)**2}
2456	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2457	dDijDict = {phfx+'D11':h2+k2+hk,phfx+'D33':l*2}
2458	elif laue in ['3R','3mR']:
2459	dDijDict = {phfx+'D11':h2+k2+l*2,phfx+'D12':hk+hl+kl}
2460	elif laue in ['4/m','4/mmm']:
2461	dDijDict = {phfx+'D11':h2+k2,phfx+'D33':l**2}
2462	elif laue in ['mmm']:
2463	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2464	elif laue in ['2/m']:
2465	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2}
2466	if uniq == 'a':
2467	dDijDict[phfx+'D23'] = k*l
2468	elif uniq == 'b':
2469	dDijDict[phfx+'D13'] = h*l
2470	elif uniq == 'c':
2471	dDijDict[phfx+'D12'] = h*k
2472	names.append()
2473	else:
2474	dDijDict = {phfx+'D11':h2,phfx+'D22':k2,phfx+'D33':l**2,
2475	phfx+'D12':hk,phfx+'D13':hl,phfx+'D23':k*l}
2476	for item in dDijDict:
2477	dDijDict[item] = -refl[4]2tand(refl[5]/2.0)
2478	return dDijDict
2479
2480	def GetFprime(controlDict,Histograms):
2481	FFtables = controlDict['FFtables']
2482	if not FFtables:
2483	return
2484	histoList = Histograms.keys()
2485	histoList.sort()
2486	for histogram in histoList:
2487	if histogram[:4] in ['PWDR','HKLF']:
2488	Histogram = Histograms[histogram]
2489	hId = Histogram['hId']
2490	hfx = ':%d:'%(hId)
2491	keV = controlDict[hfx+'keV']
2492	for El in FFtables:
2493	Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2494	FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2495	FFtables[El][hfx+'FP'] = FP
2496	FFtables[El][hfx+'FPP'] = FPP
2497
2498	def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2499	histoList = Histograms.keys()
2500	histoList.sort()
2501	for histogram in histoList:
2502	if 'PWDR' in histogram[:4]:
2503	Histogram = Histograms[histogram]
2504	hId = Histogram['hId']
2505	hfx = ':%d:'%(hId)
2506	Limits = calcControls[hfx+'Limits']
2507	shl = max(parmDict[hfx+'SH/L'],0.002)
2508	Ka2 = False
2509	kRatio = 0.0
2510	if hfx+'Lam1' in parmDict.keys():
2511	Ka2 = True
2512	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2513	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2514	x,y,w,yc,yb,yd = Histogram['Data']
2515	xB = np.searchsorted(x,Limits[0])
2516	xF = np.searchsorted(x,Limits[1])
2517	ymb = np.array(y-yb)
2518	ymb = np.where(ymb,ymb,1.0)
2519	ycmb = np.array(yc-yb)
2520	ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)
2521	refLists = Histogram['Reflection Lists']
2522	for phase in refLists:
2523	Phase = Phases[phase]
2524	pId = Phase['pId']
2525	phfx = '%d:%d:'%(pId,hId)
2526	refList = refLists[phase]
2527	sumFo = 0.0
2528	sumdF = 0.0
2529	sumFosq = 0.0
2530	sumdFsq = 0.0
2531	for refl in refList:
2532	if 'C' in calcControls[hfx+'histType']:
2533	yp = np.zeros_like(yb)
2534	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2535	iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
2536	iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
2537	iFin2 = iFin
2538	yp[iBeg:iFin] = refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
2539	if Ka2:
2540	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
2541	Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2542	iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
2543	iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
2544	if not iBeg2+iFin2: #peak below low limit - skip peak
2545	continue
2546	elif not iBeg2-iFin2: #peak above high limit - done
2547	break
2548	yp[iBeg2:iFin2] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here
2549	refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]ratio[iBeg:iFin2]/(refl[13](1.+kRatio)),0.0))
2550	elif 'T' in calcControls[hfx+'histType']:
2551	print 'TOF Undefined at present'
2552	raise Exception #no TOF yet
2553	Fo = np.sqrt(np.abs(refl[8]))
2554	Fc = np.sqrt(np.abs(refl[9]))
2555	sumFo += Fo
2556	sumFosq += refl[8]**2
2557	sumdF += np.abs(Fo-Fc)
2558	sumdFsq += (refl[8]-refl[9])**2
2559	Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2560	Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2561	Histogram[phfx+'Nref'] = len(refList)
2562
2563	def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2564
2565	def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
2566	U = parmDict[hfx+'U']
2567	V = parmDict[hfx+'V']
2568	W = parmDict[hfx+'W']
2569	X = parmDict[hfx+'X']
2570	Y = parmDict[hfx+'Y']
2571	tanPos = tand(refl[5]/2.0)
2572	Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
2573	sig = UtanPos2+VtanPos+W+Ssig #save peak sigma
2574	sig = max(0.001,sig)
2575	gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam #save peak gamma
2576	gam = max(0.001,gam)
2577	return sig,gam
2578
2579	hId = Histogram['hId']
2580	hfx = ':%d:'%(hId)
2581	bakType = calcControls[hfx+'bakType']
2582	yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2583	yc = np.zeros_like(yb)
2584
2585	if 'C' in calcControls[hfx+'histType']:
2586	shl = max(parmDict[hfx+'SH/L'],0.002)
2587	Ka2 = False
2588	if hfx+'Lam1' in parmDict.keys():
2589	wave = parmDict[hfx+'Lam1']
2590	Ka2 = True
2591	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2592	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2593	else:
2594	wave = parmDict[hfx+'Lam']
2595	else:
2596	print 'TOF Undefined at present'
2597	raise ValueError
2598	for phase in Histogram['Reflection Lists']:
2599	refList = Histogram['Reflection Lists'][phase]
2600	Phase = Phases[phase]
2601	pId = Phase['pId']
2602	pfx = '%d::'%(pId)
2603	phfx = '%d:%d:'%(pId,hId)
2604	hfx = ':%d:'%(hId)
2605	SGData = Phase['General']['SGData']
2606	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2607	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2608	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2609	Vst = np.sqrt(nl.det(G)) #V*
2610	if not Phase['General'].get('doPawley'):
2611	refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2612	for refl in refList:
2613	if 'C' in calcControls[hfx+'histType']:
2614	h,k,l = refl[:3]
2615	refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position
2616	Lorenz = 1./(2.sind(refl[5]/2.)2cosd(refl[5]/2.)) #Lorentz correction
2617	refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift
2618	refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam
2619	GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]
2620	refl[13] = VstLorenz
2621	if Phase['General'].get('doPawley'):
2622	try:
2623	pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2624	parmDict[pInd] = max(parmDict[pInd]/2.,parmDict[pInd])
2625	refl[9] = parmDict[pInd]
2626	except KeyError:
2627	# print ' *Error %d,%d,%d missing from Pawley reflection list *'%(h,k,l)
2628	continue
2629	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2630	iBeg = np.searchsorted(x,refl[5]-fmin)
2631	iFin = np.searchsorted(x,refl[5]+fmax)
2632	if not iBeg+iFin: #peak below low limit - skip peak
2633	continue
2634	elif not iBeg-iFin: #peak above high limit - done
2635	break
2636	yc[iBeg:iFin] += refl[13]refl[9]G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here
2637	if Ka2:
2638	pos2 = refl[5]+lamRatiotand(refl[5]/2.0) # + 360/pi Dlam/lam * tan(th)
2639	Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2640	iBeg = np.searchsorted(x,pos2-fmin)
2641	iFin = np.searchsorted(x,pos2+fmax)
2642	if not iBeg+iFin: #peak below low limit - skip peak
2643	continue
2644	elif not iBeg-iFin: #peak above high limit - done
2645	return yc,yb
2646	yc[iBeg:iFin] += refl[13]refl[9]kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here
2647	elif 'T' in calcControls[hfx+'histType']:
2648	print 'TOF Undefined at present'
2649	raise Exception #no TOF yet
2650	return yc,yb
2651
2652	def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2653
2654	def cellVaryDerv(pfx,SGData,dpdA):
2655	if SGData['SGLaue'] in ['-1',]:
2656	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2657	[pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2658	elif SGData['SGLaue'] in ['2/m',]:
2659	if SGData['SGUniq'] == 'a':
2660	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2661	elif SGData['SGUniq'] == 'b':
2662	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2663	else:
2664	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2665	elif SGData['SGLaue'] in ['mmm',]:
2666	return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2667	elif SGData['SGLaue'] in ['4/m','4/mmm']:
2668	# return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2669	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2670	elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2671	# return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2672	return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2673	elif SGData['SGLaue'] in ['3R', '3mR']:
2674	return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]
2675	elif SGData['SGLaue'] in ['m3m','m3']:
2676	# return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2677	return [[pfx+'A0',dpdA[0]]]
2678	# create a list of dependent variables and set up a dictionary to hold their derivatives
2679	dependentVars = G2mv.GetDependentVars()
2680	depDerivDict = {}
2681	for j in dependentVars:
2682	depDerivDict[j] = np.zeros(shape=(len(x)))
2683	#print 'dependent vars',dependentVars
2684	lenX = len(x)
2685	hId = Histogram['hId']
2686	hfx = ':%d:'%(hId)
2687	bakType = calcControls[hfx+'bakType']
2688	dMdv = np.zeros(shape=(len(varylist),len(x)))
2689	dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2690	if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2691	bBpos =varylist.index(hfx+'Back:0')
2692	dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2693	names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2694	for name in varylist:
2695	if 'Debye' in name:
2696	id = int(name.split(':')[-1])
2697	parm = name[:int(name.rindex(':'))]
2698	ip = names.index(parm)
2699	dMdv[varylist.index(name)] = dMddb[3*id+ip]
2700	names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
2701	for name in varylist:
2702	if 'BkPk' in name:
2703	id = int(name.split(':')[-1])
2704	parm = name[:int(name.rindex(':'))]
2705	ip = names.index(parm)
2706	dMdv[varylist.index(name)] = dMdpk[4*id+ip]
2707	if 'C' in calcControls[hfx+'histType']:
2708	dx = x[1]-x[0]
2709	shl = max(parmDict[hfx+'SH/L'],0.002)
2710	Ka2 = False
2711	if hfx+'Lam1' in parmDict.keys():
2712	wave = parmDict[hfx+'Lam1']
2713	Ka2 = True
2714	lamRatio = 360(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.piparmDict[hfx+'Lam1'])
2715	kRatio = parmDict[hfx+'I(L2)/I(L1)']
2716	else:
2717	wave = parmDict[hfx+'Lam']
2718	else:
2719	print 'TOF Undefined at present'
2720	raise ValueError
2721	for phase in Histogram['Reflection Lists']:
2722	refList = Histogram['Reflection Lists'][phase]
2723	Phase = Phases[phase]
2724	SGData = Phase['General']['SGData']
2725	pId = Phase['pId']
2726	pfx = '%d::'%(pId)
2727	phfx = '%d:%d:'%(pId,hId)
2728	A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2729	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2730	GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies
2731	if not Phase['General'].get('doPawley'):
2732	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2733	for iref,refl in enumerate(refList):
2734	if 'C' in calcControls[hfx+'histType']: #CW powder
2735	h,k,l = refl[:3]
2736	dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2737	Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2738	iBeg = np.searchsorted(x,refl[5]-fmin)
2739	iFin = np.searchsorted(x,refl[5]+fmax)
2740	if not iBeg+iFin: #peak below low limit - skip peak
2741	continue
2742	elif not iBeg-iFin: #peak above high limit - done
2743	break
2744	pos = refl[5]
2745	tanth = tand(pos/2.0)
2746	costh = cosd(pos/2.0)
2747	lenBF = iFin-iBeg
2748	dMdpk = np.zeros(shape=(6,lenBF))
2749	dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2750	for i in range(1,5):
2751	dMdpk[i] += 100.dxrefl[13]refl[9]dMdipk[i]
2752	dMdpk[0] += 100.dxrefl[13]refl[9]dMdipk[0]
2753	dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2754	if Ka2:
2755	pos2 = refl[5]+lamRatiotanth # + 360/pi Dlam/lam * tan(th)
2756	kdelt = int((pos2-refl[5])/dx)
2757	iBeg2 = min(lenX,iBeg+kdelt)
2758	iFin2 = min(lenX,iFin+kdelt)
2759	if iBeg2-iFin2:
2760	lenBF2 = iFin2-iBeg2
2761	dMdpk2 = np.zeros(shape=(6,lenBF2))
2762	dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2763	for i in range(1,5):
2764	dMdpk2[i] = 100.dxrefl[13]refl[9]kRatio*dMdipk2[i]
2765	dMdpk2[0] = 100.dxrefl[13]refl[9]kRatio*dMdipk2[0]
2766	dMdpk2[5] = 100.dxrefl[13]*dMdipk2[0]
2767	dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2768	if Phase['General'].get('doPawley'):
2769	dMdpw = np.zeros(len(x))
2770	try:
2771	pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
2772	idx = varylist.index(pIdx)
2773	# parmDict[pIdx] = max(parmDict[pIdx]/2.,parmDict[pIdx])
2774	# refl[9] = abs(parmDict[pIdx])
2775	dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
2776	if parmDict[pIdx] < 0.:
2777	dMdpw[iBeg:iFin] = 2.*dervDict['int']/refl[9]
2778	if Ka2:
2779	dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
2780	if parmDict[pIdx] < 0.:
2781	dMdpw[iBeg2:iFin2] += 2.*dervDict['int']/refl[9]
2782	dMdv[idx] = dMdpw
2783	except: # ValueError:
2784	pass
2785	dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2786	names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2787	hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2788	hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2789	hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2790	hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2791	hfx+'DisplaceY':[dpdY,'pos'],hfx+'Absorption':[dFdAb,'int'],}
2792	for name in names:
2793	item = names[name]
2794	if name in varylist:
2795	dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2796	if Ka2:
2797	dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2798	elif name in dependentVars:
2799	if Ka2:
2800	depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2801	depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2802	for iPO in dIdPO:
2803	if iPO in varylist:
2804	dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2805	if Ka2:
2806	dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2807	elif iPO in dependentVars:
2808	depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2809	if Ka2:
2810	depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2811	for i,name in enumerate(['omega','chi','phi']):
2812	aname = pfx+'SH '+name
2813	if aname in varylist:
2814	dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2815	if Ka2:
2816	dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2817	elif aname in dependentVars:
2818	depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2819	if Ka2:
2820	depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2821	for iSH in dFdODF:
2822	if iSH in varylist:
2823	dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2824	if Ka2:
2825	dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2826	elif iSH in dependentVars:
2827	depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2828	if Ka2:
2829	depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2830	cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2831	for name,dpdA in cellDervNames:
2832	if name in varylist:
2833	dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2834	if Ka2:
2835	dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2836	elif name in dependentVars:
2837	depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2838	if Ka2:
2839	depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2840	dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2841	for name in dDijDict:
2842	if name in varylist:
2843	dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2844	if Ka2:
2845	dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2846	elif name in dependentVars:
2847	depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2848	if Ka2:
2849	depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2850	sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2851	for name in gamDict:
2852	if name in varylist:
2853	dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2854	if Ka2:
2855	dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2856	elif name in dependentVars:
2857	depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2858	if Ka2:
2859	depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2860	for name in sigDict:
2861	if name in varylist:
2862	dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2863	if Ka2:
2864	dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2865	elif name in dependentVars:
2866	depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2867	if Ka2:
2868	depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2869
2870	elif 'T' in calcControls[hfx+'histType']:
2871	print 'TOF Undefined at present'
2872	raise Exception #no TOF yet
2873	#do atom derivatives - for F,X & U so far
2874	corr = dervDict['int']/refl[9]
2875	if Ka2:
2876	corr2 = dervDict2['int']/refl[9]
2877	for name in varylist+dependentVars:
2878	try:
2879	aname = name.split(pfx)[1][:2]
2880	if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2881	except IndexError:
2882	continue
2883	if name in varylist:
2884	dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2885	if Ka2:
2886	dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2887	elif name in dependentVars:
2888	depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2889	if Ka2:
2890	depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2891	# now process derivatives in constraints
2892	G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2893	return dMdv
2894
2895	def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2896	parmdict.update(zip(varylist,values))
2897	G2mv.Dict2Map(parmdict,varylist)
2898	Histograms,Phases,restraintDict = HistoPhases
2899	nvar = len(varylist)
2900	dMdv = np.empty(0)
2901	histoList = Histograms.keys()
2902	histoList.sort()
2903	for histogram in histoList:
2904	if 'PWDR' in histogram[:4]:
2905	Histogram = Histograms[histogram]
2906	hId = Histogram['hId']
2907	hfx = ':%d:'%(hId)
2908	wtFactor = calcControls[hfx+'wtFactor']
2909	Limits = calcControls[hfx+'Limits']
2910	x,y,w,yc,yb,yd = Histogram['Data']
2911	W = wtFactor*w
2912	xB = np.searchsorted(x,Limits[0])
2913	xF = np.searchsorted(x,Limits[1])
2914	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2915	varylist,Histogram,Phases,calcControls,pawleyLookup)
2916	elif 'HKLF' in histogram[:4]:
2917	Histogram = Histograms[histogram]
2918	nobs = Histogram['Nobs']
2919	phase = Histogram['Reflection Lists']
2920	Phase = Phases[phase]
2921	hId = Histogram['hId']
2922	hfx = ':%d:'%(hId)
2923	wtFactor = calcControls[hfx+'wtFactor']
2924	pfx = '%d::'%(Phase['pId'])
2925	phfx = '%d:%d:'%(Phase['pId'],hId)
2926	SGData = Phase['General']['SGData']
2927	A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
2928	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
2929	refList = Histogram['Data']
2930	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
2931	dMdvh = np.zeros((len(varylist),len(refList)))
2932	for iref,ref in enumerate(refList):
2933	if ref[6] > 0:
2934	if calcControls['F**2']:
2935	if ref[5]/ref[6] >= calcControls['minF/sig']:
2936	for j,var in enumerate(varylist):
2937	if var in dFdvDict:
2938	dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
2939	if phfx+'Scale' in varylist:
2940	dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]
2941	else:
2942	Fo = np.sqrt(ref[5])
2943	Fc = np.sqrt(ref[7])
2944	sig = ref[6]/(2.0*Fo)
2945	if Fo/sig >= calcControls['minF/sig']:
2946	for j,var in enumerate(varylist):
2947	if var in dFdvDict:
2948	dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
2949	if phfx+'Scale' in varylist:
2950	dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]
2951	else:
2952	continue #skip non-histogram entries
2953	if len(dMdv):
2954	dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2955	else:
2956	dMdv = dMdvh
2957
2958	# dpdv = penaltyDeriv(parmdict,varylist)
2959	# if np.any(dpdv):
2960	# dMdv = np.concatenate((dMdv.T,np.outer(dpdv,dpdv))).T
2961	return dMdv
2962
2963	def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2964	parmdict.update(zip(varylist,values))
2965	G2mv.Dict2Map(parmdict,varylist)
2966	Histograms,Phases,restraintDict = HistoPhases
2967	nvar = len(varylist)
2968	Hess = np.empty(0)
2969	histoList = Histograms.keys()
2970	histoList.sort()
2971	for histogram in histoList:
2972	if 'PWDR' in histogram[:4]:
2973	Histogram = Histograms[histogram]
2974	hId = Histogram['hId']
2975	hfx = ':%d:'%(hId)
2976	wtFactor = calcControls[hfx+'wtFactor']
2977	Limits = calcControls[hfx+'Limits']
2978	x,y,w,yc,yb,yd = Histogram['Data']
2979	W = wtFactor*w
2980	dy = y-yc
2981	xB = np.searchsorted(x,Limits[0])
2982	xF = np.searchsorted(x,Limits[1])
2983	dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2984	varylist,Histogram,Phases,calcControls,pawleyLookup)
2985	if dlg:
2986	dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d\nAll data Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
2987	if len(Hess):
2988	Vec += np.sum(dMdvhnp.sqrt(W[xB:xF])dy[xB:xF],axis=1)
2989	Hess += np.inner(dMdvh,dMdvh)
2990	else:
2991	Vec = np.sum(dMdvhnp.sqrt(W[xB:xF])dy[xB:xF],axis=1)
2992	Hess = np.inner(dMdvh,dMdvh)
2993	elif 'HKLF' in histogram[:4]:
2994	Histogram = Histograms[histogram]
2995	nobs = Histogram['Nobs']
2996	phase = Histogram['Reflection Lists']
2997	Phase = Phases[phase]
2998	hId = Histogram['hId']
2999	hfx = ':%d:'%(hId)
3000	wtFactor = calcControls[hfx+'wtFactor']
3001	pfx = '%d::'%(Phase['pId'])
3002	phfx = '%d:%d:'%(Phase['pId'],hId)
3003	SGData = Phase['General']['SGData']
3004	A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3005	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3006	refList = Histogram['Data']
3007	dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3008	dMdvh = np.zeros((len(varylist),len(refList)))
3009	wdf = np.zeros(len(refList))
3010	for iref,ref in enumerate(refList):
3011	if ref[6] > 0:
3012	if calcControls['F**2']:
3013	if ref[5]/ref[6] >= calcControls['minF/sig']:
3014	wdf[iref] = (ref[5]-ref[7])/ref[6]
3015	for j,var in enumerate(varylist):
3016	if var in dFdvDict:
3017	dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
3018	if phfx+'Scale' in varylist:
3019	dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]
3020	else:
3021	Fo = np.sqrt(ref[5])
3022	Fc = np.sqrt(ref[7])
3023	sig = ref[6]/(2.0*Fo)
3024	wdf[iref] = (Fo-Fc)/sig
3025	if Fo/sig >= calcControls['minF/sig']:
3026	for j,var in enumerate(varylist):
3027	if var in dFdvDict:
3028	dMdvh[j][iref] = dFdvDict[var][iref]/ref[6]
3029	if phfx+'Scale' in varylist:
3030	dMdvh[varylist.index(phfx+'Scale')][iref] = ref[9]/ref[6]
3031	if dlg:
3032	dlg.Update(Histogram['wR'],newmsg='Hessian for histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3033	if len(Hess):
3034	Vec += np.sum(dMdvh*wdf,axis=1)
3035	Hess += np.inner(dMdvh,dMdvh)
3036	else:
3037	Vec = np.sum(dMdvh*wdf,axis=1)
3038	Hess = np.inner(dMdvh,dMdvh)
3039	else:
3040	continue #skip non-histogram entries
3041	# dpdv = penaltyDeriv(parmdict,varylist)
3042	# if np.any(dpdv):
3043	# pVec = dpdv*penaltyFxn(parmdict,varylist)
3044	# Vec += pVec
3045	# pHess = np.zeros((len(varylist),len(varylist)))
3046	# for i,val in enumerate(dpdv):
3047	# pHess[i][i] = dpdv[i]**2
3048	# Hess += pHess
3049	return Vec,Hess
3050
3051	def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
3052	parmdict.update(zip(varylist,values))
3053	Values2Dict(parmdict, varylist, values)
3054	G2mv.Dict2Map(parmdict,varylist)
3055	Histograms,Phases,restraintDict = HistoPhases
3056	M = np.empty(0)
3057	SumwYo = 0
3058	Nobs = 0
3059	histoList = Histograms.keys()
3060	histoList.sort()
3061	for histogram in histoList:
3062	if 'PWDR' in histogram[:4]:
3063	Histogram = Histograms[histogram]
3064	hId = Histogram['hId']
3065	hfx = ':%d:'%(hId)
3066	wtFactor = calcControls[hfx+'wtFactor']
3067	Limits = calcControls[hfx+'Limits']
3068	x,y,w,yc,yb,yd = Histogram['Data']
3069	W = wtFactor*w
3070	yc *= 0.0 #zero full calcd profiles
3071	yb *= 0.0
3072	yd *= 0.0
3073	xB = np.searchsorted(x,Limits[0])
3074	xF = np.searchsorted(x,Limits[1])
3075	Histogram['Nobs'] = xF-xB
3076	Nobs += Histogram['Nobs']
3077	Histogram['sumwYo'] = np.sum(W[xB:xF]y[xB:xF]*2)
3078	SumwYo += Histogram['sumwYo']
3079	yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
3080	varylist,Histogram,Phases,calcControls,pawleyLookup)
3081	yc[xB:xF] += yb[xB:xF]
3082	yd[xB:xF] = y[xB:xF]-yc[xB:xF]
3083	Histogram['sumwYd'] = np.sum(np.sqrt(W[xB:xF])*(yd[xB:xF]))
3084	wdy = -np.sqrt(W[xB:xF])*(yd[xB:xF])
3085	Histogram['wR'] = min(100.,np.sqrt(np.sum(wdy*2)/Histogram['sumwYo'])100.)
3086	if dlg:
3087	dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3088	M = np.concatenate((M,wdy))
3089	#end of PWDR processing
3090	elif 'HKLF' in histogram[:4]:
3091	Histogram = Histograms[histogram]
3092	phase = Histogram['Reflection Lists']
3093	Phase = Phases[phase]
3094	hId = Histogram['hId']
3095	hfx = ':%d:'%(hId)
3096	wtFactor = calcControls[hfx+'wtFactor']
3097	pfx = '%d::'%(Phase['pId'])
3098	phfx = '%d:%d:'%(Phase['pId'],hId)
3099	SGData = Phase['General']['SGData']
3100	A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
3101	G,g = G2lat.A2Gmat(A) #recip & real metric tensors
3102	refList = Histogram['Data']
3103	refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmdict)
3104	df = np.zeros(len(refList))
3105	sumwYo = 0
3106	sumFo = 0
3107	sumFo2 = 0
3108	sumdF = 0
3109	sumdF2 = 0
3110	nobs = 0
3111	for i,ref in enumerate(refList):
3112	if ref[6] > 0:
3113	ref[7] = parmdict[phfx+'Scale']*ref[9]
3114	ref[8] = ref[5]/parmdict[phfx+'Scale']
3115	if calcControls['F**2']:
3116	if ref[5]/ref[6] >= calcControls['minF/sig']:
3117	sumFo2 += ref[5]
3118	Fo = np.sqrt(ref[5])
3119	sumFo += Fo
3120	sumFo2 += ref[5]
3121	sumdF += abs(Fo-np.sqrt(ref[7]))
3122	sumdF2 += abs(ref[5]-ref[7])
3123	nobs += 1
3124	df[i] = -(ref[5]-ref[7])/ref[6]
3125	sumwYo += (ref[5]/ref[6])**2
3126	else:
3127	Fo = np.sqrt(ref[5])
3128	Fc = np.sqrt(ref[7])
3129	sig = ref[6]/(2.0*Fo)
3130	if Fo/sig >= calcControls['minF/sig']:
3131	sumFo += Fo
3132	sumFo2 += ref[5]
3133	sumdF += abs(Fo-Fc)
3134	sumdF2 += abs(ref[5]-ref[7])
3135	nobs += 1
3136	df[i] = -(Fo-Fc)/sig
3137	sumwYo += (Fo/sig)**2
3138	Histogram['Nobs'] = nobs
3139	Histogram['sumwYo'] = sumwYo
3140	SumwYo += sumwYo
3141	Histogram['wR'] = min(100.,np.sqrt(np.sum(df*2)/Histogram['sumwYo'])100.)
3142	Histogram[phfx+'Rf'] = 100.*sumdF/sumFo
3143	Histogram[phfx+'Rf^2'] = 100.*sumdF2/sumFo2
3144	Histogram[phfx+'Nref'] = nobs
3145	Nobs += nobs
3146	if dlg:
3147	dlg.Update(Histogram['wR'],newmsg='For histogram %d Rw=%8.3f%s'%(hId,Histogram['wR'],'%'))[0]
3148	M = np.concatenate((M,df))
3149	# end of HKLF processing
3150	Histograms['sumwYo'] = SumwYo
3151	Histograms['Nobs'] = Nobs
3152	Rw = min(100.,np.sqrt(np.sum(M*2)/SumwYo)100.)
3153	if dlg:
3154	GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
3155	if not GoOn:
3156	parmDict['saved values'] = values
3157	raise Exception #Abort!!
3158	# pFunc = penaltyFxn(parmdict,varylist)
3159	# if np.any(pFunc):
3160	# M = np.concatenate((M,pFunc))
3161	return M
3162
3163	def Refine(GPXfile,dlg):
3164	import pytexture as ptx
3165	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
3166
3167	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
3168	ShowBanner(printFile)
3169	varyList = []
3170	parmDict = {}
3171	G2mv.InitVars()
3172	Controls = GetControls(GPXfile)
3173	ShowControls(Controls,printFile)
3174	calcControls = {}
3175	calcControls.update(Controls)
3176	constrDict,fixedList = GetConstraints(GPXfile)
3177	restraintDict = GetRestraints(GPXfile)
3178	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3179	if not Phases:
3180	print ' * ERROR - you have no phases! *'
3181	print ' * Refine aborted *'
3182	raise Exception
3183	if not Histograms:
3184	print ' * ERROR - you have no data to refine with! *'
3185	print ' * Refine aborted *'
3186	raise Exception
3187	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restraintDict,pFile=printFile)
3188	calcControls['atomIndx'] = atomIndx
3189	calcControls['Natoms'] = Natoms
3190	calcControls['FFtables'] = FFtables
3191	calcControls['BLtables'] = BLtables
3192	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
3193	calcControls.update(controlDict)
3194	histVary,histDict,controlDict = GetHistogramData(Histograms,pFile=printFile)
3195	calcControls.update(controlDict)
3196	varyList = phaseVary+hapVary+histVary
3197	parmDict.update(phaseDict)
3198	parmDict.update(hapDict)
3199	parmDict.update(histDict)
3200	GetFprime(calcControls,Histograms)
3201	# do constraint processing
3202	try:
3203	groups,parmlist = G2mv.GroupConstraints(constrDict)
3204	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3205	except:
3206	print ' * ERROR - your constraints are internally inconsistent *'
3207	# traceback for debug
3208	#print 'varyList',varyList
3209	#print 'constrDict',constrDict
3210	#print 'fixedList',fixedList
3211	#import traceback
3212	#print traceback.format_exc()
3213	raise Exception(' * Refine aborted *')
3214	# # check to see which generated parameters are fully varied
3215	# msg = G2mv.SetVaryFlags(varyList)
3216	# if msg:
3217	# print ' * ERROR - you have not set the refine flags for constraints consistently! *'
3218	# print msg
3219	# raise Exception(' * Refine aborted *')
3220	#print G2mv.VarRemapShow(varyList)
3221	G2mv.Map2Dict(parmDict,varyList)
3222	Rvals = {}
3223	while True:
3224	begin = time.time()
3225	values = np.array(Dict2Values(parmDict, varyList))
3226	Ftol = Controls['min dM/M']
3227	Factor = Controls['shift factor']
3228	maxCyc = Controls['max cyc']
3229	if 'Jacobian' in Controls['deriv type']:
3230	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3231	ftol=Ftol,col_deriv=True,factor=Factor,
3232	args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3233	ncyc = int(result[2]['nfev']/2)
3234	elif 'Hessian' in Controls['deriv type']:
3235	result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3236	args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3237	ncyc = result[2]['num cyc']+1
3238	Rvals['lamMax'] = result[2]['lamMax']
3239	else: #'numeric'
3240	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3241	args=([Histograms,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3242	ncyc = int(result[2]['nfev']/len(varyList))
3243	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
3244	# for item in table: print item,table[item] #useful debug - are things shifting?
3245	runtime = time.time()-begin
3246	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3247	Values2Dict(parmDict, varyList, result[0])
3248	G2mv.Dict2Map(parmDict,varyList)
3249
3250	Rvals['Nobs'] = Histograms['Nobs']
3251	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
3252	Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
3253	print >>printFile,'\n Refinement results:'
3254	print >>printFile,135*'-'
3255	print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
3256	print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3257	print >>printFile,' wR = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3258	try:
3259	covMatrix = result[1]*Rvals['GOF']
3260	sig = np.sqrt(np.diag(covMatrix))
3261	if np.any(np.isnan(sig)):
3262	print '* Least squares aborted - some invalid esds possible *'
3263	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
3264	# for item in table: print item,table[item] #useful debug - are things shifting?
3265	break #refinement succeeded - finish up!
3266	except TypeError: #result[1] is None on singular matrix
3267	print '** Refinement failed - singular matrix **'
3268	if 'Hessian' in Controls['deriv type']:
3269	num = len(varyList)-1
3270	for i,val in enumerate(np.flipud(result[2]['psing'])):
3271	if val:
3272	print 'Removing parameter: ',varyList[num-i]
3273	del(varyList[num-i])
3274	else:
3275	Ipvt = result[2]['ipvt']
3276	for i,ipvt in enumerate(Ipvt):
3277	if not np.sum(result[2]['fjac'],axis=1)[i]:
3278	print 'Removing parameter: ',varyList[ipvt-1]
3279	del(varyList[ipvt-1])
3280	break
3281
3282	# print 'dependentParmList: ',G2mv.dependentParmList
3283	# print 'arrayList: ',G2mv.arrayList
3284	# print 'invarrayList: ',G2mv.invarrayList
3285	# print 'indParmList: ',G2mv.indParmList
3286	# print 'fixedDict: ',G2mv.fixedDict
3287	# print 'test1'
3288	GetFobsSq(Histograms,Phases,parmDict,calcControls)
3289	# print 'test2'
3290	sigDict = dict(zip(varyList,sig))
3291	newCellDict = GetNewCellParms(parmDict,varyList)
3292	newAtomDict = ApplyXYZshifts(parmDict,varyList)
3293	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3294	'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3295	# add the uncertainties into the esd dictionary (sigDict)
3296	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3297	SetPhaseData(parmDict,sigDict,Phases,covData,printFile)
3298	SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
3299	SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
3300	G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
3301	SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3302	printFile.close()
3303	print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
3304	print ' *** Refinement successful ***'
3305
3306	#for testing purposes!!!
3307	if DEBUG:
3308	import cPickle
3309	fl = open('structTestdata.dat','wb')
3310	cPickle.dump(parmDict,fl,1)
3311	cPickle.dump(varyList,fl,1)
3312	for histogram in Histograms:
3313	if 'PWDR' in histogram[:4]:
3314	Histogram = Histograms[histogram]
3315	cPickle.dump(Histogram,fl,1)
3316	cPickle.dump(Phases,fl,1)
3317	cPickle.dump(calcControls,fl,1)
3318	cPickle.dump(pawleyLookup,fl,1)
3319	fl.close()
3320
3321	if dlg:
3322	return Rvals['Rwp']
3323
3324	def SeqRefine(GPXfile,dlg):
3325	import pytexture as ptx
3326	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
3327
3328	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
3329	print ' Sequential Refinement'
3330	ShowBanner(printFile)
3331	G2mv.InitVars()
3332	Controls = GetControls(GPXfile)
3333	ShowControls(Controls,printFile)
3334	constrDict,fixedList = GetConstraints(GPXfile)
3335	restraintDict = GetRestraints(GPXfile)
3336	Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3337	if not Phases:
3338	print ' * ERROR - you have no histograms to refine! *'
3339	print ' * Refine aborted *'
3340	raise Exception
3341	if not Histograms:
3342	print ' * ERROR - you have no data to refine with! *'
3343	print ' * Refine aborted *'
3344	raise Exception
3345	Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,restraintDict,False,printFile)
3346	if 'Seq Data' in Controls:
3347	histNames = Controls['Seq Data']
3348	else:
3349	histNames = GetHistogramNames(GPXfile,['PWDR',])
3350	if 'Reverse Seq' in Controls:
3351	if Controls['Reverse Seq']:
3352	histNames.reverse()
3353	SeqResult = {'histNames':histNames}
3354	makeBack = True
3355	for ihst,histogram in enumerate(histNames):
3356	ifPrint = False
3357	if dlg:
3358	dlg.SetTitle('Residual for histogram '+str(ihst))
3359	calcControls = {}
3360	calcControls['atomIndx'] = atomIndx
3361	calcControls['Natoms'] = Natoms
3362	calcControls['FFtables'] = FFtables
3363	calcControls['BLtables'] = BLtables
3364	varyList = []
3365	parmDict = {}
3366	Histo = {histogram:Histograms[histogram],}
3367	hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3368	calcControls.update(controlDict)
3369	histVary,histDict,controlDict = GetHistogramData(Histo,False)
3370	calcControls.update(controlDict)
3371	varyList = phaseVary+hapVary+histVary
3372	if not ihst:
3373	saveVaryList = varyList[:]
3374	for i,item in enumerate(saveVaryList):
3375	items = item.split(':')
3376	if items[1]:
3377	items[1] = ''
3378	item = ':'.join(items)
3379	saveVaryList[i] = item
3380	SeqResult['varyList'] = saveVaryList
3381	else:
3382	newVaryList = varyList[:]
3383	for i,item in enumerate(newVaryList):
3384	items = item.split(':')
3385	if items[1]:
3386	items[1] = ''
3387	item = ':'.join(items)
3388	newVaryList[i] = item
3389	if newVaryList != SeqResult['varyList']:
3390	print newVaryList
3391	print SeqResult['varyList']
3392	print '**** ERROR - variable list for this histogram does not match previous'
3393	raise Exception
3394	parmDict.update(phaseDict)
3395	parmDict.update(hapDict)
3396	parmDict.update(histDict)
3397	GetFprime(calcControls,Histo)
3398	# do constraint processing
3399	try:
3400	groups,parmlist = G2mv.GroupConstraints(constrDict)
3401	G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3402	except:
3403	print ' * ERROR - your constraints are internally inconsistent *'
3404	raise Exception(' * Refine aborted *')
3405	# check to see which generated parameters are fully varied
3406	# msg = G2mv.SetVaryFlags(varyList)
3407	# if msg:
3408	# print ' * ERROR - you have not set the refine flags for constraints consistently! *'
3409	# print msg
3410	# raise Exception(' * Refine aborted *')
3411	#print G2mv.VarRemapShow(varyList)
3412	G2mv.Map2Dict(parmDict,varyList)
3413	Rvals = {}
3414	while True:
3415	begin = time.time()
3416	values = np.array(Dict2Values(parmDict,varyList))
3417	Ftol = Controls['min dM/M']
3418	Factor = Controls['shift factor']
3419	maxCyc = Controls['max cyc']
3420
3421	if 'Jacobian' in Controls['deriv type']:
3422	result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3423	ftol=Ftol,col_deriv=True,factor=Factor,
3424	args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3425	ncyc = int(result[2]['nfev']/2)
3426	elif 'Hessian' in Controls['deriv type']:
3427	result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3428	args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3429	ncyc = result[2]['num cyc']+1
3430	else: #'numeric'
3431	result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3432	args=([Histo,Phases,restraintDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
3433	ncyc = int(result[2]['nfev']/len(varyList))
3434
3435
3436
3437	runtime = time.time()-begin
3438	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3439	Values2Dict(parmDict, varyList, result[0])
3440	G2mv.Dict2Map(parmDict,varyList)
3441
3442	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100. #to %
3443	Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
3444	Rvals['Nobs'] = Histo['Nobs']
3445	print >>printFile,'\n Refinement results for histogram: v'+histogram
3446	print >>printFile,135*'-'
3447	print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3448	print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3449	print >>printFile,' wRp = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3450	try:
3451	covMatrix = result[1]*Rvals['GOF']
3452	sig = np.sqrt(np.diag(covMatrix))
3453	if np.any(np.isnan(sig)):
3454	print '* Least squares aborted - some invalid esds possible *'
3455	ifPrint = True
3456	break #refinement succeeded - finish up!
3457	except TypeError: #result[1] is None on singular matrix
3458	print '** Refinement failed - singular matrix **'
3459	if 'Hessian' in Controls['deriv type']:
3460	num = len(varyList)-1
3461	for i,val in enumerate(np.flipud(result[2]['psing'])):
3462	if val:
3463	print 'Removing parameter: ',varyList[num-i]
3464	del(varyList[num-i])
3465	else:
3466	Ipvt = result[2]['ipvt']
3467	for i,ipvt in enumerate(Ipvt):
3468	if not np.sum(result[2]['fjac'],axis=1)[i]:
3469	print 'Removing parameter: ',varyList[ipvt-1]
3470	del(varyList[ipvt-1])
3471	break
3472
3473	GetFobsSq(Histo,Phases,parmDict,calcControls)
3474	sigDict = dict(zip(varyList,sig))
3475	newCellDict = GetNewCellParms(parmDict,varyList)
3476	newAtomDict = ApplyXYZshifts(parmDict,varyList)
3477	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3478	'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3479	SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
3480	SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
3481	SeqResult[histogram] = covData
3482	SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
3483	makeBack = False
3484	SetSeqResult(GPXfile,Histograms,SeqResult)
3485	printFile.close()
3486	print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
3487	print ' *** Sequential refinement successful ***'
3488
3489	def DistAngle(DisAglCtls,DisAglData):
3490	import numpy.ma as ma
3491
3492	def ShowBanner(name):
3493	print 80''
3494	print ' Interatomic Distances and Angles for phase '+name
3495	print 80'','\n'
3496
3497	ShowBanner(DisAglCtls['Name'])
3498	SGData = DisAglData['SGData']
3499	SGtext = G2spc.SGPrint(SGData)
3500	for line in SGtext: print line
3501	Cell = DisAglData['Cell']
3502
3503	Amat,Bmat = G2lat.cell2AB(Cell[:6])
3504	covData = {}
3505	if 'covData' in DisAglData:
3506	covData = DisAglData['covData']
3507	covMatrix = covData['covMatrix']
3508	varyList = covData['varyList']
3509	pfx = str(DisAglData['pId'])+'::'
3510	A = G2lat.cell2A(Cell[:6])
3511	cellSig = getCellEsd(pfx,SGData,A,covData)
3512	names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
3513	valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
3514	line = '\n Unit cell:'
3515	for name,vals in zip(names,valEsd):
3516	line += name+vals
3517	print line
3518	else:
3519	print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
3520	' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
3521	'%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
3522	Factor = DisAglCtls['Factors']
3523	Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
3524	indices = (-1,0,1)
3525	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3526	origAtoms = DisAglData['OrigAtoms']
3527	targAtoms = DisAglData['TargAtoms']
3528	for Oatom in origAtoms:
3529	OxyzNames = ''
3530	IndBlist = []
3531	Dist = []
3532	Vect = []
3533	VectA = []
3534	angles = []
3535	for Tatom in targAtoms:
3536	Xvcov = []
3537	TxyzNames = ''
3538	if 'covData' in DisAglData:
3539	OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
3540	TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
3541	Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
3542	result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
3543	BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
3544	AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
3545	for Txyz,Top,Tunit in result:
3546	Dx = (Txyz-np.array(Oatom[3:6]))+Units
3547	dx = np.inner(Amat,Dx)
3548	dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
3549	IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
3550	if np.any(IndB):
3551	for indb in IndB:
3552	for i in range(len(indb)):
3553	if str(dx.T[indb][i]) not in IndBlist:
3554	IndBlist.append(str(dx.T[indb][i]))
3555	unit = Units[indb][i]
3556	tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
3557	pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
3558	sig = 0.0
3559	if len(Xvcov):
3560	sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
3561	Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
3562	if (Dist[-1][-1]-AsumR) <= 0.:
3563	Vect.append(dx.T[indb][i]/Dist[-1][-2])
3564	VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
3565	else:
3566	Vect.append([0.,0.,0.])
3567	VectA.append([])
3568	Vect = np.array(Vect)
3569	angles = np.zeros((len(Vect),len(Vect)))
3570	angsig = np.zeros((len(Vect),len(Vect)))
3571	for i,veca in enumerate(Vect):
3572	if np.any(veca):
3573	for j,vecb in enumerate(Vect):
3574	if np.any(vecb):
3575	angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
3576	line = ''
3577	for i,x in enumerate(Oatom[3:6]):
3578	if len(Xvcov):
3579	line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
3580	else:
3581	line += '%12.5f'%(x)
3582	print '\n Distances & angles for ',Oatom[1],' at ',line
3583	print 80''
3584	line = ''
3585	for dist in Dist[:-1]:
3586	line += '%12s'%(dist[1].center(12))
3587	print ' To cell +(sym. op.) dist. ',line
3588	for i,dist in enumerate(Dist):
3589	line = ''
3590	for j,angle in enumerate(angles[i][0:i]):
3591	sig = angsig[i][j]
3592	if angle:
3593	if sig:
3594	line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
3595	else:
3596	val = '%.3f'%(angle)
3597	line += '%12s'%(val.center(12))
3598	else:
3599	line += 12*' '
3600	if dist[5]: #sig exists!
3601	val = G2mth.ValEsd(dist[4],dist[5])
3602	else:
3603	val = '%8.4f'%(dist[4])
3604	print ' %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
3605
3606	def DisAglTor(DATData):
3607	SGData = DATData['SGData']
3608	Cell = DATData['Cell']
3609
3610	Amat,Bmat = G2lat.cell2AB(Cell[:6])
3611	covData = {}
3612	pfx = ''
3613	if 'covData' in DATData:
3614	covData = DATData['covData']
3615	covMatrix = covData['covMatrix']
3616	varyList = covData['varyList']
3617	pfx = str(DATData['pId'])+'::'
3618	Datoms = []
3619	Oatoms = []
3620	for i,atom in enumerate(DATData['Datoms']):
3621	symop = atom[-1].split('+')
3622	if len(symop) == 1:
3623	symop.append('0,0,0')
3624	symop[0] = int(symop[0])
3625	symop[1] = eval(symop[1])
3626	atom.append(symop)
3627	Datoms.append(atom)
3628	oatom = DATData['Oatoms'][i]
3629	names = ['','','']
3630	if pfx:
3631	names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
3632	oatom += [names,]
3633	Oatoms.append(oatom)
3634	atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
3635	if DATData['Natoms'] == 4: #torsion
3636	Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3637	print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
3638	print ' NB: Atom sequence determined by selection order'
3639	return # done with torsion
3640	elif DATData['Natoms'] == 3: #angle
3641	Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3642	print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
3643	print ' NB: Atom sequence determined by selection order'
3644	else: #2 atoms - distance
3645	Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3646	print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
3647
3648	def BestPlane(PlaneData):
3649
3650	def ShowBanner(name):
3651	print 80''
3652	print ' Best plane result for phase '+name
3653	print 80'','\n'
3654
3655	ShowBanner(PlaneData['Name'])
3656
3657	Cell = PlaneData['Cell']
3658	Amat,Bmat = G2lat.cell2AB(Cell[:6])
3659	Atoms = PlaneData['Atoms']
3660	sumXYZ = np.zeros(3)
3661	XYZ = []
3662	Natoms = len(Atoms)
3663	for atom in Atoms:
3664	xyz = np.array(atom[3:6])
3665	XYZ.append(xyz)
3666	sumXYZ += xyz
3667	sumXYZ /= Natoms
3668	XYZ = np.array(XYZ)-sumXYZ
3669	XYZ = np.inner(Amat,XYZ).T
3670	Zmat = np.zeros((3,3))
3671	for i,xyz in enumerate(XYZ):
3672	Zmat += np.outer(xyz.T,xyz)
3673	print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
3674	Evec,Emat = nl.eig(Zmat)
3675	Evec = np.sqrt(Evec)/(Natoms-3)
3676	Order = np.argsort(Evec)
3677	XYZ = np.inner(XYZ,Emat.T).T
3678	XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
3679	print ' Atoms in Cartesian best plane coordinates:'
3680	print ' Name X Y Z'
3681	for i,xyz in enumerate(XYZ):
3682	print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
3683	print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])
3684
3685
3686	def main():
3687	arg = sys.argv
3688	if len(arg) > 1:
3689	GPXfile = arg[1]
3690	if not ospath.exists(GPXfile):
3691	print 'ERROR - ',GPXfile," doesn't exist!"
3692	exit()
3693	GPXpath = ospath.dirname(arg[1])
3694	Refine(GPXfile,None)
3695	else:
3696	print 'ERROR - missing filename'
3697	exit()
3698	print "Done"
3699
3700	if __name__ == '__main__':
3701	main()

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

Download in other formats: