Context Navigation

source: trunk/GSASIImath.py @ 2025

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

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

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

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

Download in other formats: