source: trunk/GSASIIpwd.py @ 5164

Last change on this file since 5164 was 5164, checked in by toby, 7 months ago

some PWD cleanups from Laue Fringe work & fullrmc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 198.9 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASIIpwd: Powder calculations module*
5==============================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2022-01-28 14:39:44 +0000 (Fri, 28 Jan 2022) $
10# $Author: toby $
11# $Revision: 5164 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 5164 2022-01-28 14:39:44Z toby $
14########### SVN repository information ###################
15from __future__ import division, print_function
16import sys
17import math
18import time
19import os
20import os.path
21import subprocess as subp
22import datetime as dt
23import copy
24
25import numpy as np
26import numpy.linalg as nl
27import numpy.ma as ma
28import random as rand
29import numpy.fft as fft
30import scipy.interpolate as si
31import scipy.stats as st
32import scipy.optimize as so
33import scipy.special as sp
34import scipy.signal as signal
35
36import GSASIIpath
37filversion = "$Revision: 5164 $"
38GSASIIpath.SetVersionNumber("$Revision: 5164 $")
39import GSASIIlattice as G2lat
40import GSASIIspc as G2spc
41import GSASIIElem as G2elem
42import GSASIImath as G2mth
43try:
44    import pypowder as pyd
45except ImportError:
46    print ('pypowder is not available - profile calcs. not allowed')
47try:
48    import pydiffax as pyx
49except ImportError:
50    print ('pydiffax is not available for this platform')
51import GSASIIfiles as G2fil
52
53   
54# trig functions in degrees
55tand = lambda x: math.tan(x*math.pi/180.)
56atand = lambda x: 180.*math.atan(x)/math.pi
57atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
58cosd = lambda x: math.cos(x*math.pi/180.)
59acosd = lambda x: 180.*math.acos(x)/math.pi
60rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
61#numpy versions
62npsind = lambda x: np.sin(x*np.pi/180.)
63npasind = lambda x: 180.*np.arcsin(x)/math.pi
64npcosd = lambda x: np.cos(x*math.pi/180.)
65npacosd = lambda x: 180.*np.arccos(x)/math.pi
66nptand = lambda x: np.tan(x*math.pi/180.)
67npatand = lambda x: 180.*np.arctan(x)/np.pi
68npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
69npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave    #=d*
70npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)    #=2pi*d*
71npq2T = lambda Q,wave: 2.0*npasind(0.25*Q*wave/np.pi)
72ateln2 = 8.0*math.log(2.0)
73sateln2 = np.sqrt(ateln2)
74nxs = np.newaxis
75
76#### Powder utilities ################################################################################
77def PhaseWtSum(G2frame,histo):
78    '''
79    Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases)
80   
81    :param G2frame: GSASII main frame structure
82    :param str histo: histogram name
83    :returns: sum(scale*mass) for phases in histo
84    '''
85    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
86    wtSum = 0.0
87    for phase in Phases:
88        if Phases[phase]['General']['Type'] != 'magnetic':
89            if histo in Phases[phase]['Histograms']:
90                if not Phases[phase]['Histograms'][histo]['Use']: continue
91                mass = Phases[phase]['General']['Mass']
92                phFr = Phases[phase]['Histograms'][histo]['Scale'][0]
93                wtSum += mass*phFr
94    return wtSum
95   
96#### GSASII pwdr & pdf calculation routines ################################################################################
97def Transmission(Geometry,Abs,Diam):
98    '''
99    Calculate sample transmission
100
101    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
102    :param float Abs: absorption coeff in cm-1
103    :param float Diam: sample thickness/diameter in mm
104    '''
105    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
106        MuR = Abs*Diam/20.0
107        if MuR <= 3.0:
108            T0 = 16/(3.*math.pi)
109            T1 = -0.045780
110            T2 = -0.02489
111            T3 = 0.003045
112            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
113            if T < -20.:
114                return 2.06e-9
115            else:
116                return math.exp(T)
117        else:
118            T1 = 1.433902
119            T2 = 0.013869+0.337894
120            T3 = 1.933433+1.163198
121            T4 = 0.044365-0.04259
122            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
123            return T/100.
124    elif 'plate' in Geometry:
125        MuR = Abs*Diam/10.
126        return math.exp(-MuR)
127    elif 'Bragg' in Geometry:
128        return 0.0
129       
130def SurfaceRough(SRA,SRB,Tth):
131    ''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction
132    :param float SRA: Suortti surface roughness parameter
133    :param float SRB: Suortti surface roughness parameter
134    :param float Tth: 2-theta(deg) - can be numpy array
135   
136    '''
137    sth = npsind(Tth/2.)
138    T1 = np.exp(-SRB/sth)
139    T2 = SRA+(1.-SRA)*np.exp(-SRB)
140    return (SRA+(1.-SRA)*T1)/T2
141   
142def SurfaceRoughDerv(SRA,SRB,Tth):
143    ''' Suortti surface roughness correction derivatives
144    :param float SRA: Suortti surface roughness parameter (dimensionless)
145    :param float SRB: Suortti surface roughness parameter (dimensionless)
146    :param float Tth: 2-theta(deg) - can be numpy array
147    :return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative
148    '''
149    sth = npsind(Tth/2.)
150    T1 = np.exp(-SRB/sth)
151    T2 = SRA+(1.-SRA)*np.exp(-SRB)
152    Trans = (SRA+(1.-SRA)*T1)/T2
153    dydSRA = ((1.-T1)*T2-(1.-np.exp(-SRB))*Trans)/T2**2
154    dydSRB = ((SRA-1.)*T1*T2/sth-Trans*(SRA-T2))/T2**2
155    return [dydSRA,dydSRB]
156
157def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
158    '''Calculate sample absorption
159    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
160    :param float MuR: absorption coeff * sample thickness/2 or radius
161    :param Tth: 2-theta scattering angle - can be numpy array
162    :param float Phi: flat plate tilt angle - future
163    :param float Psi: flat plate tilt axis - future
164    '''
165   
166    def muRunder3(MuR,Sth2):
167        T0 = 16.0/(3.*np.pi)
168        T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
169            0.109561*np.sqrt(Sth2)-26.04556
170        T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
171            1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
172        T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
173        Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
174        return np.exp(Trns)
175   
176    def muRover3(MuR,Sth2):
177        T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
178            10.02088*Sth2**3-3.36778*Sth2**4
179        T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
180            (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
181        T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
182            0.13576*np.sqrt(Sth2)+1.163198
183        T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
184        Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
185        return Trns/100.
186       
187    Sth2 = npsind(Tth/2.0)**2
188    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
189        if 'array' in str(type(MuR)):
190            MuRSTh2 = np.vstack((MuR,Sth2))
191            AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1]))
192            return AbsCr
193        else:
194            if MuR <= 3.0:
195                return muRunder3(MuR,Sth2)
196            else:
197                return muRover3(MuR,Sth2)
198    elif 'Bragg' in Geometry:
199        return 1.0
200    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
201        # and only defined for 2theta < 90
202        MuT = 2.*MuR
203        T1 = np.exp(-MuT)
204        T2 = np.exp(-MuT/npcosd(Tth))
205        Tb = MuT-MuT/npcosd(Tth)
206        return (T2-T1)/Tb
207    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
208        MuT = 2.*MuR
209        cth = npcosd(Tth/2.0)
210        return np.exp(-MuT/cth)/cth
211       
212def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
213    'needs a doc string'
214    dA = 0.001
215    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
216    if MuR:
217        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
218        return (AbsP-AbsM)/(2.0*dA)
219    else:
220        return (AbsP-1.)/dA
221       
222def Polarization(Pola,Tth,Azm=0.0):
223    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
224
225    :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
226    :param Azm: azimuthal angle e.g. 0.0 in plane of polarization - can be numpy array
227    :param Tth: 2-theta scattering angle - can be numpy array
228      which (if either) of these is "right"?
229    :return: (pola, dpdPola) - both 2-d arrays
230      * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
231        (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
232      * dpdPola: derivative needed for least squares
233
234    """
235    cazm = npcosd(Azm)**2
236    sazm = npsind(Azm)**2
237    pola = ((1.0-Pola)*cazm+Pola*sazm)*npcosd(Tth)**2+(1.0-Pola)*sazm+Pola*cazm
238    dpdPola = -npsind(Tth)**2*(sazm-cazm)
239    return pola,dpdPola
240   
241def Oblique(ObCoeff,Tth):
242    'currently assumes detector is normal to beam'
243    if ObCoeff:
244        K = (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
245        return K
246    else:
247        return 1.0
248               
249def Ruland(RulCoff,wave,Q,Compton):
250    'needs a doc string'
251    C = 2.9978e8
252    D = 1.5e-3
253    hmc = 0.024262734687    #Compton wavelength in A
254    sinth2 = (Q*wave/(4.0*np.pi))**2
255    dlam = (wave**2)*Compton*Q/C
256    dlam_c = 2.0*hmc*sinth2-D*wave**2
257    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
258
259def KleinNishina(wave,Q):
260    hmc = 0.024262734687    #Compton wavelength in A
261    TTh = npq2T(Q,wave)
262    P = 1./(1.+(1.-npcosd(TTh)*(hmc/wave)))
263    KN = (P**3-(P*npsind(TTh))**2+P)/(1.+npcosd(TTh)**2)
264    return KN
265   
266def LorchWeight(Q):
267    'needs a doc string'
268    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
269           
270def GetAsfMean(ElList,Sthl2):
271    '''Calculate various scattering factor terms for PDF calcs
272
273    :param dict ElList: element dictionary contains scattering factor coefficients, etc.
274    :param np.array Sthl2: numpy array of sin theta/lambda squared values
275    :returns: mean(f^2), mean(f)^2, mean(compton)
276    '''
277    sumNoAtoms = 0.0
278    FF = np.zeros_like(Sthl2)
279    FF2 = np.zeros_like(Sthl2)
280    CF = np.zeros_like(Sthl2)
281    for El in ElList:
282        sumNoAtoms += ElList[El]['FormulaNo']
283    for El in ElList:
284        el = ElList[El]
285        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
286        cf = G2elem.ComptonFac(el,Sthl2)
287        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
288        FF2 += ff2*el['FormulaNo']/sumNoAtoms
289        CF += cf*el['FormulaNo']/sumNoAtoms
290    return FF2,FF**2,CF
291   
292def GetNumDensity(ElList,Vol):
293    'needs a doc string'
294    sumNoAtoms = 0.0
295    for El in ElList:
296        sumNoAtoms += ElList[El]['FormulaNo']
297    return sumNoAtoms/Vol
298           
299def CalcPDF(data,inst,limits,xydata):
300    '''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
301    dict xydata; results are placed in xydata.
302    Calculation parameters are found in dicts data and inst and list limits.
303    The return value is at present an empty list.
304    '''
305    auxPlot = []
306    if 'T' in inst['Type'][0]:
307        Ibeg = 0
308        Ifin = len(xydata['Sample'][1][0])
309    else:
310        Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0])
311        Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1
312    #subtract backgrounds - if any & use PWDR limits
313    IofQ = copy.deepcopy(xydata['Sample'])
314    IofQ[1] = np.array([I[Ibeg:Ifin] for I in IofQ[1]])
315    if data['Sample Bkg.']['Name']:
316        IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult']
317    if data['Container']['Name']:
318        xycontainer = xydata['Container'][1][1]*data['Container']['Mult']
319        if data['Container Bkg.']['Name']:
320            xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult']
321        IofQ[1][1] += xycontainer[Ibeg:Ifin]
322    data['IofQmin'] = IofQ[1][1][-1]
323    IofQ[1][1] -= data.get('Flat Bkg',0.)
324    #get element data & absorption coeff.
325    ElList = data['ElList']
326    Tth = IofQ[1][0]    #2-theta or TOF!
327    if 'X' in inst['Type'][0]:
328        Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
329        #Apply angle dependent corrections
330        MuR = Abs*data['Diam']/20.0
331        IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
332        IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
333        if data['DetType'] == 'Area detector':
334            IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
335    elif 'T' in inst['Type'][0]:    #neutron TOF normalized data - needs wavelength dependent absorption
336        wave = 2.*G2lat.TOF2dsp(inst,IofQ[1][0])*npsind(inst['2-theta'][1]/2.)
337        Els = ElList.keys()
338        Isotope = {El:'Nat. abund.' for El in Els}
339        GD = {'AtomTypes':ElList,'Isotope':Isotope}
340        BLtables = G2elem.GetBLtable(GD)
341        FP,FPP = G2elem.BlenResTOF(Els,BLtables,wave)
342        Abs = np.zeros(len(wave))
343        for iel,El in enumerate(Els):
344            BL = BLtables[El][1]
345            SA = BL['SA']*wave/1.798197+4.0*np.pi*FPP[iel]**2 #+BL['SL'][1]?
346            SA *= ElList[El]['FormulaNo']/data['Form Vol']
347            Abs += SA
348        MuR = Abs*data['Diam']/2.
349        IofQ[1][1] /= Absorb(data['Geometry'],MuR,inst['2-theta'][1]*np.ones(len(wave))) 
350    # improves look of F(Q) but no impact on G(R)
351    # bBut,aBut = signal.butter(8,.5,"lowpass")
352    # IofQ[1][1] = signal.filtfilt(bBut,aBut,IofQ[1][1])
353    XY = IofQ[1]   
354    #convert to Q
355    nQpoints = 5000
356    if 'C' in inst['Type'][0]:
357        wave = G2mth.getWave(inst)
358        minQ = npT2q(Tth[0],wave)
359        maxQ = npT2q(Tth[-1],wave)   
360        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
361        dq = Qpoints[1]-Qpoints[0]
362        XY[0] = npT2q(XY[0],wave)
363        Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])    #interpolate I(Q)
364    elif 'T' in inst['Type'][0]:
365        difC = inst['difC'][1]
366        minQ = 2.*np.pi*difC/Tth[-1]
367        maxQ = 2.*np.pi*difC/Tth[0]
368        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
369        dq = Qpoints[1]-Qpoints[0]
370        XY[0] = 2.*np.pi*difC/XY[0]
371        Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][-1])    #interpolate I(Q)
372    Qdata -= np.min(Qdata)*data['BackRatio']
373   
374    qLimits = data['QScaleLim']
375    maxQ = np.searchsorted(Qpoints,min(Qpoints[-1],qLimits[1]))+1
376    minQ = np.searchsorted(Qpoints,min(qLimits[0],0.90*Qpoints[-1]))
377    qLimits = [Qpoints[minQ],Qpoints[maxQ-1]]
378    newdata = []
379    if len(IofQ) < 3:
380        xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],'']
381    else:
382        xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
383    for item in xydata['IofQ'][1]:
384        newdata.append(item[:maxQ])
385    xydata['IofQ'][1] = newdata
386   
387    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
388    if 'XC' in inst['Type'][0]:
389        FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
390    else: #TOF
391        CF = np.zeros(len(xydata['SofQ'][1][0]))
392        FFSq = np.ones(len(xydata['SofQ'][1][0]))
393        SqFF = np.ones(len(xydata['SofQ'][1][0]))
394    Q = xydata['SofQ'][1][0]
395#    auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
396    if 'XC' in inst['Type'][0]:
397#        CF *= KleinNishina(wave,Q)
398        ruland = Ruland(data['Ruland'],wave,Q,CF)
399#        auxPlot.append([Q,ruland,'Ruland'])     
400        CF *= ruland
401#    auxPlot.append([Q,CF,'CF-Corr'])
402    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
403    xydata['SofQ'][1][1] *= scale
404    if 'XC' in inst['Type'][0]:
405        xydata['SofQ'][1][1] -= CF
406    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
407    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
408    xydata['SofQ'][1][1] *= scale
409    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
410    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
411    if data['Lorch']:
412        xydata['FofQ'][1][1] *= LorchWeight(Q)   
413    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
414    xydata['gofr'] = copy.deepcopy(xydata['FofQ'])
415    nR = len(xydata['GofR'][1][1])
416    Rmax = GSASIIpath.GetConfigValue('PDF_Rmax',100.)
417    mul = int(round(2.*np.pi*nR/(Rmax*qLimits[1])))
418#    mul = int(round(2.*np.pi*nR/(data.get('Rmax',100.)*qLimits[1])))
419    R = 2.*np.pi*np.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1])
420    xydata['GofR'][1][0] = R
421    xydata['gofr'][1][0] = R
422    GR = -dq*np.imag(fft.fft(xydata['FofQ'][1][1],mul*nR)[:nR])*data.get('GR Scale',1.0)
423    xydata['GofR'][1][1] = GR
424    gr = GR/(np.pi*R)
425    xydata['gofr'][1][1] = gr
426    numbDen = 0.
427    if 'ElList' in data:
428        numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
429    if data.get('noRing',True):
430        Rmin = data['Rmin']
431        xydata['gofr'][1][1] = np.where(R<Rmin,-4.*numbDen,xydata['gofr'][1][1])
432        xydata['GofR'][1][1] = np.where(R<Rmin,-4.*R*np.pi*numbDen,xydata['GofR'][1][1])
433    return auxPlot
434   
435def PDFPeakFit(peaks,data):
436    rs2pi = 1./np.sqrt(2*np.pi)
437   
438    def MakeParms(peaks):
439        varyList = []
440        parmDict = {'slope':peaks['Background'][1][1]}
441        if peaks['Background'][2]:
442            varyList.append('slope')
443        for i,peak in enumerate(peaks['Peaks']):
444            parmDict['PDFpos;'+str(i)] = peak[0]
445            parmDict['PDFmag;'+str(i)] = peak[1]
446            parmDict['PDFsig;'+str(i)] = peak[2]
447            if 'P' in peak[3]:
448                varyList.append('PDFpos;'+str(i))
449            if 'M' in peak[3]:
450                varyList.append('PDFmag;'+str(i))
451            if 'S' in peak[3]:
452                varyList.append('PDFsig;'+str(i))
453        return parmDict,varyList
454       
455    def SetParms(peaks,parmDict,varyList):
456        if 'slope' in varyList:
457            peaks['Background'][1][1] = parmDict['slope']
458        for i,peak in enumerate(peaks['Peaks']):
459            if 'PDFpos;'+str(i) in varyList:
460                peak[0] = parmDict['PDFpos;'+str(i)]
461            if 'PDFmag;'+str(i) in varyList:
462                peak[1] = parmDict['PDFmag;'+str(i)]
463            if 'PDFsig;'+str(i) in varyList:
464                peak[2] = parmDict['PDFsig;'+str(i)]
465       
466   
467    def CalcPDFpeaks(parmdict,Xdata):
468        Z = parmDict['slope']*Xdata
469        ipeak = 0
470        while True:
471            try:
472                pos = parmdict['PDFpos;'+str(ipeak)]
473                mag = parmdict['PDFmag;'+str(ipeak)]
474                wid = parmdict['PDFsig;'+str(ipeak)]
475                wid2 = 2.*wid**2
476                Z += mag*rs2pi*np.exp(-(Xdata-pos)**2/wid2)/wid
477                ipeak += 1
478            except KeyError:        #no more peaks to process
479                return Z
480               
481    def errPDFProfile(values,xdata,ydata,parmdict,varylist):       
482        parmdict.update(zip(varylist,values))
483        M = CalcPDFpeaks(parmdict,xdata)-ydata
484        return M
485           
486    newpeaks = copy.copy(peaks)
487    iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0])
488    iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1
489    X = data[1][0][iBeg:iFin]
490    Y = data[1][1][iBeg:iFin]
491    parmDict,varyList = MakeParms(peaks)
492    if not len(varyList):
493        G2fil.G2Print (' Nothing varied')
494        return newpeaks,None,None,None,None,None
495   
496    Rvals = {}
497    values =  np.array(Dict2Values(parmDict, varyList))
498    result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001,
499           args=(X,Y,parmDict,varyList))
500    chisq = np.sum(result[2]['fvec']**2)
501    Values2Dict(parmDict, varyList, result[0])
502    SetParms(peaks,parmDict,varyList)
503    Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y**2))*100.      #to %
504    chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values))   #reduced chi^2 = M/(Nobs-Nvar)
505    sigList = list(np.sqrt(chisq*np.diag(result[1])))   
506    Z = CalcPDFpeaks(parmDict,X)
507    newpeaks['calc'] = [X,Z]
508    return newpeaks,result[0],varyList,sigList,parmDict,Rvals   
509   
510def MakeRDF(RDFcontrols,background,inst,pwddata):
511    auxPlot = []
512    if 'C' in inst['Type'][0] or 'B' in inst['Type'][0]:
513        Tth = pwddata[0]
514        wave = G2mth.getWave(inst)
515        minQ = npT2q(Tth[0],wave)
516        maxQ = npT2q(Tth[-1],wave)
517        powQ = npT2q(Tth,wave) 
518    elif 'T' in inst['Type'][0]:
519        TOF = pwddata[0]
520        difC = inst['difC'][1]
521        minQ = 2.*np.pi*difC/TOF[-1]
522        maxQ = 2.*np.pi*difC/TOF[0]
523        powQ = 2.*np.pi*difC/TOF
524    piDQ = np.pi/(maxQ-minQ)
525    Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
526    if RDFcontrols['UseObsCalc'] == 'obs-calc':
527        Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
528    elif RDFcontrols['UseObsCalc'] == 'obs-back':
529        Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
530    elif RDFcontrols['UseObsCalc'] == 'calc-back':
531        Qdata = si.griddata(powQ,pwddata[3]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
532    Qdata *= np.sin((Qpoints-minQ)*piDQ)/piDQ
533    Qdata *= 0.5*np.sqrt(Qpoints)       #Qbin normalization
534    dq = Qpoints[1]-Qpoints[0]
535    nR = len(Qdata)
536    R = 0.5*np.pi*np.linspace(0,nR,nR)/(4.*maxQ)
537    iFin = np.searchsorted(R,RDFcontrols['maxR'])+1
538    bBut,aBut = signal.butter(4,0.01)
539    Qsmooth = signal.filtfilt(bBut,aBut,Qdata)
540#    auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']])
541#    auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']])
542    DofR = dq*np.imag(fft.fft(Qsmooth,16*nR)[:nR])
543    auxPlot.append([R[:iFin],DofR[:iFin],'D(R) for '+RDFcontrols['UseObsCalc']])   
544    return auxPlot
545
546# PDF optimization =============================================================
547def OptimizePDF(data,xydata,limits,inst,showFit=True,maxCycles=25):
548    import scipy.optimize as opt
549    numbDen = GetNumDensity(data['ElList'],data['Form Vol'])
550    Min,Init,Done = SetupPDFEval(data,xydata,limits,inst,numbDen)
551    xstart = Init()
552    bakMul = data['Sample Bkg.']['Mult']
553    if showFit:
554        rms = Min(xstart)
555        G2fil.G2Print('  Optimizing corrections to improve G(r) at low r')
556        if data['Sample Bkg.'].get('Refine',False):
557#            data['Flat Bkg'] = 0.
558            G2fil.G2Print('  start: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f})'.format(
559                data['Ruland'],data['Sample Bkg.']['Mult'],rms))
560        else:
561            G2fil.G2Print('  start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.4f})'.format(
562                data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms))
563    if data['Sample Bkg.'].get('Refine',False):
564        res = opt.minimize(Min,xstart,bounds=([0.01,1.],[1.2*bakMul,0.8*bakMul]),
565                    method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
566    else:
567        res = opt.minimize(Min,xstart,bounds=([0.,None],[0,1],[0.01,1.]),
568                    method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001)
569    Done(res['x'])
570    if showFit:
571        if res['success']:
572            msg = 'Converged'
573        else:
574            msg = 'Not Converged'
575        if data['Sample Bkg.'].get('Refine',False):
576            G2fil.G2Print('  end:   Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f}) *** {} ***\n'.format(
577                data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
578        else:
579            G2fil.G2Print('  end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} RMS:{:.4f}) *** {} ***\n'.format(
580                data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
581    return res
582
583def SetupPDFEval(data,xydata,limits,inst,numbDen):
584    Data = copy.deepcopy(data)
585    BkgMax = 1.
586    def EvalLowPDF(arg):
587        '''Objective routine -- evaluates the RMS deviations in G(r)
588        from -4(pi)*#density*r for for r<Rmin
589        arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that
590        the min & max values are between 0 and 1.
591        '''
592        if Data['Sample Bkg.'].get('Refine',False):
593            R,S = arg
594            Data['Sample Bkg.']['Mult'] = S
595        else:
596            F,B,R = arg
597            Data['Flat Bkg'] = BkgMax*(2.*F-1.)
598            Data['BackRatio'] = B
599        Data['Ruland'] = R
600        CalcPDF(Data,inst,limits,xydata)
601        # test low r computation
602        g = xydata['GofR'][1][1]
603        r = xydata['GofR'][1][0]
604        g0 = g[r < Data['Rmin']] + 4*np.pi*r[r < Data['Rmin']]*numbDen
605        M = sum(g0**2)/len(g0)
606        return M
607    def GetCurrentVals():
608        '''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling
609        '''
610        if data['Sample Bkg.'].get('Refine',False):
611                return [max(data['Ruland'],.05),data['Sample']['Mult']]
612        try:
613            F = 0.5+0.5*data['Flat Bkg']/BkgMax
614        except:
615            F = 0
616        return [F,data['BackRatio'],max(data['Ruland'],.05)]
617    def SetFinalVals(arg):
618        '''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the
619        scaled, refined values and plot corrected region of G(r)
620        '''
621        if data['Sample Bkg.'].get('Refine',False):
622            R,S = arg
623            data['Sample Bkg.']['Mult'] = S
624        else:
625            F,B,R = arg
626            data['Flat Bkg'] = BkgMax*(2.*F-1.)
627            data['BackRatio'] = B
628        data['Ruland'] = R
629        CalcPDF(data,inst,limits,xydata)
630    EvalLowPDF(GetCurrentVals())
631    BkgMax = max(xydata['IofQ'][1][1])/50.
632    return EvalLowPDF,GetCurrentVals,SetFinalVals
633
634#### GSASII convolution peak fitting routines: Finger, Cox & Jephcoat model 
635def factorize(num):
636    ''' Provide prime number factors for integer num
637    :returns: dictionary of prime factors (keys) & power for each (data)
638    '''
639    factors = {}
640    orig = num
641
642    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
643    i, sqi = 2, 4
644    while sqi <= num:
645        while not num%i:
646            num /= i
647            factors[i] = factors.get(i, 0) + 1
648
649        sqi += 2*i + 1
650        i += 1
651
652    if num != 1 and num != orig:
653        factors[num] = factors.get(num, 0) + 1
654
655    if factors:
656        return factors
657    else:
658        return {num:1}          #a prime number!
659           
660def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
661    ''' Provide list of optimal data sizes for FFT calculations
662
663    :param int nmin: minimum data size >= 1
664    :param int nmax: maximum data size > nmin
665    :param int thresh: maximum prime factor allowed
666    :Returns: list of data sizes where the maximum prime factor is < thresh
667    ''' 
668    plist = []
669    nmin = max(1,nmin)
670    nmax = max(nmin+1,nmax)
671    for p in range(nmin,nmax):
672        if max(list(factorize(p).keys())) < thresh:
673            plist.append(p)
674    return plist
675
676np.seterr(divide='ignore')
677
678# Normal distribution
679
680# loc = mu, scale = std
681_norm_pdf_C = 1./math.sqrt(2*math.pi)
682class norm_gen(st.rv_continuous):
683    'needs a doc string'
684     
685    def pdf(self,x,*args,**kwds):
686        loc,scale=kwds['loc'],kwds['scale']
687        x = (x-loc)/scale
688        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
689       
690norm = norm_gen(name='norm',longname='A normal',extradoc="""
691
692Normal distribution
693
694The location (loc) keyword specifies the mean.
695The scale (scale) keyword specifies the standard deviation.
696
697normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
698""")
699
700## Cauchy
701
702# median = loc
703
704class cauchy_gen(st.rv_continuous):
705    'needs a doc string'
706
707    def pdf(self,x,*args,**kwds):
708        loc,scale=kwds['loc'],kwds['scale']
709        x = (x-loc)/scale
710        return 1.0/np.pi/(1.0+x*x) / scale
711       
712cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
713
714Cauchy distribution
715
716cauchy.pdf(x) = 1/(pi*(1+x**2))
717
718This is the t distribution with one degree of freedom.
719""")
720   
721
722class fcjde_gen(st.rv_continuous):
723    """
724    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
725    Ref: J. Appl. Cryst. (1994) 27, 892-900.
726
727    :param x: array -1 to 1
728    :param t: 2-theta position of peak
729    :param s: sum(S/L,H/L); S: sample height, H: detector opening,
730      L: sample to detector opening distance
731    :param dx: 2-theta step size in deg
732
733    :returns: for fcj.pdf
734
735     * T = x*dx+t
736     * s = S/L+H/L
737     * if x < 0::
738
739        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
740
741     * if x >= 0: fcj.pdf = 0   
742     
743    """
744    def _pdf(self,x,t,s,dx):
745        T = dx*x+t
746        ax2 = abs(npcosd(T))
747        ax = ax2**2
748        bx = npcosd(t)**2
749        bx = np.where(ax>bx,bx,ax)
750        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
751        fx = np.where(fx > 0.,fx,0.0)
752        return fx
753             
754    def pdf(self,x,*args,**kwds):
755        loc=kwds['loc']
756        return self._pdf(x-loc,*args)
757       
758fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
759               
760def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
761    '''Compute the Finger-Cox-Jepcoat modified Voigt function for a
762    CW powder peak by direct convolution. This version is not used.
763    '''
764    DX = xdata[1]-xdata[0]
765    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
766    x = np.linspace(pos-fmin,pos+fmin,256)
767    dx = x[1]-x[0]
768    Norm = norm.pdf(x,loc=pos,scale=widths[0])
769    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
770    arg = [pos,shl/57.2958,dx,]
771    FCJ = fcjde.pdf(x,*arg,loc=pos)
772    if len(np.nonzero(FCJ)[0])>5:
773        z = np.column_stack([Norm,Cauchy,FCJ]).T
774        Z = fft.fft(z)
775        Df = fft.ifft(Z.prod(axis=0)).real
776    else:
777        z = np.column_stack([Norm,Cauchy]).T
778        Z = fft.fft(z)
779        Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
780    Df /= np.sum(Df)
781    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
782    return intens*Df(xdata)*DX/dx
783   
784#### GSASII peak fitting routine: Finger, Cox & Jephcoat model       
785
786def getWidthsCW(pos,sig,gam,shl):
787    '''Compute the peak widths used for computing the range of a peak
788    for constant wavelength data. On low-angle side, 50 FWHM are used,
789    on high-angle side 75 are used, high angle side extended for axial divergence
790    (for peaks above 90 deg, these are reversed.)
791   
792    :param pos: peak position; 2-theta in degrees
793    :param sig: Gaussian peak variance in centideg^2
794    :param gam: Lorentzian peak width in centidegrees
795    :param shl: axial divergence parameter (S+H)/L
796   
797    :returns: widths; [Gaussian sigma, Lorentzian gamma] in degrees, and
798        low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position
799        reversed for 2-theta > 90 deg.
800    '''
801    widths = [np.sqrt(sig)/100.,gam/100.]
802    fwhm = 2.355*widths[0]+widths[1]
803    fmin = 50.*(fwhm+shl*abs(npcosd(pos)))
804    fmax = 75.0*fwhm
805    if pos > 90:
806        fmin,fmax = [fmax,fmin]         
807    return widths,fmin,fmax
808   
809def getWidthsED(pos,sig):
810    '''Compute the peak widths used for computing the range of a peak
811    for energy dispersive data. On low-energy side, 20 FWHM are used,
812    on high-energy side 20 are used
813   
814    :param pos: peak position; energy in keV (not used)
815    :param sig: Gaussian peak variance in keV^2
816   
817    :returns: widths; [Gaussian sigma] in keV, and
818        low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position
819    '''
820    widths = [np.sqrt(sig),.001]
821    fwhm = 2.355*widths[0]
822    fmin = 5.*fwhm
823    fmax = 5.*fwhm
824    return widths,fmin,fmax
825   
826def getWidthsTOF(pos,alp,bet,sig,gam):
827    '''Compute the peak widths used for computing the range of a peak
828    for constant wavelength data. 50 FWHM are used on both sides each
829    extended by exponential coeff.
830   
831    param pos: peak position; TOF in musec (not used)
832    param alp,bet: TOF peak exponential rise & decay parameters
833    param sig: Gaussian peak variance in musec^2
834    param gam: Lorentzian peak width in musec
835   
836    returns: widths; [Gaussian sigma, Lornetzian gamma] in musec
837    returns: low TOF, high TOF ends of peak; 50FWHM from position
838    '''
839    widths = [np.sqrt(sig),gam]
840    fwhm = 2.355*widths[0]+2.*widths[1]
841    fmin = 50.*fwhm*(1.+1./alp)   
842    fmax = 50.*fwhm*(1.+1./bet)
843    return widths,fmin,fmax
844   
845def getFWHM(pos,Inst):
846    '''Compute total FWHM from Thompson, Cox & Hastings (1987) , J. Appl. Cryst. 20, 79-83
847    via getgamFW(g,s).
848   
849    :param pos: float peak position in deg 2-theta or tof in musec
850    :param Inst: dict instrument parameters
851   
852    :returns float: total FWHM of pseudoVoigt in deg or musec
853    ''' 
854   
855    sig = lambda Th,U,V,W: np.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))
856    sigED = lambda E,A,B,C: np.sqrt(max(0.001,A*E**2+B*E+C))
857    sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1*dsp**2+S2*dsp**4+Sq*dsp)
858    gam = lambda Th,X,Y,Z: Z+X/cosd(Th)+Y*tand(Th)
859    gamTOF = lambda dsp,X,Y,Z: Z+X*dsp+Y*dsp**2
860    alpTOF = lambda dsp,alp: alp/dsp
861    betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp**4+betq/dsp**2
862    alpPink = lambda pos,alp0,alp1: alp0+alp1*tand(pos/2.)
863    betPink = lambda pos,bet0,bet1: bet0+bet1*tand(pos/2.)
864    if 'T' in Inst['Type'][0]:
865        dsp = pos/Inst['difC'][1]
866        alp = alpTOF(dsp,Inst['alpha'][1])
867        bet = betTOF(dsp,Inst['beta-0'][1],Inst['beta-1'][1],Inst['beta-q'][1])
868        s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
869        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
870        return getgamFW(g,s)+np.log(2.0)*(alp+bet)/(alp*bet)
871    elif 'C' in Inst['Type'][0]:
872        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
873        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
874        return getgamFW(g,s)/100.  #returns FWHM in deg
875    elif 'E' in Inst['Type'][0]:
876        s = sigED(pos,Inst['A'][1],Inst['B'][1],Inst['C'][1])
877        return 2.35482*s
878    else:   #'B'
879        alp = alpPink(pos,Inst['alpha-0'][1],Inst['alpha-1'][1])
880        bet = betPink(pos,Inst['beta-0'][1],Inst['beta-1'][1])
881        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])
882        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1],Inst['Z'][1])
883        return getgamFW(g,s)/100.+np.log(2.0)*(alp+bet)/(alp*bet)  #returns FWHM in deg
884   
885def getgamFW(g,s):
886    '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83
887    lambda fxn needs FWHM for both Gaussian & Lorentzian components
888   
889    :param g: float Lorentzian gamma = FWHM(L)
890    :param s: float Gaussian sig
891   
892    :returns float: total FWHM of pseudoVoigt
893    ''' 
894    gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
895    return gamFW(2.35482*s,g)   #sqrt(8ln2)*sig = FWHM(G)
896               
897def getBackground(pfx,parmDict,bakType,dataType,xdata,fixback=None):
898    '''Computes the background from vars pulled from gpx file or tree.
899    '''
900    if 'T' in dataType:
901        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
902    elif 'E' in dataType:
903        const = 4.*np.pi*npsind(parmDict[pfx+'2-theta']/2.0)
904        q = const*xdata
905    else:
906        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
907        q = npT2q(xdata,wave)
908    yb = np.zeros_like(xdata)
909    nBak = 0
910    sumBk = [0.,0.,0]
911    while True:
912        key = pfx+'Back;'+str(nBak)
913        if key in parmDict:
914            nBak += 1
915        else:
916            break
917#empirical functions
918    if bakType in ['chebyschev','cosine','chebyschev-1']:
919        dt = xdata[-1]-xdata[0]   
920        for iBak in range(nBak):
921            key = pfx+'Back;'+str(iBak)
922            if bakType == 'chebyschev':
923                ybi = parmDict[key]*(-1.+2.*(xdata-xdata[0])/dt)**iBak
924            elif bakType == 'chebyschev-1':
925                xpos = -1.+2.*(xdata-xdata[0])/dt
926                ybi = parmDict[key]*np.cos(iBak*np.arccos(xpos))
927            elif bakType == 'cosine':
928                ybi = parmDict[key]*npcosd(180.*xdata*iBak/xdata[-1])
929            yb += ybi
930        sumBk[0] = np.sum(yb)
931    elif bakType in ['Q^2 power series','Q^-2 power series']:
932        QT = 1.
933        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
934        for iBak in range(nBak-1):
935            key = pfx+'Back;'+str(iBak+1)
936            if '-2' in bakType:
937                QT *= (iBak+1)*q**-2
938            else:
939                QT *= q**2/(iBak+1)
940            yb += QT*parmDict[key]
941        sumBk[0] = np.sum(yb)
942    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
943        if nBak == 1:
944            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
945        elif nBak == 2:
946            dX = xdata[-1]-xdata[0]
947            T2 = (xdata-xdata[0])/dX
948            T1 = 1.0-T2
949            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
950        else:
951            xnomask = ma.getdata(xdata)
952            xmin,xmax = xnomask[0],xnomask[-1]
953            if bakType == 'lin interpolate':
954                bakPos = np.linspace(xmin,xmax,nBak,True)
955            elif bakType == 'inv interpolate':
956                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
957            elif bakType == 'log interpolate':
958                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
959            bakPos[0] = xmin
960            bakPos[-1] = xmax
961            bakVals = np.zeros(nBak)
962            for i in range(nBak):
963                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
964            bakInt = si.interp1d(bakPos,bakVals,'linear')
965            yb = bakInt(ma.getdata(xdata))
966        sumBk[0] = np.sum(yb)
967#Debye function
968    if pfx+'difC' in parmDict or 'E' in parmDict[pfx+'Type']:
969        ff = 1.
970    else:       
971        try:
972            wave = parmDict[pfx+'Lam']
973        except KeyError:
974            wave = parmDict[pfx+'Lam1']
975        SQ = (q/(4.*np.pi))**2
976        FF = G2elem.GetFormFactorCoeff('Si')[0]
977        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
978    iD = 0       
979    while True:
980        try:
981            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
982            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
983            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
984            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
985            yb += ybi
986            sumBk[1] += np.sum(ybi)
987            iD += 1       
988        except KeyError:
989            break
990#peaks
991    iD = 0
992    while True:
993        try:
994            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
995            pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1)
996            pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],0.01)
997            pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1)
998            if 'C' in dataType:
999                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1000            elif 'E' in dataType:
1001                Wd,fmin,fmax = getWidthsED(pkP,pkS)
1002            else: #'T'OF
1003                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1004            iBeg = np.searchsorted(xdata,pkP-fmin)
1005            iFin = np.searchsorted(xdata,pkP+fmax)
1006            lenX = len(xdata)
1007            if not iBeg:
1008                iFin = np.searchsorted(xdata,pkP+fmax)
1009            elif iBeg == lenX:
1010                iFin = iBeg
1011            else:
1012                iFin = np.searchsorted(xdata,pkP+fmax)
1013            if 'C' in dataType:
1014                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])[0]
1015                yb[iBeg:iFin] += ybi
1016            elif 'T' in dataType:
1017                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])[0]
1018                yb[iBeg:iFin] += ybi
1019            elif 'B' in dataType:
1020                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin])[0]
1021                yb[iBeg:iFin] += ybi
1022            elif 'E' in dataType:
1023                ybi = pkI*getPsVoigt(pkP,pkS*10.**4,pkG*100.,xdata[iBeg:iFin])[0]
1024                yb[iBeg:iFin] += ybi
1025            sumBk[2] += np.sum(ybi)
1026            iD += 1       
1027        except KeyError:
1028            break
1029        except ValueError:
1030            G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
1031            break
1032    if fixback is not None:   
1033        yb += parmDict[pfx+'BF mult']*fixback
1034        sumBk[0] = sum(yb)
1035    return yb,sumBk
1036   
1037def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata,fixback=None):
1038    'needs a doc string'
1039    if 'T' in dataType:
1040        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
1041    elif 'E' in dataType:
1042        const = 4.*np.pi*npsind(parmDict[hfx+'2-theta']/2.0)
1043        q = const*xdata
1044    else:
1045        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1046        q = 2.*np.pi*npsind(xdata/2.)/wave
1047    nBak = 0
1048    while True:
1049        key = hfx+'Back;'+str(nBak)
1050        if key in parmDict:
1051            nBak += 1
1052        else:
1053            break
1054    dydb = np.zeros(shape=(nBak,len(xdata)))
1055    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
1056    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
1057    dydfb = []
1058
1059    if bakType in ['chebyschev','cosine','chebyschev-1']:
1060        dt = xdata[-1]-xdata[0]   
1061        for iBak in range(nBak):   
1062            if bakType == 'chebyschev':
1063                dydb[iBak] = (-1.+2.*(xdata-xdata[0])/dt)**iBak
1064            elif bakType == 'chebyschev-1':
1065                xpos = -1.+2.*(xdata-xdata[0])/dt
1066                dydb[iBak] = np.cos(iBak*np.arccos(xpos))
1067            elif bakType == 'cosine':
1068                dydb[iBak] = npcosd(180.*xdata*iBak/xdata[-1])
1069    elif bakType in ['Q^2 power series','Q^-2 power series']:
1070        QT = 1.
1071        dydb[0] = np.ones_like(xdata)
1072        for iBak in range(nBak-1):
1073            if '-2' in bakType:
1074                QT *= (iBak+1)*q**-2
1075            else:
1076                QT *= q**2/(iBak+1)
1077            dydb[iBak+1] = QT
1078    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
1079        if nBak == 1:
1080            dydb[0] = np.ones_like(xdata)
1081        elif nBak == 2:
1082            dX = xdata[-1]-xdata[0]
1083            T2 = (xdata-xdata[0])/dX
1084            T1 = 1.0-T2
1085            dydb = [T1,T2]
1086        else:
1087            xnomask = ma.getdata(xdata)
1088            xmin,xmax = xnomask[0],xnomask[-1]
1089            if bakType == 'lin interpolate':
1090                bakPos = np.linspace(xmin,xmax,nBak,True)
1091            elif bakType == 'inv interpolate':
1092                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
1093            elif bakType == 'log interpolate':
1094                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
1095            bakPos[0] = xmin
1096            bakPos[-1] = xmax
1097            for i,pos in enumerate(bakPos):
1098                if i == 0:
1099                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
1100                elif i == len(bakPos)-1:
1101                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
1102                else:
1103                    dydb[i] = np.where(xdata>bakPos[i],
1104                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
1105                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
1106    if hfx+'difC' in parmDict:
1107        ff = 1.
1108    else:
1109        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1110        q = npT2q(xdata,wave)
1111        SQ = (q/(4*np.pi))**2
1112        FF = G2elem.GetFormFactorCoeff('Si')[0]
1113        ff = np.array(G2elem.ScatFac(FF,SQ)[0])*np.pi**2    #needs pi^2~10. for cw data (why?)
1114    iD = 0       
1115    while True:
1116        try:
1117            if hfx+'difC' in parmDict:
1118                q = 2*np.pi*parmDict[hfx+'difC']/xdata
1119            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
1120            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
1121            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
1122            sqr = np.sin(q*dbR)/(q*dbR)
1123            cqr = np.cos(q*dbR)
1124            temp = np.exp(-dbU*q**2)
1125            dyddb[3*iD] = ff*sqr*temp
1126            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(dbR)
1127            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2
1128            iD += 1
1129        except KeyError:
1130            break
1131    iD = 0
1132    while True:
1133        try:
1134            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
1135            pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1)
1136            pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],0.01)
1137            pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1)
1138            if 'C' in dataType:
1139                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1140            elif 'E' in dataType:
1141                Wd,fmin,fmax = getWidthsED(pkP,pkS)
1142            else: #'T' or 'B'
1143                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1144            iBeg = np.searchsorted(xdata,pkP-fmin)
1145            iFin = np.searchsorted(xdata,pkP+fmax)
1146            lenX = len(xdata)
1147            if not iBeg:
1148                iFin = np.searchsorted(xdata,pkP+fmax)
1149            elif iBeg == lenX:
1150                iFin = iBeg
1151            else:
1152                iFin = np.searchsorted(xdata,pkP+fmax)
1153            if 'C' in dataType:
1154                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
1155            elif 'E' in dataType:
1156                Df,dFdp,dFds,dFdg = getdPsVoigt(pkP,pkS*10.**4,pkG*100.,xdata[iBeg:iFin])
1157            else:   #'T'OF
1158                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
1159            dydpk[4*iD][iBeg:iFin] += pkI*dFdp
1160            dydpk[4*iD+1][iBeg:iFin] += Df
1161            dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
1162            dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
1163            iD += 1       
1164        except KeyError:
1165            break
1166        except ValueError:
1167            G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
1168            break       
1169    # fixed background from file
1170    if  fixback is not None:
1171        dydfb = fixback
1172    return dydb,dyddb,dydpk,dydfb
1173
1174#### Using old gsas fortran routines for powder peak shapes & derivatives
1175def getFCJVoigt3(pos,sig,gam,shl,xdata):
1176    '''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
1177    CW powder peak in external Fortran routine
1178   
1179    param pos: peak position in degrees
1180    param sig: Gaussian variance in centideg^2
1181    param gam: Lorentzian width in centideg
1182    param shl: axial divergence parameter (S+H)/L
1183    param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1184   
1185    returns: array: calculated peak function at each xdata
1186    returns: integral of peak; nominally = 1.0
1187    '''
1188    if len(xdata):
1189        cw = np.diff(xdata)
1190        cw = np.append(cw,cw[-1])
1191        Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1192        return Df,np.sum(100.*Df*cw)
1193    else:
1194        return 0.,1.
1195
1196def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1197    '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1198    function for a CW powder peak
1199   
1200    param pos: peak position in degrees
1201    param sig: Gaussian variance in centideg^2
1202    param gam: Lorentzian width in centideg
1203    param shl: axial divergence parameter (S+H)/L
1204    param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1205   
1206    returns: arrays: function and derivatives of pos, sig, gam, & shl
1207    '''
1208    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1209    return Df,dFdp,dFds,dFdg,dFdsh
1210
1211def getPsVoigt(pos,sig,gam,xdata):
1212    '''Compute the simple Pseudo-Voigt function for a
1213    small angle Bragg peak in external Fortran routine
1214   
1215    param pos: peak position in degrees
1216    param sig: Gaussian variance in centideg^2
1217    param gam: Lorentzian width in centideg
1218   
1219    returns: array: calculated peak function at each xdata
1220    returns: integral of peak; nominally = 1.0
1221    '''
1222   
1223    cw = np.diff(xdata)
1224    cw = np.append(cw,cw[-1])
1225    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
1226    return Df,np.sum(100.*Df*cw)
1227
1228def getdPsVoigt(pos,sig,gam,xdata):
1229    '''Compute the simple Pseudo-Voigt function derivatives for a
1230    small angle Bragg peak peak in external Fortran routine
1231   
1232    param pos: peak position in degrees
1233    param sig: Gaussian variance in centideg^2
1234    param gam: Lorentzian width in centideg
1235
1236    returns: arrays: function and derivatives of pos, sig & gam
1237    NB: the pos derivative has the opposite sign compared to that in other profile functions
1238    '''
1239   
1240    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1241    return Df,dFdp,dFds,dFdg
1242
1243def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1244    '''Compute the double exponential Pseudo-Voigt convolution function for a
1245    neutron TOF & CW pink peak in external Fortran routine
1246    '''
1247   
1248    cw = np.diff(xdata)
1249    cw = np.append(cw,cw[-1])
1250    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1251    return Df,np.sum(Df*cw) 
1252   
1253def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1254    '''Compute the double exponential Pseudo-Voigt convolution function derivatives for a
1255    neutron TOF & CW pink peak in external Fortran routine
1256    '''
1257   
1258    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1259    return Df,dFdp,dFda,dFdb,dFds,dFdg   
1260
1261def ellipseSize(H,Sij,GB):
1262    '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady
1263    '''
1264   
1265    HX = np.inner(H.T,GB)
1266    lenHX = np.sqrt(np.sum(HX**2))
1267    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
1268    R = np.inner(HX/lenHX,Rsize)**2*Esize         #want column length for hkl in crystal
1269    lenR = 1./np.sqrt(np.sum(R))
1270    return lenR
1271
1272def ellipseSizeDerv(H,Sij,GB):
1273    '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) derivative per note from Alexander Brady
1274    '''
1275   
1276    lenR = ellipseSize(H,Sij,GB)
1277    delt = 0.001
1278    dRdS = np.zeros(6)
1279    for i in range(6):
1280        Sij[i] -= delt
1281        lenM = ellipseSize(H,Sij,GB)
1282        Sij[i] += 2.*delt
1283        lenP = ellipseSize(H,Sij,GB)
1284        Sij[i] -= delt
1285        dRdS[i] = (lenP-lenM)/(2.*delt)
1286    return lenR,dRdS
1287
1288def getMustrain(HKL,G,SGData,muStrData):
1289    if muStrData[0] == 'isotropic':
1290        return np.ones(HKL.shape[1])*muStrData[1][0]
1291    elif muStrData[0] == 'uniaxial':
1292        H = np.array(HKL)
1293        P = np.array(muStrData[3])
1294        cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1295        Si = muStrData[1][0]
1296        Sa = muStrData[1][1]
1297        return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1298    else:       #generalized - P.W. Stephens model
1299        H = np.array(HKL)
1300        rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T])
1301        Strms = np.array(G2spc.MustrainCoeff(H,SGData))
1302        Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0)
1303        return np.sqrt(Sum)/rdsq
1304   
1305def getCrSize(HKL,G,GB,sizeData):
1306    if sizeData[0] == 'isotropic':
1307        return np.ones(HKL.shape[1])*sizeData[1][0]
1308    elif sizeData[0] == 'uniaxial':
1309        H = np.array(HKL)
1310        P = np.array(sizeData[3])
1311        cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1312        Si = sizeData[1][0]
1313        Sa = sizeData[1][1]
1314        return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1315    else:
1316        Sij =[sizeData[4][i] for i in range(6)]
1317        H = np.array(HKL)
1318        return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2
1319
1320def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1321    '''
1322    Generates allowed by symmetry reflections with d >= dmin
1323    NB: GenHKLf & checkMagextc return True for extinct reflections
1324
1325    :param dmin:  minimum d-spacing
1326    :param SGData: space group data obtained from SpcGroup
1327    :param A: lattice parameter terms A1-A6
1328    :param Inst: instrument parameter info
1329    :returns: HKLs: np.array hkl, etc for allowed reflections
1330
1331    '''
1332    HKL = G2lat.GenHLaue(dmin,SGData,A)       
1333    HKLs = []
1334    ds = []
1335    for h,k,l,d in HKL:
1336        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1337        if ext and 'MagSpGrp' in SGData:
1338            ext = G2spc.checkMagextc([h,k,l],SGData)
1339        if not ext:
1340            if nodup and int(10000*d) in ds:
1341                continue
1342            ds.append(int(10000*d))
1343            if Inst == None:
1344                HKLs.append([h,k,l,d,0,-1])
1345            else:
1346                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1347    return np.array(HKLs)
1348
1349def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1350    'needs a doc string'
1351    HKLs = []
1352    vec = np.array(Vec)
1353    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1354    dvec = 1./(maxH*vstar+1./dmin)
1355    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1356    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1357    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1358    ifMag = False
1359    if 'MagSpGrp' in SGData:
1360        ifMag = True
1361    for h,k,l,d in HKL:
1362        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1363        if not ext and d >= dmin:
1364            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1365        for dH in SSdH:
1366            if dH:
1367                DH = SSdH[dH]
1368                H = [h+DH[0],k+DH[1],l+DH[2]]
1369                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1370                if d >= dmin:
1371                    HKLM = np.array([h,k,l,dH])
1372                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1373                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1374    return G2lat.sortHKLd(HKLs,True,True,True)
1375
1376peakInstPrmMode = True
1377'''Determines the mode used for peak fitting. When peakInstPrmMode=True peak
1378width parameters are computed from the instrument parameters (UVW,... or
1379alpha,... etc) unless the individual parameter is refined. This allows the
1380instrument parameters to be refined. When peakInstPrmMode=False, the instrument
1381parameters are not used and cannot be refined.
1382The default is peakFitMode=True. This is changed only in
1383:func:`setPeakInstPrmMode`, which is called only from :mod:`GSASIIscriptable`.
1384'''
1385
1386def setPeakInstPrmMode(normal=True):
1387    '''Determines the mode used for peak fitting. If normal=True (default)
1388    peak width parameters are computed from the instrument parameters
1389    unless the individual parameter is refined. If normal=False,
1390    peak widths are used as supplied for each peak.
1391    At present this is called only in G2scriptable.
1392
1393    Note that normal=True unless this routine is called. Also,
1394    instrument parameters can only be refined with normal=True.
1395
1396    :param bool normal: setting to apply to global variable
1397      :data:`peakInstPrmMode`
1398    '''
1399    global peakInstPrmMode
1400    peakInstPrmMode = normal
1401
1402def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType):
1403    '''Computes the profiles from multiple peaks for individual peak fitting
1404    for powder patterns.
1405    NB: not used for Rietveld refinement
1406    '''
1407   
1408    yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
1409    yc = np.zeros_like(yb)
1410    # cw = np.diff(xdata)
1411    # cw = np.append(cw,cw[-1])
1412    if 'C' in dataType:
1413        shl = max(parmDict['SH/L'],0.002)
1414        Ka2 = False
1415        if 'Lam1' in parmDict.keys():
1416            Ka2 = True
1417            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1418            kRatio = parmDict['I(L2)/I(L1)']
1419        iPeak = 0
1420        while True:
1421            try:
1422                pos = parmDict['pos'+str(iPeak)]
1423                tth = (pos-parmDict['Zero'])
1424                intens = parmDict['int'+str(iPeak)]
1425                sigName = 'sig'+str(iPeak)
1426                if sigName in varyList or not peakInstPrmMode:
1427                    sig = parmDict[sigName]
1428                else:
1429                    sig = G2mth.getCWsig(parmDict,tth)
1430                sig = max(sig,0.001)          #avoid neg sigma^2
1431                gamName = 'gam'+str(iPeak)
1432                if gamName in varyList or not peakInstPrmMode:
1433                    gam = parmDict[gamName]
1434                else:
1435                    gam = G2mth.getCWgam(parmDict,tth)
1436                gam = max(gam,0.001)             #avoid neg gamma
1437                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1438                iBeg = np.searchsorted(xdata,pos-fmin)
1439                iFin = np.searchsorted(xdata,pos+fmin)
1440                if not iBeg+iFin:       #peak below low limit
1441                    iPeak += 1
1442                    continue
1443                elif not iBeg-iFin:     #peak above high limit
1444                    return yb+yc
1445                fp = getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])[0]
1446                yc[iBeg:iFin] += intens*fp
1447                if Ka2:
1448                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1449                    iBeg = np.searchsorted(xdata,pos2-fmin)
1450                    iFin = np.searchsorted(xdata,pos2+fmin)
1451                    if iBeg-iFin:
1452                        fp2 = getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])[0]
1453                        yc[iBeg:iFin] += intens*kRatio*fp2
1454                iPeak += 1
1455            except KeyError:        #no more peaks to process
1456                return yb+yc
1457    elif 'E' in dataType:
1458        iPeak = 0
1459        dsp = 1.0 #for now - fix later
1460        while True:
1461            try:
1462                pos = parmDict['pos'+str(iPeak)]
1463                intens = parmDict['int'+str(iPeak)]
1464                sigName = 'sig'+str(iPeak)
1465                if sigName in varyList or not peakInstPrmMode:
1466                    sig = parmDict[sigName]
1467                else:
1468                    sig = G2mth.getEDsig(parmDict,pos)
1469                sig = max(sig,0.001)          #avoid neg sigma^2
1470                Wd,fmin,fmax = getWidthsED(pos,sig)
1471                iBeg = np.searchsorted(xdata,pos-fmin)
1472                iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin))
1473                if not iBeg+iFin:       #peak below low limit
1474                    iPeak += 1
1475                    continue
1476                elif not iBeg-iFin:     #peak above high limit
1477                    return yb+yc
1478                yc[iBeg:iFin] += intens*getPsVoigt(pos,sig*10.**4,0.001,xdata[iBeg:iFin])[0]
1479                iPeak += 1
1480            except KeyError:        #no more peaks to process
1481                return yb+yc       
1482    elif 'B' in dataType:
1483        iPeak = 0
1484        dsp = 1.0 #for now - fix later
1485        while True:
1486            try:
1487                pos = parmDict['pos'+str(iPeak)]
1488                tth = (pos-parmDict['Zero'])
1489                intens = parmDict['int'+str(iPeak)]
1490                alpName = 'alp'+str(iPeak)
1491                if alpName in varyList or not peakInstPrmMode:
1492                    alp = parmDict[alpName]
1493                else:
1494                    alp = G2mth.getPinkalpha(parmDict,tth)
1495                alp = max(0.1,alp)
1496                betName = 'bet'+str(iPeak)
1497                if betName in varyList or not peakInstPrmMode:
1498                    bet = parmDict[betName]
1499                else:
1500                    bet = G2mth.getPinkbeta(parmDict,tth)
1501                bet = max(0.1,bet)
1502                sigName = 'sig'+str(iPeak)
1503                if sigName in varyList or not peakInstPrmMode:
1504                    sig = parmDict[sigName]
1505                else:
1506                    sig = G2mth.getCWsig(parmDict,tth)
1507                sig = max(sig,0.001)          #avoid neg sigma^2
1508                gamName = 'gam'+str(iPeak)
1509                if gamName in varyList or not peakInstPrmMode:
1510                    gam = parmDict[gamName]
1511                else:
1512                    gam = G2mth.getCWgam(parmDict,tth)
1513                gam = max(gam,0.001)             #avoid neg gamma
1514                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1515                iBeg = np.searchsorted(xdata,pos-fmin)
1516                iFin = np.searchsorted(xdata,pos+fmin)
1517                if not iBeg+iFin:       #peak below low limit
1518                    iPeak += 1
1519                    continue
1520                elif not iBeg-iFin:     #peak above high limit
1521                    return yb+yc
1522                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])[0]
1523                iPeak += 1
1524            except KeyError:        #no more peaks to process
1525                return yb+yc       
1526    else:
1527        Pdabc = parmDict['Pdabc']
1528        difC = parmDict['difC']
1529        iPeak = 0
1530        while True:
1531            try:
1532                pos = parmDict['pos'+str(iPeak)]               
1533                tof = pos-parmDict['Zero']
1534                dsp = tof/difC
1535                intens = parmDict['int'+str(iPeak)]
1536                alpName = 'alp'+str(iPeak)
1537                if alpName in varyList or not peakInstPrmMode:
1538                    alp = parmDict[alpName]
1539                else:
1540                    if len(Pdabc):
1541                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1542                    else:
1543                        alp = G2mth.getTOFalpha(parmDict,dsp)
1544                alp = max(0.1,alp)
1545                betName = 'bet'+str(iPeak)
1546                if betName in varyList or not peakInstPrmMode:
1547                    bet = parmDict[betName]
1548                else:
1549                    if len(Pdabc):
1550                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1551                    else:
1552                        bet = G2mth.getTOFbeta(parmDict,dsp)
1553                bet = max(0.0001,bet)
1554                sigName = 'sig'+str(iPeak)
1555                if sigName in varyList or not peakInstPrmMode:
1556                    sig = parmDict[sigName]
1557                else:
1558                    sig = G2mth.getTOFsig(parmDict,dsp)
1559                gamName = 'gam'+str(iPeak)
1560                if gamName in varyList or not peakInstPrmMode:
1561                    gam = parmDict[gamName]
1562                else:
1563                    gam = G2mth.getTOFgamma(parmDict,dsp)
1564                gam = max(gam,0.001)             #avoid neg gamma
1565                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1566                iBeg = np.searchsorted(xdata,pos-fmin)
1567                iFin = np.searchsorted(xdata,pos+fmax)
1568                lenX = len(xdata)
1569                if not iBeg:
1570                    iFin = np.searchsorted(xdata,pos+fmax)
1571                elif iBeg == lenX:
1572                    iFin = iBeg
1573                else:
1574                    iFin = np.searchsorted(xdata,pos+fmax)
1575                if not iBeg+iFin:       #peak below low limit
1576                    iPeak += 1
1577                    continue
1578                elif not iBeg-iFin:     #peak above high limit
1579                    return yb+yc
1580                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])[0]
1581                iPeak += 1
1582            except KeyError:        #no more peaks to process
1583                return yb+yc
1584           
1585def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType):
1586    '''Computes the profile derivatives for a powder pattern for single peak fitting
1587   
1588    return: np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1589   
1590    NB: not used for Rietveld refinement
1591    '''
1592    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1593    dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback)
1594    if 'Back;0' in varyList:            #background derivs are in front if present
1595        dMdv[0:len(dMdb)] = dMdb
1596    names = ['DebyeA','DebyeR','DebyeU']
1597    for name in varyList:
1598        if 'Debye' in name:
1599            parm,Id = name.split(';')
1600            ip = names.index(parm)
1601            dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1602    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1603    for name in varyList:
1604        if 'BkPk' in name:
1605            parm,Id = name.split(';')
1606            ip = names.index(parm)
1607            dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1608    cw = np.diff(xdata)
1609    cw = np.append(cw,cw[-1])
1610    if 'C' in dataType:
1611        shl = max(parmDict['SH/L'],0.002)
1612        Ka2 = False
1613        if 'Lam1' in parmDict.keys():
1614            Ka2 = True
1615            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1616            kRatio = parmDict['I(L2)/I(L1)']
1617        iPeak = 0
1618        while True:
1619            try:
1620                pos = parmDict['pos'+str(iPeak)]
1621                tth = (pos-parmDict['Zero'])
1622                intens = parmDict['int'+str(iPeak)]
1623                sigName = 'sig'+str(iPeak)
1624                if sigName in varyList or not peakInstPrmMode:
1625                    sig = parmDict[sigName]
1626                    dsdU = dsdV = dsdW = 0
1627                else:
1628                    sig = G2mth.getCWsig(parmDict,tth)
1629                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1630                sig = max(sig,0.001)          #avoid neg sigma
1631                gamName = 'gam'+str(iPeak)
1632                if gamName in varyList or not peakInstPrmMode:
1633                    gam = parmDict[gamName]
1634                    dgdX = dgdY = dgdZ = 0
1635                else:
1636                    gam = G2mth.getCWgam(parmDict,tth)
1637                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1638                gam = max(gam,0.001)             #avoid neg gamma
1639                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1640                iBeg = np.searchsorted(xdata,pos-fmin)
1641                iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin))
1642                if not iBeg+iFin:       #peak below low limit
1643                    iPeak += 1
1644                    continue
1645                elif not iBeg-iFin:     #peak above high limit
1646                    break
1647                dMdpk = np.zeros(shape=(6,len(xdata)))
1648                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1649                for i in range(1,5):
1650                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1651                dMdpk[0][iBeg:iFin] += dMdipk[0]
1652                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1653                if Ka2:
1654                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1655                    iBeg = np.searchsorted(xdata,pos2-fmin)
1656                    iFin = np.searchsorted(xdata,pos2+fmin)
1657                    if iBeg-iFin:
1658                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1659                        for i in range(1,5):
1660                            dMdpk[i][iBeg:iFin] += intens*kRatio*dMdipk2[i]
1661                        dMdpk[0][iBeg:iFin] += kRatio*dMdipk2[0]
1662                        dMdpk[5][iBeg:iFin] += dMdipk2[0]
1663                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1664                for parmName in ['pos','int','sig','gam']:
1665                    try:
1666                        idx = varyList.index(parmName+str(iPeak))
1667                        dMdv[idx] = dervDict[parmName]
1668                    except ValueError:
1669                        pass
1670                if 'U' in varyList:
1671                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1672                if 'V' in varyList:
1673                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1674                if 'W' in varyList:
1675                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1676                if 'X' in varyList:
1677                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1678                if 'Y' in varyList:
1679                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1680                if 'Z' in varyList:
1681                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1682                if 'SH/L' in varyList:
1683                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1684                if 'I(L2)/I(L1)' in varyList:
1685                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1686                iPeak += 1
1687            except KeyError:        #no more peaks to process
1688                break
1689    elif 'E' in dataType:
1690        iPeak = 0
1691        while True:
1692            try:
1693                pos = parmDict['pos'+str(iPeak)]
1694                intens = parmDict['int'+str(iPeak)]
1695                sigName = 'sig'+str(iPeak)
1696                if sigName in varyList or not peakInstPrmMode:
1697                    sig = parmDict[sigName]
1698                    dsdA = dsdB = dsdC = 0
1699                else:
1700                    sig = G2mth.getEDsig(parmDict,pos)
1701                    dsdA,dsdB,dsdC = G2mth.getEDsigDeriv(parmDict,pos)
1702                sig = max(sig,0.001)          #avoid neg sigma
1703                Wd,fmin,fmax = getWidthsED(pos,sig)
1704                iBeg = np.searchsorted(xdata,pos-fmin)
1705                iFin = np.searchsorted(xdata,pos+fmin)
1706                if not iBeg+iFin:       #peak below low limit
1707                    iPeak += 1
1708                    continue
1709                elif not iBeg-iFin:     #peak above high limit
1710                    break
1711                dMdpk = np.zeros(shape=(4,len(xdata)))
1712                dMdipk = getdPsVoigt(pos,sig*10.**4,0.001,xdata[iBeg:iFin])
1713                dMdpk[0][iBeg:iFin] += dMdipk[0]
1714                for i in range(1,4):
1715                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1716                dervDict = {'int':dMdpk[0],'pos':-dMdpk[1],'sig':dMdpk[2]*10**4}
1717                for parmName in ['pos','int','sig']:
1718                    try:
1719                        idx = varyList.index(parmName+str(iPeak))
1720                        dMdv[idx] = dervDict[parmName]
1721                    except ValueError:
1722                        pass
1723                if 'A' in varyList:
1724                    dMdv[varyList.index('A')] += dsdA*dervDict['sig']
1725                if 'B' in varyList:
1726                    dMdv[varyList.index('B')] += dsdB*dervDict['sig']
1727                if 'C' in varyList:
1728                    dMdv[varyList.index('C')] += dsdC*dervDict['sig']
1729                iPeak += 1
1730            except KeyError:        #no more peaks to process
1731                break       
1732       
1733    elif 'B' in dataType:
1734        iPeak = 0
1735        while True:
1736            try:
1737                pos = parmDict['pos'+str(iPeak)]
1738                tth = (pos-parmDict['Zero'])
1739                intens = parmDict['int'+str(iPeak)]
1740                alpName = 'alp'+str(iPeak)
1741                if alpName in varyList or not peakInstPrmMode:
1742                    alp = parmDict[alpName]
1743                    dada0 = dada1 = 0.0
1744                else:
1745                    alp = G2mth.getPinkalpha(parmDict,tth)
1746                    dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
1747                alp = max(0.0001,alp)
1748                betName = 'bet'+str(iPeak)
1749                if betName in varyList or not peakInstPrmMode:
1750                    bet = parmDict[betName]
1751                    dbdb0 = dbdb1 = 0.0
1752                else:
1753                    bet = G2mth.getPinkbeta(parmDict,tth)
1754                    dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
1755                bet = max(0.0001,bet)
1756                sigName = 'sig'+str(iPeak)
1757                if sigName in varyList or not peakInstPrmMode:
1758                    sig = parmDict[sigName]
1759                    dsdU = dsdV = dsdW = 0
1760                else:
1761                    sig = G2mth.getCWsig(parmDict,tth)
1762                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1763                sig = max(sig,0.001)          #avoid neg sigma
1764                gamName = 'gam'+str(iPeak)
1765                if gamName in varyList or not peakInstPrmMode:
1766                    gam = parmDict[gamName]
1767                    dgdX = dgdY = dgdZ = 0
1768                else:
1769                    gam = G2mth.getCWgam(parmDict,tth)
1770                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1771                gam = max(gam,0.001)             #avoid neg gamma
1772                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
1773                iBeg = np.searchsorted(xdata,pos-fmin)
1774                iFin = np.searchsorted(xdata,pos+fmin)
1775                if not iBeg+iFin:       #peak below low limit
1776                    iPeak += 1
1777                    continue
1778                elif not iBeg-iFin:     #peak above high limit
1779                    break
1780                dMdpk = np.zeros(shape=(7,len(xdata)))
1781                dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
1782                for i in range(1,6):
1783                    dMdpk[i][iBeg:iFin] += cw[iBeg:iFin]*intens*dMdipk[i]
1784                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1785                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
1786                for parmName in ['pos','int','alp','bet','sig','gam']:
1787                    try:
1788                        idx = varyList.index(parmName+str(iPeak))
1789                        dMdv[idx] = dervDict[parmName]
1790                    except ValueError:
1791                        pass
1792                if 'U' in varyList:
1793                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1794                if 'V' in varyList:
1795                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1796                if 'W' in varyList:
1797                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1798                if 'X' in varyList:
1799                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1800                if 'Y' in varyList:
1801                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1802                if 'Z' in varyList:
1803                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1804                if 'alpha-0' in varyList:
1805                    dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
1806                if 'alpha-1' in varyList:
1807                    dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
1808                if 'beta-0' in varyList:
1809                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1810                if 'beta-1' in varyList:
1811                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1812                iPeak += 1
1813            except KeyError:        #no more peaks to process
1814                break       
1815    else:
1816        Pdabc = parmDict['Pdabc']
1817        difC = parmDict['difC']
1818        iPeak = 0
1819        while True:
1820            try:
1821                pos = parmDict['pos'+str(iPeak)]               
1822                tof = pos-parmDict['Zero']
1823                dsp = tof/difC
1824                intens = parmDict['int'+str(iPeak)]
1825                alpName = 'alp'+str(iPeak)
1826                if alpName in varyList or not peakInstPrmMode:
1827                    alp = parmDict[alpName]
1828                else:
1829                    if len(Pdabc):
1830                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1831                        dada0 = 0
1832                    else:
1833                        alp = G2mth.getTOFalpha(parmDict,dsp)
1834                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1835                betName = 'bet'+str(iPeak)
1836                if betName in varyList or not peakInstPrmMode:
1837                    bet = parmDict[betName]
1838                else:
1839                    if len(Pdabc):
1840                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1841                        dbdb0 = dbdb1 = dbdb2 = 0
1842                    else:
1843                        bet = G2mth.getTOFbeta(parmDict,dsp)
1844                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1845                sigName = 'sig'+str(iPeak)
1846                if sigName in varyList or not peakInstPrmMode:
1847                    sig = parmDict[sigName]
1848                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1849                else:
1850                    sig = G2mth.getTOFsig(parmDict,dsp)
1851                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1852                gamName = 'gam'+str(iPeak)
1853                if gamName in varyList or not peakInstPrmMode:
1854                    gam = parmDict[gamName]
1855                    dsdX = dsdY = dsdZ = 0
1856                else:
1857                    gam = G2mth.getTOFgamma(parmDict,dsp)
1858                    dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1859                gam = max(gam,0.001)             #avoid neg gamma
1860                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1861                iBeg = np.searchsorted(xdata,pos-fmin)
1862                lenX = len(xdata)
1863                if not iBeg:
1864                    iFin = np.searchsorted(xdata,pos+fmax)
1865                elif iBeg == lenX:
1866                    iFin = iBeg
1867                else:
1868                    iFin = np.searchsorted(xdata,pos+fmax)
1869                if not iBeg+iFin:       #peak below low limit
1870                    iPeak += 1
1871                    continue
1872                elif not iBeg-iFin:     #peak above high limit
1873                    break
1874                dMdpk = np.zeros(shape=(7,len(xdata)))
1875                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1876                for i in range(1,6):
1877                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1878                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1879                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1880                for parmName in ['pos','int','alp','bet','sig','gam']:
1881                    try:
1882                        idx = varyList.index(parmName+str(iPeak))
1883                        dMdv[idx] = dervDict[parmName]
1884                    except ValueError:
1885                        pass
1886                if 'alpha' in varyList:
1887                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1888                if 'beta-0' in varyList:
1889                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1890                if 'beta-1' in varyList:
1891                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1892                if 'beta-q' in varyList:
1893                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1894                if 'sig-0' in varyList:
1895                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1896                if 'sig-1' in varyList:
1897                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1898                if 'sig-2' in varyList:
1899                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1900                if 'sig-q' in varyList:
1901                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1902                if 'X' in varyList:
1903                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1904                if 'Y' in varyList:
1905                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1906                if 'Z' in varyList:
1907                    dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1908                iPeak += 1
1909            except KeyError:        #no more peaks to process
1910                break
1911    if 'BF mult' in varyList:
1912        dMdv[varyList.index('BF mult')] = fixback
1913       
1914    return dMdv
1915       
1916def Dict2Values(parmdict, varylist):
1917    '''Use before call to leastsq to setup list of values for the parameters
1918    in parmdict, as selected by key in varylist'''
1919    return [parmdict[key] for key in varylist] 
1920   
1921def Values2Dict(parmdict, varylist, values):
1922    ''' Use after call to leastsq to update the parameter dictionary with
1923    values corresponding to keys in varylist'''
1924    parmdict.update(zip(varylist,values))
1925   
1926def SetBackgroundParms(Background):
1927    'Loads background parameters into dicts/lists to create varylist & parmdict'
1928    if len(Background) == 1:            # fix up old backgrounds
1929        Background.append({'nDebye':0,'debyeTerms':[]})
1930    bakType,bakFlag = Background[0][:2]
1931    backVals = Background[0][3:]
1932    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1933    Debye = Background[1]           #also has background peaks stuff
1934    backDict = dict(zip(backNames,backVals))
1935    backVary = []
1936    if bakFlag:
1937        backVary = backNames
1938
1939    backDict['nDebye'] = Debye['nDebye']
1940    debyeDict = {}
1941    debyeList = []
1942    for i in range(Debye['nDebye']):
1943        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1944        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1945        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1946    debyeVary = []
1947    for item in debyeList:
1948        if item[1]:
1949            debyeVary.append(item[0])
1950    backDict.update(debyeDict)
1951    backVary += debyeVary
1952
1953    backDict['nPeaks'] = Debye['nPeaks']
1954    peaksDict = {}
1955    peaksList = []
1956    for i in range(Debye['nPeaks']):
1957        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1958        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1959        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1960    peaksVary = []
1961    for item in peaksList:
1962        if item[1]:
1963            peaksVary.append(item[0])
1964    backDict.update(peaksDict)
1965    backVary += peaksVary
1966    if 'background PWDR' in Background[1]:
1967        backDict['Back File'] = Background[1]['background PWDR'][0]
1968        backDict['BF mult'] = Background[1]['background PWDR'][1]
1969        if len(Background[1]['background PWDR']) > 2:
1970            if Background[1]['background PWDR'][2]:
1971                backVary += ['BF mult',]
1972    return bakType,backDict,backVary
1973   
1974def DoCalibInst(IndexPeaks,Inst):
1975   
1976    def SetInstParms():
1977        dataType = Inst['Type'][0]
1978        insVary = []
1979        insNames = []
1980        insVals = []
1981        for parm in Inst:
1982            insNames.append(parm)
1983            insVals.append(Inst[parm][1])
1984            if parm in ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']:
1985                if Inst[parm][2]:
1986                    insVary.append(parm)
1987        instDict = dict(zip(insNames,insVals))
1988        return dataType,instDict,insVary
1989       
1990    def GetInstParms(parmDict,Inst,varyList):
1991        for name in Inst:
1992            Inst[name][1] = parmDict[name]
1993       
1994    def InstPrint(Inst,sigDict):
1995        print ('Instrument Parameters:')
1996        if 'C' in Inst['Type'][0] or 'B' in Inst['Type'][0]:
1997            ptfmt = "%12.6f"
1998        else:
1999            ptfmt = "%12.3f"
2000        ptlbls = 'names :'
2001        ptstr =  'values:'
2002        sigstr = 'esds  :'
2003        for parm in Inst:
2004            if parm in  ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']:
2005                ptlbls += "%s" % (parm.center(12))
2006                ptstr += ptfmt % (Inst[parm][1])
2007                if parm in sigDict:
2008                    sigstr += ptfmt % (sigDict[parm])
2009                else:
2010                    sigstr += 12*' '
2011        print (ptlbls)
2012        print (ptstr)
2013        print (sigstr)
2014       
2015    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
2016        parmDict.update(zip(varyList,values))
2017        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
2018
2019    peakPos = []
2020    peakDsp = []
2021    peakWt = []
2022    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
2023        if peak[2] and peak[3] and sig > 0.:
2024            peakPos.append(peak[0])
2025            peakDsp.append(peak[-1])    #d-calc
2026#            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
2027            peakWt.append(1./(sig*peak[-1]))   #
2028    peakPos = np.array(peakPos)
2029    peakDsp = np.array(peakDsp)
2030    peakWt = np.array(peakWt)
2031    dataType,insDict,insVary = SetInstParms()
2032    parmDict = {}
2033    parmDict.update(insDict)
2034    varyList = insVary
2035    if not len(varyList):
2036        G2fil.G2Print ('**** ERROR - nothing to refine! ****')
2037        return False
2038    while True:
2039        begin = time.time()
2040        values =  np.array(Dict2Values(parmDict, varyList))
2041        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
2042            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
2043        ncyc = int(result[2]['nfev']/2)
2044        runtime = time.time()-begin   
2045        chisq = np.sum(result[2]['fvec']**2)
2046        Values2Dict(parmDict, varyList, result[0])
2047        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
2048        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
2049        G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2050        G2fil.G2Print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF))
2051        try:
2052            sig = np.sqrt(np.diag(result[1])*GOF)
2053            if np.any(np.isnan(sig)):
2054                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
2055            break                   #refinement succeeded - finish up!
2056        except ValueError:          #result[1] is None on singular matrix
2057            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
2058       
2059    sigDict = dict(zip(varyList,sig))
2060    GetInstParms(parmDict,Inst,varyList)
2061    InstPrint(Inst,sigDict)
2062    return True
2063
2064def getHeaderInfo(FitPgm,dataType):
2065    '''Provide parameter name, label name and formatting information for the
2066    contents of the peak table and where used in DoPeakFit
2067    '''
2068    names = ['pos','int']
2069    lnames = ['position','intensity']
2070    if FitPgm == 'LaueFringe':
2071        names += ['sig','gam','cells']
2072        lnames += ['sigma','gamma','cells']
2073        fmt = ["%10.5f","%10.3f","%10.3f","%10.3f","%10.3f"]
2074    elif 'C' in dataType:
2075        names += ['sig','gam']
2076        lnames += ['sigma','gamma']
2077        fmt = ["%10.5f","%10.1f","%10.3f","%10.3f"]
2078    elif 'T' in dataType:
2079        names += ['alp','bet','sig','gam']
2080        lnames += ['alpha','beta','sigma','gamma']
2081        fmt = ["%10.2f","%10.4f","%8.3f","%8.5f","%10.3f","%10.3f"]
2082    elif 'E' in dataType:
2083        names += ['sig']
2084        lnames += ['sigma']
2085        fmt = ["%10.5f","%10.1f","%8.3f"]
2086    else: # 'B'
2087        names += ['alp','bet','sig','gam']
2088        lnames += ['alpha','beta','sigma','gamma']
2089        fmt = ["%10.5f","%10.1f","%8.2f","%8.4f","%10.3f","%10.3f"]
2090    return names, fmt, lnames
2091
2092def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,wtFactor=1.0,dlg=None,noFit=False):
2093    '''Called to perform a peak fit, refining the selected items in the peak
2094    table as well as selected items in the background.
2095
2096    :param str FitPgm: type of fit to perform. At present this is ignored.
2097    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
2098      four values followed by a refine flag where the values are: position, intensity,
2099      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
2100      tree entry, dict item "peaks"
2101    :param list Background: describes the background. List with two items.
2102      Item 0 specifies a background model and coefficients. Item 1 is a dict.
2103      From the Histogram/Background tree entry.
2104    :param list Limits: min and max x-value to use
2105    :param dict Inst: Instrument parameters
2106    :param dict Inst2: more Instrument parameters
2107    :param numpy.array data: a 5xn array. data[0] is the x-values,
2108      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
2109      calc, background and difference intensities, respectively.
2110    :param array fixback: fixed background array; same size as data[0-5]
2111    :param list prevVaryList: Used in sequential refinements to override the
2112      variable list. Defaults as an empty list.
2113    :param bool oneCycle: True if only one cycle of fitting should be performed
2114    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
2115      and derivType = controls['deriv type']. If None default values are used.
2116    :param float wtFactor: weight multiplier; = 1.0 by default
2117    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
2118      Defaults to None, which means no updates are done.
2119    :param bool noFit: When noFit is True, a refinement is not performed. Default
2120      is False.
2121
2122    '''
2123    def GetBackgroundParms(parmList,Background):
2124        iBak = 0
2125        while True:
2126            try:
2127                bakName = 'Back;'+str(iBak)
2128                Background[0][iBak+3] = parmList[bakName]
2129                iBak += 1
2130            except KeyError:
2131                break
2132        iDb = 0
2133        while True:
2134            names = ['DebyeA;','DebyeR;','DebyeU;']
2135            try:
2136                for i,name in enumerate(names):
2137                    val = parmList[name+str(iDb)]
2138                    Background[1]['debyeTerms'][iDb][2*i] = val
2139                iDb += 1
2140            except KeyError:
2141                break
2142        iDb = 0
2143        while True:
2144            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
2145            try:
2146                for i,name in enumerate(names):
2147                    val = parmList[name+str(iDb)]
2148                    Background[1]['peaksList'][iDb][2*i] = val
2149                iDb += 1
2150            except KeyError:
2151                break
2152        if 'BF mult' in parmList:
2153            Background[1]['background PWDR'][1] = parmList['BF mult']
2154               
2155    def BackgroundPrint(Background,sigDict):
2156        print ('Background coefficients for '+Background[0][0]+' function')
2157        ptfmt = "%12.5f"
2158        ptstr =  'value: '
2159        sigstr = 'esd  : '
2160        for i,back in enumerate(Background[0][3:]):
2161            ptstr += ptfmt % (back)
2162            if Background[0][1]:
2163                prm = 'Back;'+str(i)
2164                if prm in sigDict:
2165                    sigstr += ptfmt % (sigDict[prm])
2166                else:
2167                    sigstr += " "*12
2168            if len(ptstr) > 75:
2169                print (ptstr)
2170                if Background[0][1]: print (sigstr)
2171                ptstr =  'value: '
2172                sigstr = 'esd  : '
2173        if len(ptstr) > 8:
2174            print (ptstr)
2175            if Background[0][1]: print (sigstr)
2176
2177        if Background[1]['nDebye']:
2178            parms = ['DebyeA;','DebyeR;','DebyeU;']
2179            print ('Debye diffuse scattering coefficients')
2180            ptfmt = "%12.5f"
2181            print (' term       DebyeA       esd        DebyeR       esd        DebyeU        esd')
2182            for term in range(Background[1]['nDebye']):
2183                line = ' term %d'%(term)
2184                for ip,name in enumerate(parms):
2185                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
2186                    if name+str(term) in sigDict:
2187                        line += ptfmt%(sigDict[name+str(term)])
2188                    else:
2189                        line += " "*12
2190                print (line)
2191        if Background[1]['nPeaks']:
2192            print ('Coefficients for Background Peaks')
2193            ptfmt = "%15.3f"
2194            for j,pl in enumerate(Background[1]['peaksList']):
2195                names =  'peak %3d:'%(j+1)
2196                ptstr =  'values  :'
2197                sigstr = 'esds    :'
2198                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
2199                    val = pl[2*i]
2200                    prm = lbl+";"+str(j)
2201                    names += '%15s'%(prm)
2202                    ptstr += ptfmt%(val)
2203                    if prm in sigDict:
2204                        sigstr += ptfmt%(sigDict[prm])
2205                    else:
2206                        sigstr += " "*15
2207                print (names)
2208                print (ptstr)
2209                print (sigstr)
2210        if 'BF mult' in sigDict:
2211            print('Background file mult: %.3f(%d)'%(Background[1]['background PWDR'][1],int(1000*sigDict['BF mult'])))
2212                           
2213    def SetInstParms(Inst):
2214        dataType = Inst['Type'][0]
2215        insVary = []
2216        insNames = []
2217        insVals = []
2218        for parm in Inst:
2219            insNames.append(parm)
2220            insVals.append(Inst[parm][1])
2221            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C',
2222                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
2223                    insVary.append(parm)
2224        instDict = dict(zip(insNames,insVals))
2225        if 'SH/L' in instDict:
2226            instDict['SH/L'] = max(instDict['SH/L'],0.002)
2227        return dataType,instDict,insVary
2228       
2229    def GetInstParms(parmDict,Inst,varyList,Peaks):
2230        for name in Inst:
2231            Inst[name][1] = parmDict[name]
2232        iPeak = 0
2233        while True:
2234            try:
2235                sigName = 'sig'+str(iPeak)
2236                pos = parmDict['pos'+str(iPeak)]
2237                if sigName not in varyList and peakInstPrmMode:
2238                    if 'T' in Inst['Type'][0]:
2239                        dsp = G2lat.Pos2dsp(Inst,pos)
2240                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
2241                    if 'E' in Inst['Type'][0]:
2242                        parmDict[sigName] = G2mth.getEDsig(parmDict,pos)
2243                    else:
2244                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
2245                gamName = 'gam'+str(iPeak)
2246                if gamName not in varyList and peakInstPrmMode:
2247                    if 'T' in Inst['Type'][0]:
2248                        dsp = G2lat.Pos2dsp(Inst,pos)
2249                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
2250                    else:
2251                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
2252                iPeak += 1
2253            except KeyError:
2254                break
2255       
2256    def InstPrint(Inst,sigDict):
2257        print ('Instrument Parameters:')
2258        ptfmt = "%12.6f"
2259        ptlbls = 'names :'
2260        ptstr =  'values:'
2261        sigstr = 'esds  :'
2262        for parm in Inst:
2263            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C',
2264                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
2265                ptlbls += "%s" % (parm.center(12))
2266                ptstr += ptfmt % (Inst[parm][1])
2267                if parm in sigDict:
2268                    sigstr += ptfmt % (sigDict[parm])
2269                else:
2270                    sigstr += 12*' '
2271        print (ptlbls)
2272        print (ptstr)
2273        print (sigstr)
2274
2275    def SetPeaksParms(dataType,Peaks):
2276        peakNames = []
2277        peakVary = []
2278        peakVals = []
2279        names,_,_ = getHeaderInfo(FitPgm,dataType)
2280        for i,peak in enumerate(Peaks):
2281            for j,name in enumerate(names):
2282                peakVals.append(peak[2*j])
2283                parName = name+str(i)
2284                peakNames.append(parName)
2285                if peak[2*j+1]:
2286                    peakVary.append(parName)
2287        return dict(zip(peakNames,peakVals)),peakVary
2288               
2289    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
2290        names,_,_ = getHeaderInfo(FitPgm,Inst['Type'][0])
2291        for i,peak in enumerate(Peaks):
2292            pos = parmDict['pos'+str(i)]
2293            if 'difC' in Inst:
2294                dsp = pos/Inst['difC'][1]
2295            for j in range(len(names)):
2296                parName = names[j]+str(i)
2297                if parName in varyList or not peakInstPrmMode:
2298                    peak[2*j] = parmDict[parName]
2299                elif 'alp' in parName:
2300                    if 'T' in Inst['Type'][0]:
2301                        peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
2302                    else: #'B'
2303                        peak[2*j] = G2mth.getPinkalpha(parmDict,pos)
2304                elif 'bet' in parName:
2305                    if 'T' in Inst['Type'][0]:
2306                        peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
2307                    else:   #'B'
2308                        peak[2*j] = G2mth.getPinkbeta(parmDict,pos)
2309                elif 'sig' in parName:
2310                    if 'T' in Inst['Type'][0]:
2311                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
2312                    elif 'E' in Inst['Type'][0]:
2313                        peak[2*j] = G2mth.getEDsig(parmDict,pos)
2314                    else:   #'C' & 'B'
2315                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
2316                elif 'gam' in parName:
2317                    if 'T' in Inst['Type'][0]:
2318                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
2319                    else:   #'C' & 'B'
2320                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
2321                       
2322    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
2323        print ('Peak coefficients:')
2324        names,fmt,_ = getHeaderInfo(FitPgm,dataType)
2325        head = 13*' '
2326        for name in names:
2327            if name in ['alp','bet']:
2328                head += name.center(8)+'esd'.center(8)
2329            else:
2330                head += name.center(10)+'esd'.center(10)
2331        head += 'bins'.center(8)
2332        print (head)
2333        ptfmt = dict(zip(names,fmt))
2334        for i,peak in enumerate(Peaks):
2335            ptstr =  ':'
2336            for j in range(len(names)):
2337                name = names[j]
2338                parName = name+str(i)
2339                ptstr += ptfmt[name] % (parmDict[parName])
2340                if parName in varyList:
2341                    ptstr += ptfmt[name] % (sigDict[parName])
2342                else:
2343                    if name in ['alp','bet']:
2344                        ptstr += 8*' '
2345                    else:
2346                        ptstr += 10*' '
2347            ptstr += '%9.2f'%(ptsperFW[i])
2348            print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
2349               
2350    def devPeakProfile(values,xdata,ydata,fixback, weights,dataType,parmdict,varylist,bakType,dlg):
2351        '''Computes a matrix where each row is the derivative of the calc-obs
2352        values (see :func:`errPeakProfile`) with respect to each parameter
2353        in backVary,insVary,peakVary. Used for peak fitting.
2354        '''
2355        parmdict.update(zip(varylist,values))
2356        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,fixback,varylist,bakType)
2357           
2358    def errPeakProfile(values,xdata,ydata,fixback,weights,dataType,parmdict,varylist,bakType,dlg):
2359        '''Computes a vector with the weighted calc-obs values differences
2360        for peak fitting
2361        '''
2362        parmdict.update(zip(varylist,values))
2363        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,fixback,varylist,bakType)-ydata)
2364        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
2365        if dlg:
2366            dlg.Raise()
2367            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
2368            if not GoOn:
2369                return -M           #abort!!
2370        return M
2371
2372    #---- beginning of DoPeakFit ----------------------------------------------
2373    if controls:
2374        Ftol = controls['min dM/M']
2375    else:
2376        Ftol = 0.0001
2377    if oneCycle:
2378        Ftol = 1.0
2379    x,y,w,yc,yb,yd = data   #these are numpy arrays - remove masks!
2380    if fixback is None:
2381        fixback = np.zeros_like(y)
2382    yc *= 0.                            #set calcd ones to zero
2383    yb *= 0.
2384    yd *= 0.
2385    xBeg = np.searchsorted(x,Limits[0])
2386    xFin = np.searchsorted(x,Limits[1])+1
2387    bakType,bakDict,bakVary = SetBackgroundParms(Background)
2388    dataType,insDict,insVary = SetInstParms(Inst)
2389    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
2390    parmDict = {}
2391    parmDict.update(bakDict)
2392    parmDict.update(insDict)
2393    parmDict.update(peakDict)
2394    parmDict['Pdabc'] = []      #dummy Pdabc
2395    parmDict.update(Inst2)      #put in real one if there
2396    if prevVaryList:
2397        varyList = prevVaryList[:]
2398    else:
2399        varyList = bakVary+insVary+peakVary
2400    fullvaryList = varyList[:]
2401    if not peakInstPrmMode:
2402        for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1','A','B','C',
2403            'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
2404            if v in varyList:
2405                raise Exception('Instrumental profile terms cannot be varied '+
2406                                    'after setPeakInstPrmMode(False) is used')
2407    while not noFit:
2408        begin = time.time()
2409        values =  np.array(Dict2Values(parmDict, varyList))
2410        Rvals = {}
2411        badVary = []
2412        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
2413               args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
2414        ncyc = int(result[2]['nfev']/2)
2415        runtime = time.time()-begin   
2416        chisq = np.sum(result[2]['fvec']**2)
2417        Values2Dict(parmDict, varyList, result[0])
2418        Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactor*w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
2419        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
2420        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
2421        if ncyc:
2422            G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2423        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2424        sig = [0]*len(varyList)
2425        if len(varyList) == 0: break  # if nothing was refined
2426        try:
2427            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2428            if np.any(np.isnan(sig)):
2429                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
2430            break                   #refinement succeeded - finish up!
2431        except ValueError:          #result[1] is None on singular matrix
2432            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
2433            Ipvt = result[2]['ipvt']
2434            for i,ipvt in enumerate(Ipvt):
2435                if not np.sum(result[2]['fjac'],axis=1)[i]:
2436                    G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2437                    badVary.append(varyList[ipvt-1])
2438                    del(varyList[ipvt-1])
2439                    break
2440            else: # nothing removed
2441                break
2442    if dlg: dlg.Destroy()
2443    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin],fixback[xBeg:xFin])[0]
2444    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],fixback[xBeg:xFin],varyList,bakType)
2445    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2446    if noFit:
2447        return
2448    sigDict = dict(zip(varyList,sig))
2449    GetBackgroundParms(parmDict,Background)
2450    if bakVary: BackgroundPrint(Background,sigDict)
2451    GetInstParms(parmDict,Inst,varyList,Peaks)
2452    if insVary: InstPrint(Inst,sigDict)
2453    GetPeaksParms(Inst,parmDict,Peaks,varyList)
2454    binsperFWHM = []
2455    for peak in Peaks:
2456        FWHM = getFWHM(peak[0],Inst)
2457        try:
2458            xpk = x.searchsorted(peak[0])
2459            cw = x[xpk]-x[xpk-1]
2460            binsperFWHM.append(FWHM/cw)
2461        except IndexError:
2462            binsperFWHM.append(0.)
2463    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2464    if len(binsperFWHM):
2465        if min(binsperFWHM) < 1.:
2466            G2fil.G2Print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***')
2467            if 'T' in Inst['Type'][0]:
2468                G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2469            else:
2470                G2fil.G2Print (' Manually increase W in Instrument Parameters')
2471        elif min(binsperFWHM) < 4.:
2472            G2fil.G2Print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***')
2473            G2fil.G2Print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM)))
2474    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2475   
2476def calcIncident(Iparm,xdata):
2477    'needs a doc string'
2478
2479    def IfunAdv(Iparm,xdata):
2480        Itype = Iparm['Itype']
2481        Icoef = Iparm['Icoeff']
2482        DYI = np.ones((12,xdata.shape[0]))
2483        YI = np.ones_like(xdata)*Icoef[0]
2484       
2485        x = xdata/1000.                 #expressions are in ms
2486        if Itype == 'Exponential':
2487            for i in [1,3,5,7,9]:
2488                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2489                YI += Icoef[i]*Eterm
2490                DYI[i] *= Eterm
2491                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
2492        elif 'Maxwell'in Itype:
2493            Eterm = np.exp(-Icoef[2]/x**2)
2494            DYI[1] = Eterm/x**5
2495            DYI[2] = -Icoef[1]*DYI[1]/x**2
2496            YI += (Icoef[1]*Eterm/x**5)
2497            if 'Exponential' in Itype:
2498                for i in range(3,11,2):
2499                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2500                    YI += Icoef[i]*Eterm
2501                    DYI[i] *= Eterm
2502                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
2503            else:   #Chebyschev
2504                T = (2./x)-1.
2505                Ccof = np.ones((12,xdata.shape[0]))
2506                Ccof[1] = T
2507                for i in range(2,12):
2508                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
2509                for i in range(1,10):
2510                    YI += Ccof[i]*Icoef[i+2]
2511                    DYI[i+2] =Ccof[i]
2512        return YI,DYI
2513       
2514    Iesd = np.array(Iparm['Iesd'])
2515    Icovar = Iparm['Icovar']
2516    YI,DYI = IfunAdv(Iparm,xdata)
2517    YI = np.where(YI>0,YI,1.)
2518    WYI = np.zeros_like(xdata)
2519    vcov = np.zeros((12,12))
2520    k = 0
2521    for i in range(12):
2522        for j in range(i,12):
2523            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
2524            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
2525            k += 1
2526    M = np.inner(vcov,DYI.T)
2527    WYI = np.sum(M*DYI,axis=0)
2528    WYI = np.where(WYI>0.,WYI,0.)
2529    return YI,WYI
2530
2531#### RMCutilities ################################################################################
2532def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd):
2533    inst = PWDdata['Instrument Parameters'][0]
2534    Xsb = 0.
2535    Ysb = 0.
2536    if 'T' in inst['Type'][1]:
2537        difC = inst['difC'][1]
2538        if useSamBrd[0]:
2539            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2540                Xsb = 1.e-4*difC/Size[1][0]
2541        if useSamBrd[1]:
2542            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2543                Ysb = 1.e-6*difC*Mustrain[1][0]
2544        prms = ['Bank',
2545                'difC','difA','Zero','2-theta',
2546                'alpha','beta-0','beta-1',
2547                'sig-0','sig-1','sig-2',
2548                'Z','X','Y']
2549        fname = Name+'.inst'
2550        fl = open(fname,'w')
2551        fl.write('1\n')
2552        fl.write('%d\n'%int(inst[prms[0]][1]))
2553        fl.write('%19.11f%19.11f%19.11f%19.11f\n'%(inst[prms[1]][1],inst[prms[2]][1],inst[prms[3]][1],inst[prms[4]][1]))
2554        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2555        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))   
2556        fl.write('%12.6e%14.6e%14.6e%14.6e%14.6e\n'%(inst[prms[11]][1],inst[prms[12]][1]+Ysb,inst[prms[13]][1]+Xsb,0.0,0.0))
2557        fl.write('\n\n\n')
2558        fl.close()
2559    else:
2560        if useSamBrd[0]:
2561            wave = G2mth.getWave(inst)
2562            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2563                Xsb = 1.8*wave/(np.pi*Size[1][0])
2564        if useSamBrd[1]:
2565            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2566                Ysb = 0.0180*Mustrain[1][0]/np.pi
2567        prms = ['Bank',
2568                'Lam','Zero','Polariz.',
2569                'U','V','W',
2570                'X','Y']
2571        fname = Name+'.inst'
2572        fl = open(fname,'w')
2573        fl.write('1\n')
2574        fl.write('%d\n'%int(inst[prms[0]][1]))
2575        fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2576        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2577        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))   
2578        fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2579        fl.write('\n\n\n')
2580        fl.close()
2581    return fname
2582   
2583def MakeBack(PWDdata,Name):
2584    Back = PWDdata['Background'][0]
2585    inst = PWDdata['Instrument Parameters'][0]
2586    if 'chebyschev-1' != Back[0]:
2587        return None
2588    Nback = Back[2]
2589    BackVals = Back[3:]
2590    fname = Name+'.back'
2591    fl = open(fname,'w')
2592    fl.write('%10d\n'%Nback)
2593    for val in BackVals:
2594        if 'T' in inst['Type'][1]:
2595            fl.write('%12.6g\n'%(float(val)))
2596        else:
2597            fl.write('%12.6g\n'%val)
2598    fl.close()
2599    return fname
2600
2601def findDup(Atoms):
2602    Dup = []
2603    Fracs = []
2604    for iat1,at1 in enumerate(Atoms):
2605        if any([at1[0] in dup for dup in Dup]):
2606            continue
2607        else:
2608            Dup.append([at1[0],])
2609            Fracs.append([at1[6],])
2610        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2611            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2612                Dup[-1] += [at2[0],]
2613                Fracs[-1]+= [at2[6],]
2614    return Dup,Fracs
2615
2616def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):   
2617   
2618    Meta = RMCPdict['metadata']
2619    Atseq = RMCPdict['atSeq']
2620    Supercell =  RMCPdict['SuperCell']
2621    generalData = Phase['General']
2622    Dups,Fracs = findDup(Phase['Atoms'])
2623    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2624    ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2625    Sample = PWDdata['Sample Parameters']
2626    Meta['temperature'] = Sample['Temperature']
2627    Meta['pressure'] = Sample['Pressure']
2628    Cell = generalData['Cell'][1:7]
2629    Trans = np.eye(3)*np.array(Supercell)
2630    newPhase = copy.deepcopy(Phase)
2631    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2632    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2633    GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2634    RMCPdict['Rmax'] = np.min(np.sqrt(np.array([1./G2lat.calc_rDsq2(H,GB) for H in [[1,0,0],[0,1,0],[0,0,1]]])))/2.
2635    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2636    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2637    Natm = np.count_nonzero(Natm-1)
2638    Atoms = newPhase['Atoms']
2639    reset = False
2640   
2641    if ifSfracs:
2642        Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2643        Natm = np.count_nonzero(Natm-1)
2644        Satoms = []
2645        for i in range(len(Atoms)//Natm):
2646            ind = i*Natm
2647            Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2648        Natoms = []
2649        for satoms in Satoms:
2650            for idup,dup in enumerate(Dups):
2651                ldup = len(dup)
2652                natm = len(satoms)
2653                i = 0
2654                while i < natm:
2655                    if satoms[i][0] in dup:
2656                        atoms = satoms[i:i+ldup]
2657                        try:
2658                            atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2659                            Natoms.append(atom)
2660                        except IndexError:      #what about vacancies?
2661                            if 'Va' not in Atseq:
2662                                reset = True
2663                                Atseq.append('Va')
2664                                RMCPdict['aTypes']['Va'] = 0.0
2665                            atom = atoms[0]
2666                            atom[1] = 'Va'
2667                            Natoms.append(atom)
2668                        i += ldup
2669                    else:
2670                       i += 1
2671    else:
2672        Natoms = Atoms
2673   
2674    NAtype = np.zeros(len(Atseq))
2675    for atom in Natoms:
2676        NAtype[Atseq.index(atom[1])] += 1
2677    NAstr = ['%6d'%i for i in NAtype]
2678    Cell = newPhase['General']['Cell'][1:7]
2679    if os.path.exists(Name+'.his6f'):
2680        os.remove(Name+'.his6f')
2681    if os.path.exists(Name+'.neigh'):
2682        os.remove(Name+'.neigh')
2683    fname = Name+'.rmc6f'
2684    fl = open(fname,'w')
2685    fl.write('(Version 6f format configuration file)\n')
2686    for item in Meta:
2687        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2688    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2689    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2690    fl.write('Number of atoms:                %d\n'%len(Natoms))
2691    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2692    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2693            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2694    A,B = G2lat.cell2AB(Cell,True)
2695    fl.write('Lattice vectors (Ang):\n')   
2696    for i in [0,1,2]:
2697        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2698    fl.write('Atoms (fractional coordinates):\n')
2699    nat = 0
2700    for atm in Atseq:
2701        for iat,atom in enumerate(Natoms):
2702            if atom[1] == atm:
2703                nat += 1
2704                atcode = Atcodes[iat].split(':')
2705                cell = [0,0,0]
2706                if '+' in atcode[1]:
2707                    cell = eval(atcode[1].split('+')[1])
2708                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2709                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2710    fl.close()
2711    return fname,reset
2712
2713def MakeBragg(PWDdata,Name,Phase):
2714    generalData = Phase['General']
2715    Vol = generalData['Cell'][7]
2716    Data = PWDdata['Data']
2717    Inst = PWDdata['Instrument Parameters'][0]
2718    Bank = int(Inst['Bank'][1])
2719    Sample = PWDdata['Sample Parameters']
2720    Scale = Sample['Scale'][0]
2721    Limits = PWDdata['Limits'][1]
2722    Ibeg = np.searchsorted(Data[0],Limits[0])
2723    Ifin = np.searchsorted(Data[0],Limits[1])+1
2724    fname = Name+'.bragg'
2725    fl = open(fname,'w')
2726    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2727    if 'T' in Inst['Type'][0]:
2728        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2729        for i in range(Ibeg,Ifin-1):
2730            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2731    else:
2732        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2733        for i in range(Ibeg,Ifin-1):
2734            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2735    fl.close()
2736    return fname
2737
2738def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2739    Meta = RMCPdict['metadata']
2740    Times = RMCPdict['runTimes']
2741    Atseq = RMCPdict['atSeq']
2742    Natoms = RMCPdict['NoAtoms']
2743    sumatms = np.sum(np.array([Natoms[iatm] for iatm in Natoms]))
2744    Isotope = RMCPdict['Isotope']
2745    Isotopes = RMCPdict['Isotopes']
2746    Atypes = RMCPdict['aTypes']
2747    if 'Va' in Atypes:
2748        Isotope['Va'] = 'Nat. Abund.'
2749        Isotopes['Va'] = {'Nat. Abund.':{'SL':[0.0,0.0]}}
2750    atPairs = RMCPdict['Pairs']
2751    Files = RMCPdict['files']
2752    BraggWt = RMCPdict['histogram'][1]
2753    inst = PWDdata['Instrument Parameters'][0]
2754    try:
2755        refList = PWDdata['Reflection Lists'][Name]['RefList']
2756    except KeyError:
2757        return 'Error - missing reflection list; you must do Refine first'
2758    dMin = refList[-1][4]
2759    gsasType = 'xray2'
2760    if 'T' in inst['Type'][1]:
2761        gsasType = 'gsas3'
2762    elif 'X' in inst['Type'][1]:
2763        XFF = G2elem.GetFFtable(Atseq)
2764        Xfl = open(Name+'.xray','w')
2765        for atm in Atseq:
2766            fa = XFF[atm]['fa']
2767            fb = XFF[atm]['fb']
2768            fc = XFF[atm]['fc']
2769            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2770                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2771        Xfl.close()
2772    lenA = len(Atseq)
2773    Pairs = []
2774    Ncoeff = []
2775    Nblen = [Isotopes[at][Isotope[at]]['SL'][0] for at in Atypes]
2776    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2777        Pairs += pair
2778    for pair in Pairs:
2779        pair = pair.replace(' ','')
2780        at1,at2 = pair.split('-')
2781        if at1 == 'Va' or at2 == 'Va':
2782            ncoef = 0.0
2783        else:
2784            ncoef = Isotopes[at1][Isotope[at1]]['SL'][0]*Natoms[at1]/sumatms
2785            ncoef *= Isotopes[at2][Isotope[at2]]['SL'][0]*Natoms[at2]/sumatms
2786        if at1 != at2:
2787            ncoef *= 2.
2788        Ncoeff += [ncoef,]
2789    pairMin = [atPairs[pair] if pair in atPairs else [0.0,0.,0.] for pair in Pairs ]
2790    maxMoves = [Atypes[atm] if atm in Atypes else 0.0 for atm in Atseq ]
2791    fname = Name+'.dat'
2792    fl = open(fname,'w')
2793    fl.write(' %% Hand edit the following as needed\n')
2794    fl.write('TITLE :: '+Name+'\n')
2795    fl.write('MATERIAL :: '+Meta['material']+'\n')
2796    fl.write('PHASE :: '+Meta['phase']+'\n')
2797    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2798    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2799    if RMCPdict.get('useGPU',False):
2800        fl.write('GPU_ACCELERATOR :: 0\n')
2801    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2802    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2803    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2804    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2805    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2806    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2807    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2808    fl.write('PRINT_PERIOD :: 100\n')
2809    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2810    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2811    fl.write('\n')
2812    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2813    fl.write('\n')
2814    fl.write('FLAGS ::\n')
2815    fl.write('  > NO_MOVEOUT\n')
2816    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2817    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2818    fl.write('\n')
2819    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2820    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2821    fl.write('IGNORE_HISTORY_FILE ::\n')
2822    fl.write('\n')
2823    if 'T' in inst['Type'][1]:
2824        fl.write('NEUTRON_COEFFICIENTS :: '+''.join(['%9.5f'%coeff for coeff in Ncoeff])+'\n')
2825    fl.write('DISTANCE_WINDOW ::\n')
2826    fl.write('  > MNDIST :: %s\n'%minD)
2827    fl.write('  > MXDIST :: %s\n'%maxD)
2828    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2829        fl.write('\n')
2830        fl.write('POTENTIALS ::\n')
2831        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2832        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2833        if len(RMCPdict['Potentials']['Stretch']):
2834            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2835            for bond in RMCPdict['Potentials']['Stretch']:
2836                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2837        if len(RMCPdict['Potentials']['Angles']):
2838            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2839            for angle in RMCPdict['Potentials']['Angles']:
2840                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2841                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2842    if RMCPdict['useBVS']:
2843        fl.write('BVS ::\n')
2844        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2845        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2846        oxid = []
2847        for val in RMCPdict['Oxid']:
2848            if len(val) == 3:
2849                oxid.append(val[0][1:])
2850            else:
2851                oxid.append(val[0][2:])
2852        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2853        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2854        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2855        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))       
2856        fl.write('  > SAVE :: 100000\n')
2857        fl.write('  > UPDATE :: 100000\n')
2858        if len(RMCPdict['Swap']):
2859            fl.write('\n')
2860            fl.write('SWAP_MULTI ::\n')
2861            for swap in RMCPdict['Swap']:
2862                try:
2863                    at1 = Atseq.index(swap[0])
2864                    at2 = Atseq.index(swap[1])
2865                except ValueError:
2866                    break
2867                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2868       
2869    if len(RMCPdict['FxCN']):
2870        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2871        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2872            try:
2873                at1 = Atseq.index(fxcn[0])
2874                at2 = Atseq.index(fxcn[1])
2875            except ValueError:
2876                break
2877            fl.write('  > CSTR%d ::   %d %d %.2f %.2f %.2f %.2f %.6f\n'%(ifx+1,at1+1,at2+1,fxcn[2],fxcn[3],fxcn[4],fxcn[5],fxcn[6]))
2878    if len(RMCPdict['AveCN']):
2879        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2880        for iav,avcn in enumerate(RMCPdict['AveCN']):
2881            try:
2882                at1 = Atseq.index(avcn[0])
2883                at2 = Atseq.index(avcn[1])
2884            except ValueError:
2885                break
2886            fl.write('  > CAVSTR%d ::   %d %d %.2f %.2f %.2f %.6f\n'%(iav+1,at1+1,at2+1,avcn[2],avcn[3],avcn[4],avcn[5]))
2887    for File in Files:
2888        if Files[File][0] and Files[File][0] != 'Select':
2889            if 'Xray' in File and 'F(Q)' in File:
2890                fqdata = open(Files[File][0],'r')
2891                lines = int(fqdata.readline()[:-1])
2892                fqdata.close()
2893            fl.write('\n')
2894            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2895            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2896            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2897            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2898            if 'Xray' not in File:
2899                fl.write('  > START_POINT :: 1\n')
2900                fl.write('  > END_POINT :: 3000\n')
2901                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2902            fl.write('  > CONSTANT_OFFSET 0.000\n')
2903            fl.write('  > NO_FITTED_OFFSET\n')
2904            if RMCPdict['FitScale']:
2905                fl.write('  > FITTED_SCALE\n')
2906            else:
2907                fl.write('  > NO_FITTED_SCALE\n')
2908            if Files[File][3] !='RMC':
2909                fl.write('  > %s\n'%Files[File][3])
2910            if 'reciprocal' in File:
2911                fl.write('  > CONVOLVE ::\n')
2912                if 'Xray' in File:
2913                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2914                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2915                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2916                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2917    fl.write('\n')
2918    fl.write('BRAGG ::\n')
2919    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2920    fl.write('  > RECALCUATE\n')
2921    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2922    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2923    if 'T' in inst['Type'][1]:
2924        fl.write('  > SCATTERING LENGTH :: '+''.join(['%8.4f'%blen for blen in Nblen])+'\n')
2925    fl.write('\n')
2926    fl.write('END  ::\n')
2927    fl.close()
2928    return fname
2929
2930# def FindBonds(Phase,RMCPdict):
2931#     generalData = Phase['General']
2932#     cx,ct,cs,cia = generalData['AtomPtrs']
2933#     atomData = Phase['Atoms']
2934#     Res = 'RMC'
2935#     if 'macro' in generalData['Type']:
2936#         Res = atomData[0][ct-3]
2937#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2938#     Pairs = RMCPdict['Pairs']   #dict!
2939#     BondList = []
2940#     notNames = []
2941#     for FrstName in AtDict:
2942#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2943#         Atyp1 = AtDict[FrstName]
2944#         if 'Va' in Atyp1:
2945#             continue
2946#         for nbr in nbrs:
2947#             Atyp2 = AtDict[nbr[0]]
2948#             if 'Va' in Atyp2:
2949#                 continue
2950#             try:
2951#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2952#             except KeyError:
2953#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2954#             if any(bndData):
2955#                 if bndData[0] <= nbr[1] <= bndData[1]:
2956#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2957#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2958#                     if bondStr not in BondList and revbondStr not in BondList:
2959#                         BondList.append(bondStr)
2960#         notNames.append(FrstName)
2961#     return Res,BondList
2962
2963# def FindAngles(Phase,RMCPdict):
2964#     generalData = Phase['General']
2965#     Cell = generalData['Cell'][1:7]
2966#     Amat = G2lat.cell2AB(Cell)[0]
2967#     cx,ct,cs,cia = generalData['AtomPtrs']
2968#     atomData = Phase['Atoms']
2969#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2970#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2971#     Angles = RMCPdict['Angles']
2972#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2973#     AngleList = []
2974#     for MidName in AtDict:
2975#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2976#         if len(nbrs) < 2: #need 2 neighbors to make an angle
2977#             continue
2978#         Atyp2 = AtDict[MidName]
2979#         for i,nbr1 in enumerate(nbrs):
2980#             Atyp1 = AtDict[nbr1[0]]
2981#             for j,nbr3 in enumerate(nbrs[i+1:]):
2982#                 Atyp3 = AtDict[nbr3[0]]
2983#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2984#                 try:
2985#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2986#                 except KeyError:
2987#                     try:
2988#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2989#                     except KeyError:
2990#                         continue
2991#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2992#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
2993#                 if angData[0] <= calAngle <= angData[1]:
2994#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
2995#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
2996#                     if angStr not in AngleList and revangStr not in AngleList:
2997#                         AngleList.append(angStr)
2998#     return AngleList
2999
3000# def GetSqConvolution(XY,d):
3001
3002#     n = XY.shape[1]
3003#     snew = np.zeros(n)
3004#     dq = np.zeros(n)
3005#     sold = XY[1]
3006#     q = XY[0]
3007#     dq[1:] = np.diff(q)
3008#     dq[0] = dq[1]
3009   
3010#     for j in range(n):
3011#         for i in range(n):
3012#             b = abs(q[i]-q[j])
3013#             t = q[i]+q[j]
3014#             if j == i:
3015#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
3016#             else:
3017#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
3018#         snew[j] /= np.pi*q[j]
3019   
3020#     snew[0] = snew[1]
3021#     return snew
3022
3023# def GetMaxSphere(pdbName):
3024#     try:
3025#         pFil = open(pdbName,'r')
3026#     except FileNotFoundError:
3027#         return None
3028#     while True:
3029#         line = pFil.readline()
3030#         if 'Boundary' in line:
3031#             line = line.split()[3:]
3032#             G = np.array([float(item) for item in line])
3033#             G = np.reshape(G,(3,3))**2
3034#             G = nl.inv(G)
3035#             pFil.close()
3036#             break
3037#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
3038#     return min(dspaces)
3039
3040def findfullrmc():
3041    '''Find where fullrmc is installed. Tries the following:
3042   
3043         1. Returns the Config var 'fullrmc_exec', if defined. No check
3044            is done that the interpreter has fullrmc
3045         2. The current Python interpreter if fullrmc can be imported
3046            and fullrmc is version 5+
3047         3. The path is checked for a fullrmc image as named by Bachir
3048
3049    :returns: the full path to a python executable that is assumed to
3050      have fullrmc installed or None, if it was not found.
3051    '''
3052    is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
3053    if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
3054            GSASIIpath.GetConfigValue('fullrmc_exec')):
3055        return GSASIIpath.GetConfigValue('fullrmc_exec')
3056    try:
3057        import fullrmc
3058        if int(fullrmc.__version__.split('.')[0]) >= 5:
3059            return sys.executable
3060    except:
3061        pass
3062    pathlist = os.environ["PATH"].split(os.pathsep)
3063    for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
3064                  os.path.split(sys.executable)[0]):
3065        if p not in pathlist: pathlist.insert(0,p)
3066    import glob
3067    for p in pathlist:
3068        if sys.platform == "darwin":
3069            lookfor = "fullrmc*macOS*i386-64bit"
3070        elif sys.platform == "win32":
3071            lookfor = "fullrmc*.exe"
3072        else:
3073            lookfor = "fullrmc*"
3074        fl = glob.glob(lookfor)
3075        if len(fl) > 0:
3076            return os.path.abspath(sorted(fl)[0])
3077       
3078def findPDFfit():
3079    '''Find if PDFfit2 is installed (may be local to GSAS-II). Does the following:
3080    :returns: the full path to a python executable or None, if it was not found.
3081    '''
3082    try:
3083        sys.path.append('%s'%GSASIIpath.path2GSAS2) 
3084        from diffpy.pdffit2 import PdfFit
3085        return sys.executable
3086    except:
3087        return None
3088   
3089def GetPDFfitAtomVar(Phase,RMCPdict):
3090    ''' Find dict of independent "@n" variables for PDFfit in atom constraints
3091    '''
3092    General = Phase['General']
3093    Atoms = Phase['Atoms']
3094    cx,ct,cs,cia = General['AtomPtrs']
3095    AtomVar = RMCPdict['AtomVar']
3096    varnames = []
3097    for iat,atom in enumerate(RMCPdict['AtomConstr']):
3098        for it,item in enumerate(atom):
3099            if it > 1 and item:
3100                itms = item.split('@')
3101                for itm in itms[1:]:
3102                    itnum = itm[:2]
3103                    varname = '@%s'%itnum
3104                    varnames.append(varname)
3105                    if it < 6:
3106                        if varname not in AtomVar:
3107                            AtomVar[varname] = 0.0      #put ISODISTORT mode displ here?
3108                    else:
3109                        for i in range(3):
3110                            if varname not in AtomVar:
3111                                AtomVar[varname] = Atoms[iat][cia+i+2]
3112    varnames = set(varnames)
3113    for name in list(AtomVar.keys()):       #clear out unused parameters
3114        if name not in varnames:
3115            del AtomVar[name]
3116   
3117def MakePDFfitAtomsFile(Phase,RMCPdict):
3118    '''Make the PDFfit atoms file
3119    '''
3120    General = Phase['General']
3121    fName = General['Name']+'-PDFfit.stru'
3122    fName = fName.replace(' ','_')
3123    fatm = open(fName,'w')
3124    fatm.write('title  structure of '+General['Name']+'\n')
3125    fatm.write('format pdffit\n')
3126    fatm.write('scale   1.000000\n')    #fixed
3127    # if RMCPdict['shape'] == 'sphere':
3128    #     sharp = '%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0])
3129    # else:
3130    #     sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3131    sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3132    fatm.write('sharp '+sharp)
3133    shape = ''
3134    if RMCPdict['shape'] == 'sphere' and RMCPdict['spdiameter'][0] > 0.:
3135        shape = '   sphere, %10.6f\n'%RMCPdict['spdiameter'][0]
3136    elif RMCPdict['stepcut'] > 0.:
3137        shape = 'stepcut, %10.6f\n'%RMCPdict['stepcut']
3138    if shape:
3139        fatm.write('shape  '+shape)
3140    fatm.write('spcgr   %s\n'%General['SGData']['SpGrp'].replace(' ',''))
3141    cell = General['Cell'][1:7]
3142    fatm.write('cell  %10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n'%(
3143        cell[0],cell[1],cell[2],cell[3],cell[4],cell[5]))
3144    fatm.write('dcell '+5*'  0.000000,'+'  0.000000\n') 
3145    Atoms = Phase['Atoms']
3146    fatm.write('ncell %8d,%8d,%8d,%10d\n'%(1,1,1,len(Atoms)))
3147    fatm.write('atoms\n')
3148    cx,ct,cs,cia = General['AtomPtrs']
3149    for atom in Atoms:
3150        fatm.write('%4s%18.8f%18.8f%18.8f%13.4f\n'%(atom[ct][:2].ljust(2),atom[cx],atom[cx+1],atom[cx+2],atom[cx+3]))
3151        fatm.write('    '+'%18.8f%18.8f%18.8f%13.4f\n'%(0.,0.,0.,0.))
3152        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(atom[cia+2],atom[cia+3],atom[cia+4]))
3153        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.,))
3154        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(atom[cia+5],atom[cia+6],atom[cia+7]))
3155        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.))
3156    fatm.close()
3157   
3158def MakePDFfitRunFile(Phase,RMCPdict):
3159    '''Make the PDFfit python run file
3160    '''
3161   
3162    def GetCellConstr(SGData):
3163        if SGData['SGLaue'] in ['m3', 'm3m']:
3164            return [1,1,1,0,0,0]
3165        elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3166            return [1,1,2,0,0,0]
3167        elif SGData['SGLaue'] in ['3R','3mR']:
3168            return [1,1,1,2,2,2]
3169        elif SGData['SGLaue'] == 'mmm':
3170            return [1,2,3,0,0,0]
3171        elif SGData['SGLaue'] == '2/m':
3172            if SGData['SGUniq'] == 'a':
3173                return [1,2,3,4,0,0]
3174            elif SGData['SGUniq'] == 'b':
3175                return [1,2,3,0,4,0]
3176            elif SGData['SGUniq'] == 'c':
3177                return [1,2,3,0,0,4]
3178        else:
3179            return [1,2,3,4,5,6]
3180       
3181    General = Phase['General']
3182    Cell = General['Cell'][1:7]
3183    G2path = GSASIIpath.path2GSAS2
3184    rundata = '''
3185#!/usr/bin/env python
3186# -*- coding: utf-8 -*-
3187import sys
3188sys.path.append('%s')
3189from diffpy.pdffit2 import PdfFit
3190pf = PdfFit()
3191'''%G2path
3192    Nd = 0
3193    Np = 0
3194    parms = {}
3195    if 'sequential' in RMCPdict['refinement']:
3196        Np = 3
3197        rundata += '#sequential data here\n'
3198    else:
3199        for file in RMCPdict['files']:
3200            if 'Select' in RMCPdict['files'][file][0]:
3201                continue
3202            if 'Neutron' in file:
3203                Nd += 1
3204                dType = 'Ndata'
3205            else:
3206                Nd += 1
3207                dType = 'Xdata'
3208            rundata += "pf.read_data('%s', '%s', 30.0, %.4f)\n"%(RMCPdict['files'][file][0],dType[0],RMCPdict[dType]['qdamp'][0])
3209            rundata += 'pf.setdata(%d)\n'%Nd
3210            rundata += 'pf.pdfrange(%d, %6.2f, %6.2f)\n'%(Nd,RMCPdict[dType]['Fitrange'][0],RMCPdict[dType]['Fitrange'][1])
3211            for item in ['dscale','qdamp','qbroad']:
3212                if RMCPdict[dType][item][1]:
3213                    Np += 1
3214                    rundata += 'pf.constrain(pf.%s(),"@%d")\n'%(item,Np)
3215                    parms[Np] = [RMCPdict[dType][item][0],item]
3216    fName = General['Name']+'-PDFfit.stru'
3217    fName = fName.replace(' ','_')
3218    Np = 9
3219    rundata += "pf.read_struct('%s')\n"%(fName)
3220    for item in ['delta1','delta2','sratio']:
3221        if RMCPdict[item][1]:
3222            Np += 1
3223            rundata += 'pf.constrain(pf.%s,"@%d")\n'%(item,Np)
3224            parms[Np] = [RMCPdict[item][0],item]
3225    if 'sphere' in RMCPdict['shape'][0] and RMCPdict['spdiameter'][1]:
3226        Np += 1
3227        rundata += 'pf.constrain(pf.spdiameter,"@%d")\n'%Np
3228        parms[Np] = [RMCPdict['spdiameter'][0],'spdiameter']
3229   
3230    if RMCPdict['cellref']:
3231        cellconst = GetCellConstr(RMCPdict['SGData'])
3232        used = []
3233        cellNames = ['a','b','c','alpha','beta','gamma']
3234        for ic in range(6):
3235            if cellconst[ic]:
3236                rundata += 'pf.constrain(pf.lat(%d), "@%d")\n'%(ic+1,Np+cellconst[ic])
3237                if cellconst[ic] not in used:
3238                    parms[Np+cellconst[ic]] = [Cell[ic],cellNames[ic]]
3239                used.append(cellconst[ic])
3240#Atom constraints here -------------------------------------------------------
3241    AtomVar = RMCPdict['AtomVar']
3242    used = []
3243    for iat,atom in enumerate(RMCPdict['AtomConstr']):
3244        for it,item in enumerate(atom):
3245            names = ['pf.x(%d)'%(iat+1),'pf.y(%d)'%(iat+1),'pf.z(%d)'%(iat+1),'pf.occ(%d)'%(iat+1)]
3246            if it > 1 and item:
3247                itms = item.split('@')
3248                for itm in itms[1:]:
3249                    itnum = int(itm[:2])
3250                    if it < 6:
3251                        rundata += 'pf.constrain(%s,"%s")\n'%(names[it-2],item)
3252                        if itnum not in used:
3253                            parms[itnum] = [AtomVar['@%d'%itnum],names[it-2].split('.')[1]]
3254                            used.append(itnum)
3255                    else:
3256                        uijs = ['pf.u11(%d)'%(iat+1),'pf.u22(%d)'%(iat+1),'pf.u33(%d)'%(iat+1)]     
3257                        for i in range(3):
3258                            rundata += 'pf.constrain(%s,"%s")\n'%(uijs[i],item)
3259                            if itnum not in used:
3260                                parms[itnum] = [AtomVar['@%d'%itnum],uijs[i].split('.')[1]]
3261                                used.append(itnum)
3262                           
3263    if 'sequential' in RMCPdict['refinement']:
3264        rundata += '#parameters here\n'
3265    else:       
3266# set parameter values
3267        RMCPdict['Parms'] = parms
3268        for iprm in parms:
3269            rundata += 'pf.setpar(%d,%.6f)\n'%(iprm,parms[iprm][0])
3270                       
3271# Refine & Save results ---------------------------------------------------------------   
3272    rundata += 'pf.refine()\n'
3273    if 'sequential' in RMCPdict['refinement']:
3274        fName = 'Sequential_PDFfit'
3275    else:
3276        fName = General['Name'].replace(' ','_')+'-PDFfit'
3277    Nd = 0   
3278    for file in RMCPdict['files']:
3279        if 'Select' in RMCPdict['files'][file][0]:
3280            continue
3281        Nd += 1
3282        rundata += 'pf.save_pdf(%d, "%s")\n'%(Nd,fName+file[0]+'.fgr')
3283       
3284    rundata += 'pf.save_struct(1, "%s")\n'%(fName+'.rstr')
3285    rundata += 'pf.save_res("%s")\n'%(fName+'.res')
3286 
3287    if 'sequential' in RMCPdict['refinement']:
3288        rfile = open('Seq_PDFfit_template.py','w')
3289    else:
3290        rfile = open(fName+'.py','w')
3291    rfile.writelines(rundata)
3292    rfile.close()
3293   
3294    return fName+'.py'
3295
3296def GetSeqCell(SGData,parmDict):
3297    ''' For use in processing PDFfit sequential results
3298    '''
3299    try:
3300        if SGData['SGLaue'] in ['m3', 'm3m']:
3301            return [parmDict['11'][0],parmDict['11'][0],parmDict['11'][0],90.,90.,90.]
3302        elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3303            return [parmDict['11'][0],parmDict['11'][0],parmDict['12'][0],90.,90.,90.]
3304        elif SGData['SGLaue'] in ['3R','3mR']:
3305            return [parmDict['11'][0],parmDict['11'][0],parmDict['11'][0],
3306                parmDict['12'][0],parmDict['12'][0],parmDict['12'][0]]
3307        elif SGData['SGLaue'] == 'mmm':
3308            return [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,90.,90.]
3309        elif SGData['SGLaue'] == '2/m':
3310            if SGData['SGUniq'] == 'a':
3311                return [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],parmDict['14'][0],90.,90.]
3312            elif SGData['SGUniq'] == 'b':
3313                return [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,parmDict['14'][0],90.]
3314            elif SGData['SGUniq'] == 'c':
3315                return [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,90.,parmDict['14'][0]]
3316        else:
3317            return [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],
3318                parmDict['14'][0],parmDict['15'][0],parmDict['16'][0]]
3319    except IndexError:
3320        return None
3321   
3322def UpdatePDFfit(Phase,RMCPdict):
3323    ''' Updates various PDFfit parameters held in GSAS-II
3324    '''
3325   
3326    General = Phase['General']
3327    if RMCPdict['refinement'] == 'normal':
3328        fName = General['Name']+'-PDFfit.rstr'
3329        try:
3330            rstr = open(fName.replace(' ','_'),'r')
3331        except FileNotFoundError:
3332            return [fName,'Not found - PDFfit failed']
3333        lines = rstr.readlines()
3334        rstr.close()
3335        header = [line[:-1].split(' ',1) for line in lines[:7]]
3336        resdict = dict(header)
3337        for item in ['scale','sharp','cell']:
3338            resdict[item] = [float(val) for val in resdict[item].split(',')]
3339        General['Cell'][1:7] = resdict['cell']
3340        for inam,name in enumerate(['delta2','delta1','sratio']):
3341            RMCPdict[name][0] = resdict['sharp'][inam]
3342        if 'shape' in resdict:
3343            if 'sphere' in resdict['shape']:
3344                RMCPdict['spdiameter'][0] = resdict['shape'][-1]
3345            else:
3346                RMCPdict['stepcut'][0] = resdict['shape'][-1]
3347        cx,ct,cs,ci = G2mth.getAtomPtrs(Phase)     
3348        Atoms = Phase['Atoms']
3349        atmBeg = 0
3350        for line in lines:
3351            atmBeg += 1
3352            if 'atoms' in line:
3353                break
3354        for atom in Atoms:
3355            atstr = lines[atmBeg][:-1].split()
3356            Uiistr = lines[atmBeg+2][:-1].split()
3357            Uijstr = lines[atmBeg+4][:-1].split()
3358            atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])]
3359            atom[ci] = 'A'
3360            atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])]
3361            atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])]
3362            atmBeg += 6
3363           
3364    fName = 'Sequential_PDFfit.res'
3365    if RMCPdict['refinement'] == 'normal':
3366        fName = General['Name']+'-PDFfit.res'
3367    try:
3368        res = open(fName.replace(' ','_'),'r')
3369    except FileNotFoundError:
3370        return [fName,'Not found - PDFfit failed']
3371    lines = res.readlines()
3372    res.close()
3373    Ibeg = False
3374    resline = ''
3375    for iline,line in enumerate(lines):
3376        if 'Refinement parameters' in line:
3377            Ibeg = True
3378            continue
3379        if Ibeg:
3380            if '---------' in line:
3381                break
3382            resline += line[:-1]
3383    for iline,line in enumerate(lines):
3384        if 'Rw - ' in line:
3385            if 'nan' in line:
3386                Rwp = 100.0
3387            else:
3388                Rwp = float(line.split(':')[1])
3389    results = resline.replace('(','').split(')')[:-1]
3390    results = ['@'+result.lstrip() for result in results]
3391    results = [item.split() for item in results]
3392    if RMCPdict['refinement'] == 'normal':
3393        results = dict([[item[0][:-1],float(item[1])] for item in results if item[0][:-1] in RMCPdict['AtomVar']])
3394        RMCPdict['AtomVar'].update(results)
3395        return None
3396    else:   #sequential
3397        newParms = dict([[item[0][1:-1],[float(item[1]),float(item[2])]] for item in results])
3398        return newParms,Rwp
3399       
3400def MakefullrmcRun(pName,Phase,RMCPdict):
3401    '''Creates a script to run fullrmc. Returns the name of the file that was
3402    created.
3403    '''
3404    BondList = {}
3405    for k in RMCPdict['Pairs']:
3406        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
3407            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
3408    AngleList = []
3409    for angle in RMCPdict['Angles']:
3410        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
3411            continue
3412        for i in (0,1,2):
3413            angle[i] = angle[i].strip()
3414        AngleList.append(angle)
3415    # rmin = RMCPdict['min Contact']
3416    cell = Phase['General']['Cell'][1:7]
3417    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
3418    cx,ct,cs,cia = Phase['General']['AtomPtrs']
3419    atomsList = []
3420    for atom in Phase['Atoms']:
3421        el = ''.join([i for i in atom[ct] if i.isalpha()])
3422        atomsList.append([el] + atom[cx:cx+4])
3423    projDir,projName = os.path.split(os.path.abspath(pName))
3424    scrname = pName+'-fullrmc.py'
3425    restart = '%s_restart.pdb'%pName
3426    Files = RMCPdict['files']
3427    rundata = ''
3428    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
3429    rundata += '# created in '+__file__+" v"+filversion.split()[1]
3430    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
3431    rundata += '''
3432# fullrmc imports (all that are potentially useful)
3433import os,glob
3434import time
3435import pickle
3436import numpy as np
3437from fullrmc.Core import Collection
3438from fullrmc.Engine import Engine
3439import fullrmc.Constraints.PairDistributionConstraints as fPDF
3440from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
3441from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
3442from fullrmc.Constraints.BondConstraints import BondConstraint
3443from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
3444from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
3445from fullrmc.Generators.Swaps import SwapPositionsGenerator
3446# utility routines
3447def writeHeader(ENGINE,statFP):
3448    'header for stats file'
3449    statFP.write('generated-steps, total-error, ')
3450    for c in ENGINE.constraints:
3451        statFP.write(c.constraintName)
3452        statFP.write(', ')
3453    statFP.write('\\n')
3454    statFP.flush()
3455   
3456def writeCurrentStatus(ENGINE,statFP,plotF):
3457    'line in stats file & current constraint plots'
3458    statFP.write(str(ENGINE.generated))
3459    statFP.write(', ')
3460    statFP.write(str(ENGINE.totalStandardError))
3461    statFP.write(', ')
3462    for c in ENGINE.constraints:
3463        statFP.write(str(c.standardError))
3464        statFP.write(', ')
3465    statFP.write('\\n')
3466    statFP.flush()
3467    mpl.use('agg')
3468    fp = open(plotF,'wb')
3469    for c in ENGINE.constraints:
3470        p = c.plot(show=False)
3471        p[0].canvas.draw()
3472        image = p[0].canvas.buffer_rgba()
3473        pickle.dump(c.constraintName,fp)
3474        pickle.dump(np.array(image),fp)
3475    fp.close()
3476
3477def calcRmax(ENGINE):
3478    'from Bachir, works for non-orthorhombic cells'
3479    a,b,c = ENGINE.basisVectors
3480    lens = []
3481    ts    = np.linalg.norm(np.cross(a,b))/2
3482    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3483    ts = np.linalg.norm(np.cross(b,c))/2
3484    lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3485    ts = np.linalg.norm(np.cross(a,c))/2
3486    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3487    return min(lens)
3488'''
3489    rundata += '''
3490### When True, erases an existing engine to provide a fresh start
3491FRESH_START = {:}
3492dirName = "{:}"
3493prefix = "{:}"
3494project = prefix + "-fullrmc"
3495time0 = time.time()
3496'''.format(RMCPdict['ReStart'][0],projDir,projName)
3497   
3498    rundata += '# setup structure\n'
3499    rundata += 'cell = ' + str(cell) + '\n'
3500    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3501    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
3502    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3503
3504    rundata += '\n# initialize engine\n'
3505    rundata += '''
3506engineFileName = os.path.join(dirName, project + '.rmc')
3507projectStats = os.path.join(dirName, project + '.stats')
3508projectPlots = os.path.join(dirName, project + '.plots')
3509pdbFile = os.path.join(dirName, project + '_restart.pdb')
3510# check Engine exists if so (and not FRESH_START) load it
3511# otherwise build it
3512ENGINE = Engine(path=None)
3513if not ENGINE.is_engine(engineFileName) or FRESH_START:
3514    ## create structure
3515    ENGINE = Engine(path=engineFileName, freshStart=True)
3516    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
3517                                 atoms      = atomList,
3518                                 unitcellBC = cell,
3519                                 supercell  = supercell)
3520'''   
3521    import atmdata
3522    # rundata += '    # conversion factors (may be needed)\n'
3523    # rundata += '    sumCiBi2 = 0.\n'
3524    # for elem in Phase['General']['AtomTypes']:
3525    #     rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3526    #     rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3527    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3528    # settings that require a new Engine
3529    for File in Files:
3530        filDat = RMCPdict['files'][File]
3531        if not os.path.exists(filDat[0]): continue
3532        sfwt = 'neutronCohb'
3533        if 'Xray' in File:
3534            sfwt = 'atomicNumber'
3535        if 'G(r)' in File:
3536            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3537            if filDat[3] == 0:
3538                #rundata += '''    # read and xform G(r) as defined in RMCProfile
3539    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3540                #rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
3541                #rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3542                rundata += '    # G(r) as defined in RMCProfile\n'
3543                rundata += '    GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3544            elif filDat[3] == 1:
3545                rundata += '    # This is G(r) as defined in PDFFIT\n'
3546                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3547            elif filDat[3] == 2:
3548                rundata += '    # This is g(r)\n'
3549                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3550            else:
3551                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3552            rundata += '    ENGINE.add_constraints([GofR])\n'
3553            rundata += '    GofR.set_limits((None, calcRmax(ENGINE)))\n'
3554        elif '(Q)' in File:
3555            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3556            if filDat[3] == 0:
3557                rundata += '    # F(Q) as defined in RMCProfile\n'
3558                #rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
3559                if filDat[4]:
3560                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3561                rundata += '    SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3562            elif filDat[3] == 1:
3563                rundata += '    # S(Q) as defined in PDFFIT\n'
3564                rundata += '    SOQ[1] -= 1\n'
3565                if filDat[4]:
3566                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3567                rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3568            else:
3569                raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3570            rundata += '    ENGINE.add_constraints([SofQ])\n'
3571        else:
3572            print('What is this?')
3573    minDists = ''
3574    if BondList:
3575        rundata += '''    B_CONSTRAINT   = BondConstraint()
3576    ENGINE.add_constraints(B_CONSTRAINT)
3577    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3578'''
3579        for pair in BondList:
3580            e1,e2 = pair.split('-')
3581            d1,d2 = BondList[pair]
3582            if d1 == 0: continue
3583            if d2 == 0:
3584                minDists += '("element","{}","{}",{}),'.format(e1.strip(),e2.strip(),d1)
3585            else:
3586                rundata += '            ("element","{}","{}",{},{}),\n'.format(
3587                                        e1.strip(),e2.strip(),d1,d2)
3588        rundata += '             ])\n'
3589    rundata += '    D_CONSTRAINT = DistanceConstraint(defaultLowerDistance={})\n'.format(RMCPdict['min Contact'])
3590    if minDists:
3591        rundata += "    D_CONSTRAINT.set_pairs_definition( {'inter':[" + minDists + "]})\n"
3592    rundata += '    ENGINE.add_constraints(D_CONSTRAINT)\n'
3593   
3594    if AngleList:
3595        rundata += '''    A_CONSTRAINT   = BondsAngleConstraint()
3596    ENGINE.add_constraints(A_CONSTRAINT)
3597    A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3598'''
3599        for item in AngleList:
3600            rundata += ('            '+
3601               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3602        rundata += '             ])\n'
3603    rundata += '''
3604    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3605    ENGINE.save()
3606else:
3607    ENGINE = ENGINE.load(path=engineFileName)
3608
3609ENGINE.set_log_file(os.path.join(dirName,prefix))
3610'''
3611    if RMCPdict['Swaps']:
3612        rundata += '\n#set up for site swaps\n'
3613        rundata += 'aN = ENGINE.allNames\n'
3614        rundata += 'SwapGen = {}\n'
3615        for swap in RMCPdict['Swaps']:
3616            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3617            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3618            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3619        rundata += '    for swaps in SwapGen:\n'
3620        rundata += '        AB = swaps.split("-")\n'
3621        rundata += '        ENGINE.set_groups_as_atoms()\n'
3622        rundata += '        for g in ENGINE.groups:\n'
3623        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
3624        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
3625        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
3626        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
3627        rundata += '            sProb = SwapGen[swaps][2]\n'
3628    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3629    # rundata += 'wtDict = {}\n'
3630    # for File in Files:
3631    #     filDat = RMCPdict['files'][File]
3632    #     if not os.path.exists(filDat[0]): continue
3633    #     if 'Xray' in File:
3634    #         sfwt = 'atomicNumber'
3635    #     else:
3636    #         sfwt = 'neutronCohb'
3637    #     if 'G(r)' in File:
3638    #         typ = 'Pair'
3639    #     elif '(Q)' in File:
3640    #         typ = 'Struct'
3641    #     rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3642    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
3643    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
3644    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3645    rundata += '        c.set_limits((None,calcRmax(ENGINE)))\n'
3646    if RMCPdict['FitScale']:
3647        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3648    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3649    if RMCPdict['FitScale']:
3650        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
3651        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3652    # torsions difficult to implement, must be internal to cell & named with
3653    # fullrmc atom names
3654    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
3655    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
3656    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
3657    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3658    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3659    #     for torsion in RMCPdict['Torsions']:
3660    #         rundata += '    %s\n'%str(tuple(torsion))
3661    #     rundata += '        ]})\n'           
3662    rundata += '''
3663if FRESH_START:
3664    # initialize engine with one step to get starting config energetics
3665    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3666    statFP = open(projectStats,'w')
3667    writeHeader(ENGINE,statFP)
3668    writeCurrentStatus(ENGINE,statFP,projectPlots)
3669else:
3670    statFP = open(projectStats,'a')
3671
3672# setup runs for fullrmc
3673'''
3674    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3675    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3676    rundata += '    ENGINE.set_groups_as_atoms()\n'
3677    rundata += '    expected = ENGINE.generated+steps\n'
3678   
3679    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3680    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3681    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
3682    rundata += 'statFP.close()\n'
3683    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3684    rfile = open(scrname,'w')
3685    rfile.writelines(rundata)
3686    rfile.close()
3687    return scrname
3688   
3689def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3690    bondDist = []
3691    Cell = general['Cell'][1:7]
3692    Supercell =  RMCPdict['SuperCell']
3693    Trans = np.eye(3)*np.array(Supercell)
3694    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3695    Amat,Bmat = G2lat.cell2AB(Cell)
3696    indices = (-1,0,1)
3697    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3698    for bonds in bondList:
3699        Oxyz = np.array(Atoms[bonds[0]][1:])
3700        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
3701        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3702        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3703        for dx in Dx.T:
3704            bondDist.append(np.min(dx))
3705    return np.array(bondDist)
3706   
3707def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3708    bondAngles = []
3709    Cell = general['Cell'][1:7]
3710    Supercell =  RMCPdict['SuperCell']
3711    Trans = np.eye(3)*np.array(Supercell)
3712    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3713    Amat,Bmat = G2lat.cell2AB(Cell)
3714    indices = (-1,0,1)
3715    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3716    for angle in angleList:
3717        Oxyz = np.array(Atoms[angle[0]][1:])
3718        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
3719        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])
3720        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3721        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3722        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3723        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3724        iDAx = np.argmin(DAx,axis=0)
3725        iDBx = np.argmin(DBx,axis=0)
3726        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3727            DAv = DAxV[iA,i]/DAx[iA,i]
3728            DBv = DBxV[iB,i]/DBx[iB,i]
3729            bondAngles.append(npacosd(np.sum(DAv*DBv)))
3730    return np.array(bondAngles)
3731   
3732#### Reflectometry calculations ################################################################################
3733def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3734    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3735   
3736    class RandomDisplacementBounds(object):
3737        """random displacement with bounds"""
3738        def __init__(self, xmin, xmax, stepsize=0.5):
3739            self.xmin = xmin
3740            self.xmax = xmax
3741            self.stepsize = stepsize
3742   
3743        def __call__(self, x):
3744            """take a random step but ensure the new position is within the bounds"""
3745            while True:
3746                # this could be done in a much more clever way, but it will work for example purposes
3747                steps = self.xmax-self.xmin
3748                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
3749                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3750                    break
3751            return xnew
3752   
3753    def GetModelParms():
3754        parmDict = {}
3755        varyList = []
3756        values = []
3757        bounds = []
3758        parmDict['dQ type'] = data['dQ type']
3759        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
3760        for parm in ['Scale','FltBack']:
3761            parmDict[parm] = data[parm][0]
3762            if data[parm][1]:
3763                varyList.append(parm)
3764                values.append(data[parm][0])
3765                bounds.append(Bounds[parm])
3766        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3767        parmDict['nLayers'] = len(parmDict['Layer Seq'])
3768        for ilay,layer in enumerate(data['Layers']):
3769            name = layer['Name']
3770            cid = str(ilay)+';'
3771            parmDict[cid+'Name'] = name
3772            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3773                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3774                if layer.get(parm,[0.,False])[1]:
3775                    varyList.append(cid+parm)
3776                    value = layer[parm][0]
3777                    values.append(value)
3778                    if value:
3779                        bound = [value*Bfac,value/Bfac]
3780                    else:
3781                        bound = [0.,10.]
3782                    bounds.append(bound)
3783            if name not in ['vacuum','unit scatter']:
3784                parmDict[cid+'rho'] = Substances[name]['Scatt density']
3785                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3786        return parmDict,varyList,values,bounds
3787   
3788    def SetModelParms():
3789        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3790        if 'Scale' in varyList:
3791            data['Scale'][0] = parmDict['Scale']
3792            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
3793        G2fil.G2Print (line)
3794        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3795        if 'FltBack' in varyList:
3796            data['FltBack'][0] = parmDict['FltBack']
3797            line += ' esd: %15.3g'%(sigDict['FltBack'])
3798        G2fil.G2Print (line)
3799        for ilay,layer in enumerate(data['Layers']):
3800            name = layer['Name']
3801            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3802            cid = str(ilay)+';'
3803            line = ' '
3804            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3805            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3806            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3807                if parm in layer:
3808                    if parm == 'Rough':
3809                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
3810                    else:
3811                        layer[parm][0] = parmDict[cid+parm]
3812                    line += ' %s: %.3f'%(parm,layer[parm][0])
3813                    if cid+parm in varyList:
3814                        line += ' esd: %.3g'%(sigDict[cid+parm])
3815            G2fil.G2Print (line)
3816            G2fil.G2Print (line2)
3817   
3818    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3819        parmDict.update(zip(varyList,values))
3820        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3821        return M
3822   
3823    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3824        parmDict.update(zip(varyList,values))
3825        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3826        return np.sum(M**2)
3827   
3828    def getREFD(Q,Qsig,parmDict):
3829        Ic = np.ones_like(Q)*parmDict['FltBack']
3830        Scale = parmDict['Scale']
3831        Nlayers = parmDict['nLayers']
3832        Res = parmDict['Res']
3833        depth = np.zeros(Nlayers)
3834        rho = np.zeros(Nlayers)
3835        irho = np.zeros(Nlayers)
3836        sigma = np.zeros(Nlayers)
3837        for ilay,lay in enumerate(parmDict['Layer Seq']):
3838            cid = str(lay)+';'
3839            depth[ilay] = parmDict[cid+'Thick']
3840            sigma[ilay] = parmDict[cid+'Rough']
3841            if parmDict[cid+'Name'] == u'unit scatter':
3842                rho[ilay] = parmDict[cid+'DenMul']
3843                irho[ilay] = parmDict[cid+'iDenMul']
3844            elif 'vacuum' != parmDict[cid+'Name']:
3845                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3846                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3847            if cid+'Mag SLD' in parmDict:
3848                rho[ilay] += parmDict[cid+'Mag SLD']
3849        if parmDict['dQ type'] == 'None':
3850            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3851        elif 'const' in parmDict['dQ type']:
3852            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
3853        else:       #dQ/Q in data
3854            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3855        Ic += AB*Scale
3856        return Ic
3857       
3858    def estimateT0(takestep):
3859        Mmax = -1.e-10
3860        Mmin = 1.e10
3861        for i in range(100):
3862            x0 = takestep(values)
3863            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3864            Mmin = min(M,Mmin)
3865            MMax = max(M,Mmax)
3866        return 1.5*(MMax-Mmin)
3867
3868    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3869    if data.get('2% weight'):
3870        wt = 1./(0.02*Io)**2
3871    Qmin = Limits[1][0]
3872    Qmax = Limits[1][1]
3873    wtFactor = ProfDict['wtFactor']
3874    Bfac = data['Toler']
3875    Ibeg = np.searchsorted(Q,Qmin)
3876    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
3877    Ic[:] = 0
3878    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3879              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3880    parmDict,varyList,values,bounds = GetModelParms()
3881    Msg = 'Failed to converge'
3882    if varyList:
3883        if data['Minimizer'] == 'LMLS': 
3884            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3885                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3886            parmDict.update(zip(varyList,result[0]))
3887            chisq = np.sum(result[2]['fvec']**2)
3888            ncalc = result[2]['nfev']
3889            covM = result[1]
3890            newVals = result[0]
3891        elif data['Minimizer'] == 'Basin Hopping':
3892            xyrng = np.array(bounds).T
3893            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3894            T0 = estimateT0(take_step)
3895            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3896            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3897                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3898                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3899            chisq = result.fun
3900            ncalc = result.nfev
3901            newVals = result.x
3902            covM = []
3903        elif data['Minimizer'] == 'MC/SA Anneal':
3904            xyrng = np.array(bounds).T
3905            result = G2mth.anneal(sumREFD, values, 
3906                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3907                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3908                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3909                ranRange=0.20,autoRan=False,dlg=None)
3910            newVals = result[0]
3911            parmDict.update(zip(varyList,newVals))
3912            chisq = result[1]
3913            ncalc = result[3]
3914            covM = []
3915            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3916        elif data['Minimizer'] == 'L-BFGS-B':
3917            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
3918                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3919            parmDict.update(zip(varyList,result['x']))
3920            chisq = result.fun
3921            ncalc = result.nfev
3922            newVals = result.x
3923            covM = []
3924    else:   #nothing varied
3925        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3926        chisq = np.sum(M**2)
3927        ncalc = 0
3928        covM = []
3929        sig = []
3930        sigDict = {}
3931        result = []
3932    Rvals = {}
3933    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
3934    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
3935    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3936    Ib[Ibeg:Ifin] = parmDict['FltBack']
3937    try:
3938        if not len(varyList):
3939            Msg += ' - nothing refined'
3940            raise ValueError
3941        Nans = np.isnan(newVals)
3942        if np.any(Nans):
3943            name = varyList[Nans.nonzero(True)[0]]
3944            Msg += ' Nan result for '+name+'!'
3945            raise ValueError
3946        Negs = np.less_equal(newVals,0.)
3947        if np.any(Negs):
3948            indx = Negs.nonzero()
3949            name = varyList[indx[0][0]]
3950            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3951                Msg += ' negative coefficient for '+name+'!'
3952                raise ValueError
3953        if len(covM):
3954            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3955            covMatrix = covM*Rvals['GOF']
3956        else:
3957            sig = np.zeros(len(varyList))
3958            covMatrix = []
3959        sigDict = dict(zip(varyList,sig))
3960        G2fil.G2Print (' Results of reflectometry data modelling fit:')
3961        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3962        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3963        SetModelParms()
3964        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3965    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
3966        G2fil.G2Print (Msg)
3967        return False,0,0,0,0,0,0,Msg
3968       
3969def makeSLDprofile(data,Substances):
3970   
3971    sq2 = np.sqrt(2.)
3972    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3973    Nlayers = len(laySeq)
3974    laySeq = np.array(laySeq,dtype=int)
3975    interfaces = np.zeros(Nlayers)
3976    rho = np.zeros(Nlayers)
3977    sigma = np.zeros(Nlayers)
3978    name = data['Layers'][0]['Name']
3979    thick = 0.
3980    for ilay,lay in enumerate(laySeq):
3981        layer = data['Layers'][lay]
3982        name = layer['Name']
3983        if 'Thick' in layer:
3984            thick += layer['Thick'][0]
3985            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3986        if 'Rough' in layer:
3987            sigma[ilay] = max(0.001,layer['Rough'][0])
3988        if name != 'vacuum':
3989            if name == 'unit scatter':
3990                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
3991            else:
3992                rrho = Substances[name]['Scatt density']
3993                irho = Substances[name]['XImag density']
3994                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
3995        if 'Mag SLD' in layer:
3996            rho[ilay] += layer['Mag SLD'][0]
3997    name = data['Layers'][-1]['Name']
3998    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
3999    xr = np.flipud(x)
4000    interfaces[-1] = x[-1]
4001    y = np.ones_like(x)*rho[0]
4002    iBeg = 0
4003    for ilayer in range(Nlayers-1):
4004        delt = rho[ilayer+1]-rho[ilayer]
4005        iPos = np.searchsorted(x,interfaces[ilayer])
4006        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
4007        iBeg = iPos
4008    return x,xr,y   
4009
4010def REFDModelFxn(Profile,Inst,Limits,Substances,data):
4011   
4012    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
4013    Qmin = Limits[1][0]
4014    Qmax = Limits[1][1]
4015    iBeg = np.searchsorted(Q,Qmin)
4016    iFin = np.searchsorted(Q,Qmax)+1    #include last point
4017    Ib[:] = data['FltBack'][0]
4018    Ic[:] = 0
4019    Scale = data['Scale'][0]
4020    if data['Layer Seq'] == []:
4021        return
4022    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
4023    Nlayers = len(laySeq)
4024    depth = np.zeros(Nlayers)
4025    rho = np.zeros(Nlayers)
4026    irho = np.zeros(Nlayers)
4027    sigma = np.zeros(Nlayers)
4028    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
4029        layer = data['Layers'][lay]
4030        name = layer['Name']
4031        if 'Thick' in layer:    #skips first & last layers
4032            depth[ilay] = layer['Thick'][0]
4033        if 'Rough' in layer:    #skips first layer
4034            sigma[ilay] = layer['Rough'][0]
4035        if 'unit scatter' == name:
4036            rho[ilay] = layer['DenMul'][0]
4037            irho[ilay] = layer['iDenMul'][0]
4038        else:
4039            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
4040            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
4041        if 'Mag SLD' in layer:
4042            rho[ilay] += layer['Mag SLD'][0]
4043    if data['dQ type'] == 'None':
4044        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
4045    elif 'const' in data['dQ type']:
4046        res = data['Resolution'][0]/(100.*sateln2)
4047        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
4048    else:       #dQ/Q in data
4049        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
4050    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
4051
4052def abeles(kz, depth, rho, irho=0, sigma=0):
4053    """
4054    Optical matrix form of the reflectivity calculation.
4055    O.S. Heavens, Optical Properties of Thin Solid Films
4056   
4057    Reflectometry as a function of kz for a set of slabs.
4058
4059    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
4060        This is :math:`\\tfrac12 Q_z`.       
4061    :param depth:  float[m] (Ang).
4062        thickness of each layer.  The thickness of the incident medium
4063        and substrate are ignored.
4064    :param rho:  float[n,k] (1e-6/Ang^2)
4065        Real scattering length density for each layer for each kz
4066    :param irho:  float[n,k] (1e-6/Ang^2)
4067        Imaginary scattering length density for each layer for each kz
4068        Note: absorption cross section mu = 2 irho/lambda for neutrons
4069    :param sigma: float[m-1] (Ang)
4070        interfacial roughness.  This is the roughness between a layer
4071        and the previous layer. The sigma array should have m-1 entries.
4072
4073    Slabs are ordered with the surface SLD at index 0 and substrate at
4074    index -1, or reversed if kz < 0.
4075    """
4076    def calc(kz, depth, rho, irho, sigma):
4077        if len(kz) == 0: return kz
4078   
4079        # Complex index of refraction is relative to the incident medium.
4080        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
4081        # in place of kz^2, and ignoring rho_o
4082        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
4083        k = kz
4084   
4085        # According to Heavens, the initial matrix should be [ 1 F; F 1],
4086        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
4087        # multiply versus some coding convenience.
4088        B11 = 1
4089        B22 = 1
4090        B21 = 0
4091        B12 = 0
4092        for i in range(0, len(depth)-1):
4093            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
4094            F = (k - k_next) / (k + k_next)
4095            F *= np.exp(-2*k*k_next*sigma[i]**2)
4096            #print "==== layer",i
4097            #print "kz:", kz
4098            #print "k:", k
4099            #print "k_next:",k_next
4100            #print "F:",F
4101            #print "rho:",rho[:,i+1]
4102            #print "irho:",irho[:,i+1]
4103            #print "d:",depth[i],"sigma:",sigma[i]
4104            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
4105            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
4106            M21 = F*M11
4107            M12 = F*M22
4108            C1 = B11*M11 + B21*M12
4109            C2 = B11*M21 + B21*M22
4110            B11 = C1
4111            B21 = C2
4112            C1 = B12*M11 + B22*M12
4113            C2 = B12*M21 + B22*M22
4114            B12 = C1
4115            B22 = C2
4116            k = k_next
4117   
4118        r = B12/B11
4119        return np.absolute(r)**2
4120
4121    if np.isscalar(kz): kz = np.asarray([kz], 'd')
4122
4123    m = len(depth)
4124
4125    # Make everything into arrays
4126    depth = np.asarray(depth,'d')
4127    rho = np.asarray(rho,'d')
4128    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
4129    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
4130
4131    # Repeat rho,irho columns as needed
4132    if len(rho.shape) == 1:
4133        rho = rho[None,:]
4134        irho = irho[None,:]
4135
4136    return calc(kz, depth, rho, irho, sigma)
4137   
4138def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
4139    y = abeles(kz, depth, rho, irho, sigma)
4140    s = dq/2.
4141    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
4142    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
4143    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
4144    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
4145    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
4146    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
4147    y *= 0.137023
4148    return y
4149       
4150def makeRefdFFT(Limits,Profile):
4151    G2fil.G2Print ('make fft')
4152    Q,Io = Profile[:2]
4153    Qmin = Limits[1][0]
4154    Qmax = Limits[1][1]
4155    iBeg = np.searchsorted(Q,Qmin)
4156    iFin = np.searchsorted(Q,Qmax)+1    #include last point
4157    Qf = np.linspace(0.,Q[iFin-1],5000)
4158    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
4159    If = QI(Qf)*Qf**4
4160    R = np.linspace(0.,5000.,5000)
4161    F = fft.rfft(If)
4162    return R,F
4163
4164   
4165#### Stacking fault simulation codes ################################################################################
4166def GetStackParms(Layers):
4167   
4168    Parms = []
4169#cell parms
4170    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
4171        Parms.append('cellA')
4172        Parms.append('cellC')
4173    else:
4174        Parms.append('cellA')
4175        Parms.append('cellB')
4176        Parms.append('cellC')
4177        if Layers['Laue'] != 'mmm':
4178            Parms.append('cellG')
4179#Transition parms
4180    for iY in range(len(Layers['Layers'])):
4181        for iX in range(len(Layers['Layers'])):
4182            Parms.append('TransP;%d;%d'%(iY,iX))
4183            Parms.append('TransX;%d;%d'%(iY,iX))
4184            Parms.append('TransY;%d;%d'%(iY,iX))
4185            Parms.append('TransZ;%d;%d'%(iY,iX))
4186    return Parms
4187
4188def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
4189    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
4190   
4191    :param dict Layers: dict with following items
4192
4193      ::
4194
4195       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
4196       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
4197       'Layers':[],'Stacking':[],'Transitions':[]}
4198       
4199    :param str ctrls: controls string to be written on DIFFaX controls.dif file
4200    :param float scale: scale factor
4201    :param dict background: background parameters
4202    :param list limits: min/max 2-theta to be calculated
4203    :param dict inst: instrument parameters dictionary
4204    :param list profile: powder pattern data
4205   
4206    Note that parameters all updated in place   
4207    '''
4208    import atmdata
4209    path = sys.path
4210    for name in path:
4211        if 'bin' in name:
4212            DIFFaX = name+'/DIFFaX.exe'
4213            G2fil.G2Print (' Execute '+DIFFaX)
4214            break
4215    # make form factor file that DIFFaX wants - atom types are GSASII style
4216    sf = open('data.sfc','w')
4217    sf.write('GSASII special form factor file for DIFFaX\n\n')
4218    atTypes = list(Layers['AtInfo'].keys())
4219    if 'H' not in atTypes:
4220        atTypes.insert(0,'H')
4221    for atType in atTypes:
4222        if atType == 'H': 
4223            blen = -.3741
4224        else:
4225            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
4226        Adat = atmdata.XrayFF[atType]
4227        text = '%4s'%(atType.ljust(4))
4228        for i in range(4):
4229            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
4230        text += '%11.6f%11.6f'%(Adat['fc'],blen)
4231        text += '%3d\n'%(Adat['Z'])
4232        sf.write(text)
4233    sf.close()
4234    #make DIFFaX control.dif file - future use GUI to set some of these flags
4235    cf = open('control.dif','w')
4236    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4237        x0 = profile[0]
4238        iBeg = np.searchsorted(x0,limits[0])
4239        iFin = np.searchsorted(x0,limits[1])+1
4240        if iFin-iBeg > 20000:
4241            iFin = iBeg+20000
4242        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
4243        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4244        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
4245    else:
4246        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4247        inst = {'Type':['XSC','XSC',]}
4248    cf.close()
4249    #make DIFFaX data file
4250    df = open('GSASII-DIFFaX.dat','w')
4251    df.write('INSTRUMENTAL\n')
4252    if 'X' in inst['Type'][0]:
4253        df.write('X-RAY\n')
4254    elif 'N' in inst['Type'][0]:
4255        df.write('NEUTRON\n')
4256    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4257        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
4258        U = ateln2*inst['U'][1]/10000.
4259        V = ateln2*inst['V'][1]/10000.
4260        W = ateln2*inst['W'][1]/10000.
4261        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
4262        HW = np.sqrt(np.mean(HWHM))
4263    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
4264        if 'Mean' in Layers['selInst']:
4265            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
4266        elif 'Gaussian' in Layers['selInst']:
4267            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
4268        else:
4269            df.write('None\n')
4270    else:
4271        df.write('0.10\nNone\n')
4272    df.write('STRUCTURAL\n')
4273    a,b,c = Layers['Cell'][1:4]
4274    gam = Layers['Cell'][6]
4275    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
4276    laue = Layers['Laue']
4277    if laue == '2/m(ab)':
4278        laue = '2/m(1)'
4279    elif laue == '2/m(c)':
4280        laue = '2/m(2)'
4281    if 'unknown' in Layers['Laue']:
4282        df.write('%s %.3f\n'%(laue,Layers['Toler']))
4283    else:   
4284        df.write('%s\n'%(laue))
4285    df.write('%d\n'%(len(Layers['Layers'])))
4286    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
4287        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
4288    layerNames = []
4289    for layer in Layers['Layers']:
4290        layerNames.append(layer['Name'])
4291    for il,layer in enumerate(Layers['Layers']):
4292        if layer['SameAs']:
4293            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
4294            continue
4295        df.write('LAYER %d\n'%(il+1))
4296        if '-1' in layer['Symm']:
4297            df.write('CENTROSYMMETRIC\n')
4298        else:
4299            df.write('NONE\n')
4300        for ia,atom in enumerate(layer['Atoms']):
4301            [name,atype,x,y,z,frac,Uiso] = atom
4302            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
4303                frac /= 2.
4304            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
4305    df.write('STACKING\n')
4306    df.write('%s\n'%(Layers['Stacking'][0]))
4307    if 'recursive' in Layers['Stacking'][0]:
4308        df.write('%s\n'%Layers['Stacking'][1])
4309    else:
4310        if 'list' in Layers['Stacking'][1]:
4311            Slen = len(Layers['Stacking'][2])
4312            iB = 0
4313            iF = 0
4314            while True:
4315                iF += 68
4316                if iF >= Slen:
4317                    break
4318                iF = min(iF,Slen)
4319                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))<