source: trunk/GSASIIpwd.py @ 5144

Last change on this file since 5144 was 5144, checked in by vondreele, 7 months ago

G2pwd - close file left unclosed (fqdata)
skip generating SCATTERING-LENGTH & NEUTRON-COEFFICIENT records for x-ray data for RMCProfile
G2phsGUI - revise RMC FileSizer? routine; shift FitScale? box back to original position (before files)
remove Proc.wait for RMCProfile & fullrmc runs - no point & prevented examination of plots. Still used for PDFfit runs.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 196.4 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-17 21:37:39 +0000 (Mon, 17 Jan 2022) $
10# $Author: vondreele $
11# $Revision: 5144 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 5144 2022-01-17 21:37:39Z vondreele $
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: 5144 $"
38GSASIIpath.SetVersionNumber("$Revision: 5144 $")
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    else:
903        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
904        q = npT2q(xdata,wave)
905    yb = np.zeros_like(xdata)
906    nBak = 0
907    sumBk = [0.,0.,0]
908    while True:
909        key = pfx+'Back;'+str(nBak)
910        if key in parmDict:
911            nBak += 1
912        else:
913            break
914#empirical functions
915    if bakType in ['chebyschev','cosine','chebyschev-1']:
916        dt = xdata[-1]-xdata[0]   
917        for iBak in range(nBak):
918            key = pfx+'Back;'+str(iBak)
919            if bakType == 'chebyschev':
920                ybi = parmDict[key]*(-1.+2.*(xdata-xdata[0])/dt)**iBak
921            elif bakType == 'chebyschev-1':
922                xpos = -1.+2.*(xdata-xdata[0])/dt
923                ybi = parmDict[key]*np.cos(iBak*np.arccos(xpos))
924            elif bakType == 'cosine':
925                ybi = parmDict[key]*npcosd(180.*xdata*iBak/xdata[-1])
926            yb += ybi
927        sumBk[0] = np.sum(yb)
928    elif bakType in ['Q^2 power series','Q^-2 power series']:
929        QT = 1.
930        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
931        for iBak in range(nBak-1):
932            key = pfx+'Back;'+str(iBak+1)
933            if '-2' in bakType:
934                QT *= (iBak+1)*q**-2
935            else:
936                QT *= q**2/(iBak+1)
937            yb += QT*parmDict[key]
938        sumBk[0] = np.sum(yb)
939    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
940        if nBak == 1:
941            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
942        elif nBak == 2:
943            dX = xdata[-1]-xdata[0]
944            T2 = (xdata-xdata[0])/dX
945            T1 = 1.0-T2
946            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
947        else:
948            xnomask = ma.getdata(xdata)
949            xmin,xmax = xnomask[0],xnomask[-1]
950            if bakType == 'lin interpolate':
951                bakPos = np.linspace(xmin,xmax,nBak,True)
952            elif bakType == 'inv interpolate':
953                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
954            elif bakType == 'log interpolate':
955                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
956            bakPos[0] = xmin
957            bakPos[-1] = xmax
958            bakVals = np.zeros(nBak)
959            for i in range(nBak):
960                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
961            bakInt = si.interp1d(bakPos,bakVals,'linear')
962            yb = bakInt(ma.getdata(xdata))
963        sumBk[0] = np.sum(yb)
964#Debye function
965    if pfx+'difC' in parmDict or 'E' in parmDict[pfx+'Type']:
966        ff = 1.
967    else:       
968        try:
969            wave = parmDict[pfx+'Lam']
970        except KeyError:
971            wave = parmDict[pfx+'Lam1']
972        SQ = (q/(4.*np.pi))**2
973        FF = G2elem.GetFormFactorCoeff('Si')[0]
974        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
975    iD = 0       
976    while True:
977        try:
978            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
979            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
980            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
981            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
982            yb += ybi
983            sumBk[1] += np.sum(ybi)
984            iD += 1       
985        except KeyError:
986            break
987#peaks
988    iD = 0
989    while True:
990        try:
991            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
992            pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1)
993            pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],1.)
994            pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1)
995            if 'C' in dataType:
996                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
997            else: #'T'OF
998                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
999            iBeg = np.searchsorted(xdata,pkP-fmin)
1000            iFin = np.searchsorted(xdata,pkP+fmax)
1001            lenX = len(xdata)
1002            if not iBeg:
1003                iFin = np.searchsorted(xdata,pkP+fmax)
1004            elif iBeg == lenX:
1005                iFin = iBeg
1006            else:
1007                iFin = np.searchsorted(xdata,pkP+fmax)
1008            if 'C' in dataType:
1009                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])[0]
1010                yb[iBeg:iFin] += ybi
1011            elif 'T' in dataType:
1012                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])[0]
1013                yb[iBeg:iFin] += ybi
1014            elif 'B' in dataType:
1015                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin])[0]
1016                yb[iBeg:iFin] += ybi
1017            sumBk[2] += np.sum(ybi)
1018            iD += 1       
1019        except KeyError:
1020            break
1021        except ValueError:
1022            G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
1023            break
1024    if fixback is not None:   
1025        yb += parmDict[pfx+'BF mult']*fixback
1026        sumBk[0] = sum(yb)
1027    return yb,sumBk
1028   
1029def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata,fixback=None):
1030    'needs a doc string'
1031    if 'T' in dataType:
1032        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
1033    else:
1034        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1035        q = 2.*np.pi*npsind(xdata/2.)/wave
1036    nBak = 0
1037    while True:
1038        key = hfx+'Back;'+str(nBak)
1039        if key in parmDict:
1040            nBak += 1
1041        else:
1042            break
1043    dydb = np.zeros(shape=(nBak,len(xdata)))
1044    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
1045    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
1046    dydfb = []
1047
1048    if bakType in ['chebyschev','cosine','chebyschev-1']:
1049        dt = xdata[-1]-xdata[0]   
1050        for iBak in range(nBak):   
1051            if bakType == 'chebyschev':
1052                dydb[iBak] = (-1.+2.*(xdata-xdata[0])/dt)**iBak
1053            elif bakType == 'chebyschev-1':
1054                xpos = -1.+2.*(xdata-xdata[0])/dt
1055                dydb[iBak] = np.cos(iBak*np.arccos(xpos))
1056            elif bakType == 'cosine':
1057                dydb[iBak] = npcosd(180.*xdata*iBak/xdata[-1])
1058    elif bakType in ['Q^2 power series','Q^-2 power series']:
1059        QT = 1.
1060        dydb[0] = np.ones_like(xdata)
1061        for iBak in range(nBak-1):
1062            if '-2' in bakType:
1063                QT *= (iBak+1)*q**-2
1064            else:
1065                QT *= q**2/(iBak+1)
1066            dydb[iBak+1] = QT
1067    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
1068        if nBak == 1:
1069            dydb[0] = np.ones_like(xdata)
1070        elif nBak == 2:
1071            dX = xdata[-1]-xdata[0]
1072            T2 = (xdata-xdata[0])/dX
1073            T1 = 1.0-T2
1074            dydb = [T1,T2]
1075        else:
1076            xnomask = ma.getdata(xdata)
1077            xmin,xmax = xnomask[0],xnomask[-1]
1078            if bakType == 'lin interpolate':
1079                bakPos = np.linspace(xmin,xmax,nBak,True)
1080            elif bakType == 'inv interpolate':
1081                bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True)
1082            elif bakType == 'log interpolate':
1083                bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True))
1084            bakPos[0] = xmin
1085            bakPos[-1] = xmax
1086            for i,pos in enumerate(bakPos):
1087                if i == 0:
1088                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
1089                elif i == len(bakPos)-1:
1090                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
1091                else:
1092                    dydb[i] = np.where(xdata>bakPos[i],
1093                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
1094                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
1095    if hfx+'difC' in parmDict:
1096        ff = 1.
1097    else:
1098        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
1099        q = npT2q(xdata,wave)
1100        SQ = (q/(4*np.pi))**2
1101        FF = G2elem.GetFormFactorCoeff('Si')[0]
1102        ff = np.array(G2elem.ScatFac(FF,SQ)[0])*np.pi**2    #needs pi^2~10. for cw data (why?)
1103    iD = 0       
1104    while True:
1105        try:
1106            if hfx+'difC' in parmDict:
1107                q = 2*np.pi*parmDict[hfx+'difC']/xdata
1108            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
1109            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
1110            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
1111            sqr = np.sin(q*dbR)/(q*dbR)
1112            cqr = np.cos(q*dbR)
1113            temp = np.exp(-dbU*q**2)
1114            dyddb[3*iD] = ff*sqr*temp
1115            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(dbR)
1116            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2
1117            iD += 1
1118        except KeyError:
1119            break
1120    iD = 0
1121    while True:
1122        try:
1123            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
1124            pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1)
1125            pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],1.0)
1126            pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1)
1127            if 'C' in dataType:
1128                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
1129            else: #'T' or 'B'
1130                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
1131            iBeg = np.searchsorted(xdata,pkP-fmin)
1132            iFin = np.searchsorted(xdata,pkP+fmax)
1133            lenX = len(xdata)
1134            if not iBeg:
1135                iFin = np.searchsorted(xdata,pkP+fmax)
1136            elif iBeg == lenX:
1137                iFin = iBeg
1138            else:
1139                iFin = np.searchsorted(xdata,pkP+fmax)
1140            if 'C' in dataType:
1141                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
1142            else:   #'T'OF
1143                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
1144            dydpk[4*iD][iBeg:iFin] += pkI*dFdp
1145            dydpk[4*iD+1][iBeg:iFin] += Df
1146            dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
1147            dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
1148            iD += 1       
1149        except KeyError:
1150            break
1151        except ValueError:
1152            G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****')
1153            break       
1154    # fixed background from file
1155    if  fixback is not None:
1156        dydfb = fixback
1157    return dydb,dyddb,dydpk,dydfb
1158
1159#### Using old gsas fortran routines for powder peak shapes & derivatives
1160def getFCJVoigt3(pos,sig,gam,shl,xdata):
1161    '''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
1162    CW powder peak in external Fortran routine
1163   
1164    param pos: peak position in degrees
1165    param sig: Gaussian variance in centideg^2
1166    param gam: Lorentzian width in centideg
1167    param shl: axial divergence parameter (S+H)/L
1168    param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1169   
1170    returns: array: calculated peak function at each xdata
1171    returns: integral of peak; nominally = 1.0
1172    '''
1173    if len(xdata):
1174        cw = np.diff(xdata)
1175        cw = np.append(cw,cw[-1])
1176        Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1177        return Df,np.sum(100.*Df*cw)
1178    else:
1179        return 0.,1.
1180
1181def getdFCJVoigt3(pos,sig,gam,shl,xdata):
1182    '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
1183    function for a CW powder peak
1184   
1185    param pos: peak position in degrees
1186    param sig: Gaussian variance in centideg^2
1187    param gam: Lorentzian width in centideg
1188    param shl: axial divergence parameter (S+H)/L
1189    param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv)
1190   
1191    returns: arrays: function and derivatives of pos, sig, gam, & shl
1192    '''
1193    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
1194    return Df,dFdp,dFds,dFdg,dFdsh
1195
1196def getPsVoigt(pos,sig,gam,xdata):
1197    '''Compute the simple Pseudo-Voigt function for a
1198    small angle Bragg peak in external Fortran routine
1199   
1200    param pos: peak position in degrees
1201    param sig: Gaussian variance in centideg^2
1202    param gam: Lorentzian width in centideg
1203   
1204    returns: array: calculated peak function at each xdata
1205    returns: integral of peak; nominally = 1.0
1206    '''
1207   
1208    cw = np.diff(xdata)
1209    cw = np.append(cw,cw[-1])
1210    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
1211    return Df,np.sum(100.*Df*cw)
1212
1213def getdPsVoigt(pos,sig,gam,xdata):
1214    '''Compute the simple Pseudo-Voigt function derivatives for a
1215    small angle Bragg peak peak in external Fortran routine
1216   
1217    param pos: peak position in degrees
1218    param sig: Gaussian variance in centideg^2
1219    param gam: Lorentzian width in centideg
1220
1221    returns: arrays: function and derivatives of pos, sig & gam
1222    NB: the pos derivative has the opposite sign compared to that in other profile functions
1223    '''
1224   
1225    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
1226    return Df,dFdp,dFds,dFdg
1227
1228def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
1229    '''Compute the double exponential Pseudo-Voigt convolution function for a
1230    neutron TOF & CW pink peak in external Fortran routine
1231    '''
1232   
1233    cw = np.diff(xdata)
1234    cw = np.append(cw,cw[-1])
1235    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1236    return Df,np.sum(Df*cw) 
1237   
1238def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
1239    '''Compute the double exponential Pseudo-Voigt convolution function derivatives for a
1240    neutron TOF & CW pink peak in external Fortran routine
1241    '''
1242   
1243    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
1244    return Df,dFdp,dFda,dFdb,dFds,dFdg   
1245
1246def ellipseSize(H,Sij,GB):
1247    '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady
1248    '''
1249   
1250    HX = np.inner(H.T,GB)
1251    lenHX = np.sqrt(np.sum(HX**2))
1252    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
1253    R = np.inner(HX/lenHX,Rsize)**2*Esize         #want column length for hkl in crystal
1254    lenR = 1./np.sqrt(np.sum(R))
1255    return lenR
1256
1257def ellipseSizeDerv(H,Sij,GB):
1258    '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) derivative per note from Alexander Brady
1259    '''
1260   
1261    lenR = ellipseSize(H,Sij,GB)
1262    delt = 0.001
1263    dRdS = np.zeros(6)
1264    for i in range(6):
1265        Sij[i] -= delt
1266        lenM = ellipseSize(H,Sij,GB)
1267        Sij[i] += 2.*delt
1268        lenP = ellipseSize(H,Sij,GB)
1269        Sij[i] -= delt
1270        dRdS[i] = (lenP-lenM)/(2.*delt)
1271    return lenR,dRdS
1272
1273def getMustrain(HKL,G,SGData,muStrData):
1274    if muStrData[0] == 'isotropic':
1275        return np.ones(HKL.shape[1])*muStrData[1][0]
1276    elif muStrData[0] == 'uniaxial':
1277        H = np.array(HKL)
1278        P = np.array(muStrData[3])
1279        cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1280        Si = muStrData[1][0]
1281        Sa = muStrData[1][1]
1282        return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1283    else:       #generalized - P.W. Stephens model
1284        H = np.array(HKL)
1285        rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T])
1286        Strms = np.array(G2spc.MustrainCoeff(H,SGData))
1287        Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0)
1288        return np.sqrt(Sum)/rdsq
1289   
1290def getCrSize(HKL,G,GB,sizeData):
1291    if sizeData[0] == 'isotropic':
1292        return np.ones(HKL.shape[1])*sizeData[1][0]
1293    elif sizeData[0] == 'uniaxial':
1294        H = np.array(HKL)
1295        P = np.array(sizeData[3])
1296        cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T
1297        Si = sizeData[1][0]
1298        Sa = sizeData[1][1]
1299        return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1300    else:
1301        Sij =[sizeData[4][i] for i in range(6)]
1302        H = np.array(HKL)
1303        return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2
1304
1305def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False):
1306    '''
1307    Generates allowed by symmetry reflections with d >= dmin
1308    NB: GenHKLf & checkMagextc return True for extinct reflections
1309
1310    :param dmin:  minimum d-spacing
1311    :param SGData: space group data obtained from SpcGroup
1312    :param A: lattice parameter terms A1-A6
1313    :param Inst: instrument parameter info
1314    :returns: HKLs: np.array hkl, etc for allowed reflections
1315
1316    '''
1317    HKL = G2lat.GenHLaue(dmin,SGData,A)       
1318    HKLs = []
1319    ds = []
1320    for h,k,l,d in HKL:
1321        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1322        if ext and 'MagSpGrp' in SGData:
1323            ext = G2spc.checkMagextc([h,k,l],SGData)
1324        if not ext:
1325            if nodup and int(10000*d) in ds:
1326                continue
1327            ds.append(int(10000*d))
1328            if Inst == None:
1329                HKLs.append([h,k,l,d,0,-1])
1330            else:
1331                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
1332    return np.array(HKLs)
1333
1334def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
1335    'needs a doc string'
1336    HKLs = []
1337    vec = np.array(Vec)
1338    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1339    dvec = 1./(maxH*vstar+1./dmin)
1340    HKL = G2lat.GenHLaue(dvec,SGData,A)       
1341    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1342    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1343    ifMag = False
1344    if 'MagSpGrp' in SGData:
1345        ifMag = True
1346    for h,k,l,d in HKL:
1347        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
1348        if not ext and d >= dmin:
1349            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
1350        for dH in SSdH:
1351            if dH:
1352                DH = SSdH[dH]
1353                H = [h+DH[0],k+DH[1],l+DH[2]]
1354                d = float(1/np.sqrt(G2lat.calc_rDsq(H,A)))
1355                if d >= dmin:
1356                    HKLM = np.array([h,k,l,dH])
1357                    if G2spc.checkSSextc(HKLM,SSGData) or ifMag:
1358                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
1359    return G2lat.sortHKLd(HKLs,True,True,True)
1360
1361peakInstPrmMode = True
1362'''Determines the mode used for peak fitting. When peakInstPrmMode=True peak
1363width parameters are computed from the instrument parameters (UVW,... or
1364alpha,... etc) unless the individual parameter is refined. This allows the
1365instrument parameters to be refined. When peakInstPrmMode=False, the instrument
1366parameters are not used and cannot be refined.
1367The default is peakFitMode=True. This is changed only in
1368:func:`setPeakInstPrmMode`, which is called only from :mod:`GSASIIscriptable`.
1369'''
1370
1371def setPeakInstPrmMode(normal=True):
1372    '''Determines the mode used for peak fitting. If normal=True (default)
1373    peak width parameters are computed from the instrument parameters
1374    unless the individual parameter is refined. If normal=False,
1375    peak widths are used as supplied for each peak.
1376    At present this is called only in G2scriptable.
1377
1378    Note that normal=True unless this routine is called. Also,
1379    instrument parameters can only be refined with normal=True.
1380
1381    :param bool normal: setting to apply to global variable
1382      :data:`peakInstPrmMode`
1383    '''
1384    global peakInstPrmMode
1385    peakInstPrmMode = normal
1386
1387def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType):
1388    '''Computes the profiles from multiple peaks for individual peak fitting
1389    for powder patterns.
1390    NB: not used for Rietveld refinement
1391    '''
1392   
1393    yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
1394    yc = np.zeros_like(yb)
1395    # cw = np.diff(xdata)
1396    # cw = np.append(cw,cw[-1])
1397    if 'C' in dataType:
1398        shl = max(parmDict['SH/L'],0.002)
1399        Ka2 = False
1400        if 'Lam1' in parmDict.keys():
1401            Ka2 = True
1402            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1403            kRatio = parmDict['I(L2)/I(L1)']
1404        iPeak = 0
1405        while True:
1406            try:
1407                pos = parmDict['pos'+str(iPeak)]
1408                tth = (pos-parmDict['Zero'])
1409                intens = parmDict['int'+str(iPeak)]
1410                sigName = 'sig'+str(iPeak)
1411                if sigName in varyList or not peakInstPrmMode:
1412                    sig = parmDict[sigName]
1413                else:
1414                    sig = G2mth.getCWsig(parmDict,tth)
1415                sig = max(sig,0.001)          #avoid neg sigma^2
1416                gamName = 'gam'+str(iPeak)
1417                if gamName in varyList or not peakInstPrmMode:
1418                    gam = parmDict[gamName]
1419                else:
1420                    gam = G2mth.getCWgam(parmDict,tth)
1421                gam = max(gam,0.001)             #avoid neg gamma
1422                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1423                iBeg = np.searchsorted(xdata,pos-fmin)
1424                iFin = np.searchsorted(xdata,pos+fmin)
1425                if not iBeg+iFin:       #peak below low limit
1426                    iPeak += 1
1427                    continue
1428                elif not iBeg-iFin:     #peak above high limit
1429                    return yb+yc
1430                fp = getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])[0]
1431                yc[iBeg:iFin] += intens*fp
1432                if Ka2:
1433                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1434                    iBeg = np.searchsorted(xdata,pos2-fmin)
1435                    iFin = np.searchsorted(xdata,pos2+fmin)
1436                    if iBeg-iFin:
1437                        fp2 = getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])[0]
1438                        yc[iBeg:iFin] += intens*kRatio*fp2
1439                iPeak += 1
1440            except KeyError:        #no more peaks to process
1441                return yb+yc
1442    elif 'E' in dataType:
1443        iPeak = 0
1444        dsp = 1.0 #for now - fix later
1445        while True:
1446            try:
1447                pos = parmDict['pos'+str(iPeak)]
1448                intens = parmDict['int'+str(iPeak)]
1449                sigName = 'sig'+str(iPeak)
1450                if sigName in varyList or not peakInstPrmMode:
1451                    sig = parmDict[sigName]
1452                else:
1453                    sig = G2mth.getEDsig(parmDict,pos)
1454                sig = max(sig,0.001)          #avoid neg sigma^2
1455                Wd,fmin,fmax = getWidthsED(pos,sig)
1456                iBeg = np.searchsorted(xdata,pos-fmin)
1457                iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin))
1458                if not iBeg+iFin:       #peak below low limit
1459                    iPeak += 1
1460                    continue
1461                elif not iBeg-iFin:     #peak above high limit
1462                    return yb+yc
1463                yc[iBeg:iFin] += intens*getPsVoigt(pos,sig*10.**4,0.001,xdata[iBeg:iFin])[0]
1464                iPeak += 1
1465            except KeyError:        #no more peaks to process
1466                return yb+yc       
1467    elif 'B' in dataType:
1468        iPeak = 0
1469        dsp = 1.0 #for now - fix later
1470        while True:
1471            try:
1472                pos = parmDict['pos'+str(iPeak)]
1473                tth = (pos-parmDict['Zero'])
1474                intens = parmDict['int'+str(iPeak)]
1475                alpName = 'alp'+str(iPeak)
1476                if alpName in varyList or not peakInstPrmMode:
1477                    alp = parmDict[alpName]
1478                else:
1479                    alp = G2mth.getPinkalpha(parmDict,tth)
1480                alp = max(0.1,alp)
1481                betName = 'bet'+str(iPeak)
1482                if betName in varyList or not peakInstPrmMode:
1483                    bet = parmDict[betName]
1484                else:
1485                    bet = G2mth.getPinkbeta(parmDict,tth)
1486                bet = max(0.1,bet)
1487                sigName = 'sig'+str(iPeak)
1488                if sigName in varyList or not peakInstPrmMode:
1489                    sig = parmDict[sigName]
1490                else:
1491                    sig = G2mth.getCWsig(parmDict,tth)
1492                sig = max(sig,0.001)          #avoid neg sigma^2
1493                gamName = 'gam'+str(iPeak)
1494                if gamName in varyList or not peakInstPrmMode:
1495                    gam = parmDict[gamName]
1496                else:
1497                    gam = G2mth.getCWgam(parmDict,tth)
1498                gam = max(gam,0.001)             #avoid neg gamma
1499                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1500                iBeg = np.searchsorted(xdata,pos-fmin)
1501                iFin = np.searchsorted(xdata,pos+fmin)
1502                if not iBeg+iFin:       #peak below low limit
1503                    iPeak += 1
1504                    continue
1505                elif not iBeg-iFin:     #peak above high limit
1506                    return yb+yc
1507                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])[0]
1508                iPeak += 1
1509            except KeyError:        #no more peaks to process
1510                return yb+yc       
1511    else:
1512        Pdabc = parmDict['Pdabc']
1513        difC = parmDict['difC']
1514        iPeak = 0
1515        while True:
1516            try:
1517                pos = parmDict['pos'+str(iPeak)]               
1518                tof = pos-parmDict['Zero']
1519                dsp = tof/difC
1520                intens = parmDict['int'+str(iPeak)]
1521                alpName = 'alp'+str(iPeak)
1522                if alpName in varyList or not peakInstPrmMode:
1523                    alp = parmDict[alpName]
1524                else:
1525                    if len(Pdabc):
1526                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1527                    else:
1528                        alp = G2mth.getTOFalpha(parmDict,dsp)
1529                alp = max(0.1,alp)
1530                betName = 'bet'+str(iPeak)
1531                if betName in varyList or not peakInstPrmMode:
1532                    bet = parmDict[betName]
1533                else:
1534                    if len(Pdabc):
1535                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1536                    else:
1537                        bet = G2mth.getTOFbeta(parmDict,dsp)
1538                bet = max(0.0001,bet)
1539                sigName = 'sig'+str(iPeak)
1540                if sigName in varyList or not peakInstPrmMode:
1541                    sig = parmDict[sigName]
1542                else:
1543                    sig = G2mth.getTOFsig(parmDict,dsp)
1544                gamName = 'gam'+str(iPeak)
1545                if gamName in varyList or not peakInstPrmMode:
1546                    gam = parmDict[gamName]
1547                else:
1548                    gam = G2mth.getTOFgamma(parmDict,dsp)
1549                gam = max(gam,0.001)             #avoid neg gamma
1550                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1551                iBeg = np.searchsorted(xdata,pos-fmin)
1552                iFin = np.searchsorted(xdata,pos+fmax)
1553                lenX = len(xdata)
1554                if not iBeg:
1555                    iFin = np.searchsorted(xdata,pos+fmax)
1556                elif iBeg == lenX:
1557                    iFin = iBeg
1558                else:
1559                    iFin = np.searchsorted(xdata,pos+fmax)
1560                if not iBeg+iFin:       #peak below low limit
1561                    iPeak += 1
1562                    continue
1563                elif not iBeg-iFin:     #peak above high limit
1564                    return yb+yc
1565                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])[0]
1566                iPeak += 1
1567            except KeyError:        #no more peaks to process
1568                return yb+yc
1569           
1570def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType):
1571    '''Computes the profile derivatives for a powder pattern for single peak fitting
1572   
1573    return: np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
1574   
1575    NB: not used for Rietveld refinement
1576    '''
1577    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
1578    dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback)
1579    if 'Back;0' in varyList:            #background derivs are in front if present
1580        dMdv[0:len(dMdb)] = dMdb
1581    names = ['DebyeA','DebyeR','DebyeU']
1582    for name in varyList:
1583        if 'Debye' in name:
1584            parm,Id = name.split(';')
1585            ip = names.index(parm)
1586            dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip]
1587    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1588    for name in varyList:
1589        if 'BkPk' in name:
1590            parm,Id = name.split(';')
1591            ip = names.index(parm)
1592            dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip]
1593    cw = np.diff(xdata)
1594    cw = np.append(cw,cw[-1])
1595    if 'C' in dataType:
1596        shl = max(parmDict['SH/L'],0.002)
1597        Ka2 = False
1598        if 'Lam1' in parmDict.keys():
1599            Ka2 = True
1600            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1601            kRatio = parmDict['I(L2)/I(L1)']
1602        iPeak = 0
1603        while True:
1604            try:
1605                pos = parmDict['pos'+str(iPeak)]
1606                tth = (pos-parmDict['Zero'])
1607                intens = parmDict['int'+str(iPeak)]
1608                sigName = 'sig'+str(iPeak)
1609                if sigName in varyList or not peakInstPrmMode:
1610                    sig = parmDict[sigName]
1611                    dsdU = dsdV = dsdW = 0
1612                else:
1613                    sig = G2mth.getCWsig(parmDict,tth)
1614                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1615                sig = max(sig,0.001)          #avoid neg sigma
1616                gamName = 'gam'+str(iPeak)
1617                if gamName in varyList or not peakInstPrmMode:
1618                    gam = parmDict[gamName]
1619                    dgdX = dgdY = dgdZ = 0
1620                else:
1621                    gam = G2mth.getCWgam(parmDict,tth)
1622                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1623                gam = max(gam,0.001)             #avoid neg gamma
1624                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1625                iBeg = np.searchsorted(xdata,pos-fmin)
1626                iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin))
1627                if not iBeg+iFin:       #peak below low limit
1628                    iPeak += 1
1629                    continue
1630                elif not iBeg-iFin:     #peak above high limit
1631                    break
1632                dMdpk = np.zeros(shape=(6,len(xdata)))
1633                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1634                for i in range(1,5):
1635                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1636                dMdpk[0][iBeg:iFin] += dMdipk[0]
1637                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1638                if Ka2:
1639                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1640                    iBeg = np.searchsorted(xdata,pos2-fmin)
1641                    iFin = np.searchsorted(xdata,pos2+fmin)
1642                    if iBeg-iFin:
1643                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1644                        for i in range(1,5):
1645                            dMdpk[i][iBeg:iFin] += intens*kRatio*dMdipk2[i]
1646                        dMdpk[0][iBeg:iFin] += kRatio*dMdipk2[0]
1647                        dMdpk[5][iBeg:iFin] += dMdipk2[0]
1648                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1649                for parmName in ['pos','int','sig','gam']:
1650                    try:
1651                        idx = varyList.index(parmName+str(iPeak))
1652                        dMdv[idx] = dervDict[parmName]
1653                    except ValueError:
1654                        pass
1655                if 'U' in varyList:
1656                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1657                if 'V' in varyList:
1658                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1659                if 'W' in varyList:
1660                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1661                if 'X' in varyList:
1662                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1663                if 'Y' in varyList:
1664                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1665                if 'Z' in varyList:
1666                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1667                if 'SH/L' in varyList:
1668                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1669                if 'I(L2)/I(L1)' in varyList:
1670                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1671                iPeak += 1
1672            except KeyError:        #no more peaks to process
1673                break
1674    elif 'E' in dataType:
1675        iPeak = 0
1676        while True:
1677            try:
1678                pos = parmDict['pos'+str(iPeak)]
1679                intens = parmDict['int'+str(iPeak)]
1680                sigName = 'sig'+str(iPeak)
1681                if sigName in varyList or not peakInstPrmMode:
1682                    sig = parmDict[sigName]
1683                    dsdA = dsdB = dsdC = 0
1684                else:
1685                    sig = G2mth.getEDsig(parmDict,pos)
1686                    dsdA,dsdB,dsdC = G2mth.getEDsigDeriv(parmDict,pos)
1687                sig = max(sig,0.001)          #avoid neg sigma
1688                Wd,fmin,fmax = getWidthsED(pos,sig)
1689                iBeg = np.searchsorted(xdata,pos-fmin)
1690                iFin = np.searchsorted(xdata,pos+fmin)
1691                if not iBeg+iFin:       #peak below low limit
1692                    iPeak += 1
1693                    continue
1694                elif not iBeg-iFin:     #peak above high limit
1695                    break
1696                dMdpk = np.zeros(shape=(4,len(xdata)))
1697                dMdipk = getdPsVoigt(pos,sig*10.**4,0.001,xdata[iBeg:iFin])
1698                dMdpk[0][iBeg:iFin] += dMdipk[0]
1699                for i in range(1,4):
1700                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1701                dervDict = {'int':dMdpk[0],'pos':-dMdpk[1],'sig':dMdpk[2]*10**4}
1702                for parmName in ['pos','int','sig']:
1703                    try:
1704                        idx = varyList.index(parmName+str(iPeak))
1705                        dMdv[idx] = dervDict[parmName]
1706                    except ValueError:
1707                        pass
1708                if 'A' in varyList:
1709                    dMdv[varyList.index('A')] += dsdA*dervDict['sig']
1710                if 'B' in varyList:
1711                    dMdv[varyList.index('B')] += dsdB*dervDict['sig']
1712                if 'C' in varyList:
1713                    dMdv[varyList.index('C')] += dsdC*dervDict['sig']
1714                iPeak += 1
1715            except KeyError:        #no more peaks to process
1716                break       
1717       
1718    elif 'B' in dataType:
1719        iPeak = 0
1720        while True:
1721            try:
1722                pos = parmDict['pos'+str(iPeak)]
1723                tth = (pos-parmDict['Zero'])
1724                intens = parmDict['int'+str(iPeak)]
1725                alpName = 'alp'+str(iPeak)
1726                if alpName in varyList or not peakInstPrmMode:
1727                    alp = parmDict[alpName]
1728                    dada0 = dada1 = 0.0
1729                else:
1730                    alp = G2mth.getPinkalpha(parmDict,tth)
1731                    dada0,dada1 = G2mth.getPinkalphaDeriv(tth)
1732                alp = max(0.0001,alp)
1733                betName = 'bet'+str(iPeak)
1734                if betName in varyList or not peakInstPrmMode:
1735                    bet = parmDict[betName]
1736                    dbdb0 = dbdb1 = 0.0
1737                else:
1738                    bet = G2mth.getPinkbeta(parmDict,tth)
1739                    dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth)
1740                bet = max(0.0001,bet)
1741                sigName = 'sig'+str(iPeak)
1742                if sigName in varyList or not peakInstPrmMode:
1743                    sig = parmDict[sigName]
1744                    dsdU = dsdV = dsdW = 0
1745                else:
1746                    sig = G2mth.getCWsig(parmDict,tth)
1747                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1748                sig = max(sig,0.001)          #avoid neg sigma
1749                gamName = 'gam'+str(iPeak)
1750                if gamName in varyList or not peakInstPrmMode:
1751                    gam = parmDict[gamName]
1752                    dgdX = dgdY = dgdZ = 0
1753                else:
1754                    gam = G2mth.getCWgam(parmDict,tth)
1755                    dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth)
1756                gam = max(gam,0.001)             #avoid neg gamma
1757                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.)
1758                iBeg = np.searchsorted(xdata,pos-fmin)
1759                iFin = np.searchsorted(xdata,pos+fmin)
1760                if not iBeg+iFin:       #peak below low limit
1761                    iPeak += 1
1762                    continue
1763                elif not iBeg-iFin:     #peak above high limit
1764                    break
1765                dMdpk = np.zeros(shape=(7,len(xdata)))
1766                dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])
1767                for i in range(1,6):
1768                    dMdpk[i][iBeg:iFin] += cw[iBeg:iFin]*intens*dMdipk[i]
1769                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1770                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.}
1771                for parmName in ['pos','int','alp','bet','sig','gam']:
1772                    try:
1773                        idx = varyList.index(parmName+str(iPeak))
1774                        dMdv[idx] = dervDict[parmName]
1775                    except ValueError:
1776                        pass
1777                if 'U' in varyList:
1778                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1779                if 'V' in varyList:
1780                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1781                if 'W' in varyList:
1782                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1783                if 'X' in varyList:
1784                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1785                if 'Y' in varyList:
1786                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1787                if 'Z' in varyList:
1788                    dMdv[varyList.index('Z')] += dgdZ*dervDict['gam']
1789                if 'alpha-0' in varyList:
1790                    dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp']
1791                if 'alpha-1' in varyList:
1792                    dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp']
1793                if 'beta-0' in varyList:
1794                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1795                if 'beta-1' in varyList:
1796                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1797                iPeak += 1
1798            except KeyError:        #no more peaks to process
1799                break       
1800    else:
1801        Pdabc = parmDict['Pdabc']
1802        difC = parmDict['difC']
1803        iPeak = 0
1804        while True:
1805            try:
1806                pos = parmDict['pos'+str(iPeak)]               
1807                tof = pos-parmDict['Zero']
1808                dsp = tof/difC
1809                intens = parmDict['int'+str(iPeak)]
1810                alpName = 'alp'+str(iPeak)
1811                if alpName in varyList or not peakInstPrmMode:
1812                    alp = parmDict[alpName]
1813                else:
1814                    if len(Pdabc):
1815                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1816                        dada0 = 0
1817                    else:
1818                        alp = G2mth.getTOFalpha(parmDict,dsp)
1819                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1820                betName = 'bet'+str(iPeak)
1821                if betName in varyList or not peakInstPrmMode:
1822                    bet = parmDict[betName]
1823                else:
1824                    if len(Pdabc):
1825                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1826                        dbdb0 = dbdb1 = dbdb2 = 0
1827                    else:
1828                        bet = G2mth.getTOFbeta(parmDict,dsp)
1829                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1830                sigName = 'sig'+str(iPeak)
1831                if sigName in varyList or not peakInstPrmMode:
1832                    sig = parmDict[sigName]
1833                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1834                else:
1835                    sig = G2mth.getTOFsig(parmDict,dsp)
1836                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1837                gamName = 'gam'+str(iPeak)
1838                if gamName in varyList or not peakInstPrmMode:
1839                    gam = parmDict[gamName]
1840                    dsdX = dsdY = dsdZ = 0
1841                else:
1842                    gam = G2mth.getTOFgamma(parmDict,dsp)
1843                    dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp)
1844                gam = max(gam,0.001)             #avoid neg gamma
1845                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1846                iBeg = np.searchsorted(xdata,pos-fmin)
1847                lenX = len(xdata)
1848                if not iBeg:
1849                    iFin = np.searchsorted(xdata,pos+fmax)
1850                elif iBeg == lenX:
1851                    iFin = iBeg
1852                else:
1853                    iFin = np.searchsorted(xdata,pos+fmax)
1854                if not iBeg+iFin:       #peak below low limit
1855                    iPeak += 1
1856                    continue
1857                elif not iBeg-iFin:     #peak above high limit
1858                    break
1859                dMdpk = np.zeros(shape=(7,len(xdata)))
1860                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1861                for i in range(1,6):
1862                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1863                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1864                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1865                for parmName in ['pos','int','alp','bet','sig','gam']:
1866                    try:
1867                        idx = varyList.index(parmName+str(iPeak))
1868                        dMdv[idx] = dervDict[parmName]
1869                    except ValueError:
1870                        pass
1871                if 'alpha' in varyList:
1872                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1873                if 'beta-0' in varyList:
1874                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1875                if 'beta-1' in varyList:
1876                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1877                if 'beta-q' in varyList:
1878                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1879                if 'sig-0' in varyList:
1880                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1881                if 'sig-1' in varyList:
1882                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1883                if 'sig-2' in varyList:
1884                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1885                if 'sig-q' in varyList:
1886                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1887                if 'X' in varyList:
1888                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1889                if 'Y' in varyList:
1890                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']
1891                if 'Z' in varyList:
1892                    dMdv[varyList.index('Z')] += dsdZ*dervDict['gam']
1893                iPeak += 1
1894            except KeyError:        #no more peaks to process
1895                break
1896    if 'BF mult' in varyList:
1897        dMdv[varyList.index('BF mult')] = fixback
1898       
1899    return dMdv
1900       
1901def Dict2Values(parmdict, varylist):
1902    '''Use before call to leastsq to setup list of values for the parameters
1903    in parmdict, as selected by key in varylist'''
1904    return [parmdict[key] for key in varylist] 
1905   
1906def Values2Dict(parmdict, varylist, values):
1907    ''' Use after call to leastsq to update the parameter dictionary with
1908    values corresponding to keys in varylist'''
1909    parmdict.update(zip(varylist,values))
1910   
1911def SetBackgroundParms(Background):
1912    'Loads background parameters into dicts/lists to create varylist & parmdict'
1913    if len(Background) == 1:            # fix up old backgrounds
1914        Background.append({'nDebye':0,'debyeTerms':[]})
1915    bakType,bakFlag = Background[0][:2]
1916    backVals = Background[0][3:]
1917    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1918    Debye = Background[1]           #also has background peaks stuff
1919    backDict = dict(zip(backNames,backVals))
1920    backVary = []
1921    if bakFlag:
1922        backVary = backNames
1923
1924    backDict['nDebye'] = Debye['nDebye']
1925    debyeDict = {}
1926    debyeList = []
1927    for i in range(Debye['nDebye']):
1928        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1929        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1930        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1931    debyeVary = []
1932    for item in debyeList:
1933        if item[1]:
1934            debyeVary.append(item[0])
1935    backDict.update(debyeDict)
1936    backVary += debyeVary
1937
1938    backDict['nPeaks'] = Debye['nPeaks']
1939    peaksDict = {}
1940    peaksList = []
1941    for i in range(Debye['nPeaks']):
1942        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1943        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1944        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1945    peaksVary = []
1946    for item in peaksList:
1947        if item[1]:
1948            peaksVary.append(item[0])
1949    backDict.update(peaksDict)
1950    backVary += peaksVary
1951    if 'background PWDR' in Background[1]:
1952        backDict['Back File'] = Background[1]['background PWDR'][0]
1953        backDict['BF mult'] = Background[1]['background PWDR'][1]
1954        if len(Background[1]['background PWDR']) > 2:
1955            if Background[1]['background PWDR'][2]:
1956                backVary += ['BF mult',]
1957    return bakType,backDict,backVary
1958   
1959def DoCalibInst(IndexPeaks,Inst):
1960   
1961    def SetInstParms():
1962        dataType = Inst['Type'][0]
1963        insVary = []
1964        insNames = []
1965        insVals = []
1966        for parm in Inst:
1967            insNames.append(parm)
1968            insVals.append(Inst[parm][1])
1969            if parm in ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']:
1970                if Inst[parm][2]:
1971                    insVary.append(parm)
1972        instDict = dict(zip(insNames,insVals))
1973        return dataType,instDict,insVary
1974       
1975    def GetInstParms(parmDict,Inst,varyList):
1976        for name in Inst:
1977            Inst[name][1] = parmDict[name]
1978       
1979    def InstPrint(Inst,sigDict):
1980        print ('Instrument Parameters:')
1981        if 'C' in Inst['Type'][0] or 'B' in Inst['Type'][0]:
1982            ptfmt = "%12.6f"
1983        else:
1984            ptfmt = "%12.3f"
1985        ptlbls = 'names :'
1986        ptstr =  'values:'
1987        sigstr = 'esds  :'
1988        for parm in Inst:
1989            if parm in  ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']:
1990                ptlbls += "%s" % (parm.center(12))
1991                ptstr += ptfmt % (Inst[parm][1])
1992                if parm in sigDict:
1993                    sigstr += ptfmt % (sigDict[parm])
1994                else:
1995                    sigstr += 12*' '
1996        print (ptlbls)
1997        print (ptstr)
1998        print (sigstr)
1999       
2000    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
2001        parmDict.update(zip(varyList,values))
2002        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
2003
2004    peakPos = []
2005    peakDsp = []
2006    peakWt = []
2007    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
2008        if peak[2] and peak[3] and sig > 0.:
2009            peakPos.append(peak[0])
2010            peakDsp.append(peak[-1])    #d-calc
2011#            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
2012            peakWt.append(1./(sig*peak[-1]))   #
2013    peakPos = np.array(peakPos)
2014    peakDsp = np.array(peakDsp)
2015    peakWt = np.array(peakWt)
2016    dataType,insDict,insVary = SetInstParms()
2017    parmDict = {}
2018    parmDict.update(insDict)
2019    varyList = insVary
2020    if not len(varyList):
2021        G2fil.G2Print ('**** ERROR - nothing to refine! ****')
2022        return False
2023    while True:
2024        begin = time.time()
2025        values =  np.array(Dict2Values(parmDict, varyList))
2026        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
2027            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
2028        ncyc = int(result[2]['nfev']/2)
2029        runtime = time.time()-begin   
2030        chisq = np.sum(result[2]['fvec']**2)
2031        Values2Dict(parmDict, varyList, result[0])
2032        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
2033        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList)))
2034        G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2035        G2fil.G2Print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF))
2036        try:
2037            sig = np.sqrt(np.diag(result[1])*GOF)
2038            if np.any(np.isnan(sig)):
2039                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
2040            break                   #refinement succeeded - finish up!
2041        except ValueError:          #result[1] is None on singular matrix
2042            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
2043       
2044    sigDict = dict(zip(varyList,sig))
2045    GetInstParms(parmDict,Inst,varyList)
2046    InstPrint(Inst,sigDict)
2047    return True
2048
2049def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[],oneCycle=False,controls=None,wtFactor=1.0,dlg=None,noFit=False):
2050    '''Called to perform a peak fit, refining the selected items in the peak
2051    table as well as selected items in the background.
2052
2053    :param str FitPgm: type of fit to perform. At present this is ignored.
2054    :param list Peaks: a list of peaks. Each peak entry is a list with 8 values:
2055      four values followed by a refine flag where the values are: position, intensity,
2056      sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List"
2057      tree entry, dict item "peaks"
2058    :param list Background: describes the background. List with two items.
2059      Item 0 specifies a background model and coefficients. Item 1 is a dict.
2060      From the Histogram/Background tree entry.
2061    :param list Limits: min and max x-value to use
2062    :param dict Inst: Instrument parameters
2063    :param dict Inst2: more Instrument parameters
2064    :param numpy.array data: a 5xn array. data[0] is the x-values,
2065      data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are
2066      calc, background and difference intensities, respectively.
2067    :param array fixback: fixed background array; same size as data[0-5]
2068    :param list prevVaryList: Used in sequential refinements to override the
2069      variable list. Defaults as an empty list.
2070    :param bool oneCycle: True if only one cycle of fitting should be performed
2071    :param dict controls: a dict specifying two values, Ftol = controls['min dM/M']
2072      and derivType = controls['deriv type']. If None default values are used.
2073    :param float wtFactor: weight multiplier; = 1.0 by default
2074    :param wx.Dialog dlg: A dialog box that is updated with progress from the fit.
2075      Defaults to None, which means no updates are done.
2076    :param bool noFit: When noFit is True, a refinement is not performed. Default
2077      is False.
2078
2079    '''
2080    def GetBackgroundParms(parmList,Background):
2081        iBak = 0
2082        while True:
2083            try:
2084                bakName = 'Back;'+str(iBak)
2085                Background[0][iBak+3] = parmList[bakName]
2086                iBak += 1
2087            except KeyError:
2088                break
2089        iDb = 0
2090        while True:
2091            names = ['DebyeA;','DebyeR;','DebyeU;']
2092            try:
2093                for i,name in enumerate(names):
2094                    val = parmList[name+str(iDb)]
2095                    Background[1]['debyeTerms'][iDb][2*i] = val
2096                iDb += 1
2097            except KeyError:
2098                break
2099        iDb = 0
2100        while True:
2101            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
2102            try:
2103                for i,name in enumerate(names):
2104                    val = parmList[name+str(iDb)]
2105                    Background[1]['peaksList'][iDb][2*i] = val
2106                iDb += 1
2107            except KeyError:
2108                break
2109        if 'BF mult' in parmList:
2110            Background[1]['background PWDR'][1] = parmList['BF mult']
2111               
2112    def BackgroundPrint(Background,sigDict):
2113        print ('Background coefficients for '+Background[0][0]+' function')
2114        ptfmt = "%12.5f"
2115        ptstr =  'value: '
2116        sigstr = 'esd  : '
2117        for i,back in enumerate(Background[0][3:]):
2118            ptstr += ptfmt % (back)
2119            if Background[0][1]:
2120                prm = 'Back;'+str(i)
2121                if prm in sigDict:
2122                    sigstr += ptfmt % (sigDict[prm])
2123                else:
2124                    sigstr += " "*12
2125            if len(ptstr) > 75:
2126                print (ptstr)
2127                if Background[0][1]: print (sigstr)
2128                ptstr =  'value: '
2129                sigstr = 'esd  : '
2130        if len(ptstr) > 8:
2131            print (ptstr)
2132            if Background[0][1]: print (sigstr)
2133
2134        if Background[1]['nDebye']:
2135            parms = ['DebyeA;','DebyeR;','DebyeU;']
2136            print ('Debye diffuse scattering coefficients')
2137            ptfmt = "%12.5f"
2138            print (' term       DebyeA       esd        DebyeR       esd        DebyeU        esd')
2139            for term in range(Background[1]['nDebye']):
2140                line = ' term %d'%(term)
2141                for ip,name in enumerate(parms):
2142                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
2143                    if name+str(term) in sigDict:
2144                        line += ptfmt%(sigDict[name+str(term)])
2145                    else:
2146                        line += " "*12
2147                print (line)
2148        if Background[1]['nPeaks']:
2149            print ('Coefficients for Background Peaks')
2150            ptfmt = "%15.3f"
2151            for j,pl in enumerate(Background[1]['peaksList']):
2152                names =  'peak %3d:'%(j+1)
2153                ptstr =  'values  :'
2154                sigstr = 'esds    :'
2155                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
2156                    val = pl[2*i]
2157                    prm = lbl+";"+str(j)
2158                    names += '%15s'%(prm)
2159                    ptstr += ptfmt%(val)
2160                    if prm in sigDict:
2161                        sigstr += ptfmt%(sigDict[prm])
2162                    else:
2163                        sigstr += " "*15
2164                print (names)
2165                print (ptstr)
2166                print (sigstr)
2167        if 'BF mult' in sigDict:
2168            print('Background file mult: %.3f(%d)'%(Background[1]['background PWDR'][1],int(1000*sigDict['BF mult'])))
2169                           
2170    def SetInstParms(Inst):
2171        dataType = Inst['Type'][0]
2172        insVary = []
2173        insNames = []
2174        insVals = []
2175        for parm in Inst:
2176            insNames.append(parm)
2177            insVals.append(Inst[parm][1])
2178            if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C',
2179                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]:
2180                    insVary.append(parm)
2181        instDict = dict(zip(insNames,insVals))
2182        if 'SH/L' in instDict:
2183            instDict['SH/L'] = max(instDict['SH/L'],0.002)
2184        return dataType,instDict,insVary
2185       
2186    def GetInstParms(parmDict,Inst,varyList,Peaks):
2187        for name in Inst:
2188            Inst[name][1] = parmDict[name]
2189        iPeak = 0
2190        while True:
2191            try:
2192                sigName = 'sig'+str(iPeak)
2193                pos = parmDict['pos'+str(iPeak)]
2194                if sigName not in varyList and peakInstPrmMode:
2195                    if 'T' in Inst['Type'][0]:
2196                        dsp = G2lat.Pos2dsp(Inst,pos)
2197                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
2198                    if 'E' in Inst['Type'][0]:
2199                        parmDict[sigName] = G2mth.getEDsig(parmDict,pos)
2200                    else:
2201                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
2202                gamName = 'gam'+str(iPeak)
2203                if gamName not in varyList and peakInstPrmMode:
2204                    if 'T' in Inst['Type'][0]:
2205                        dsp = G2lat.Pos2dsp(Inst,pos)
2206                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
2207                    else:
2208                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
2209                iPeak += 1
2210            except KeyError:
2211                break
2212       
2213    def InstPrint(Inst,sigDict):
2214        print ('Instrument Parameters:')
2215        ptfmt = "%12.6f"
2216        ptlbls = 'names :'
2217        ptstr =  'values:'
2218        sigstr = 'esds  :'
2219        for parm in Inst:
2220            if parm in  ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C',
2221                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']:
2222                ptlbls += "%s" % (parm.center(12))
2223                ptstr += ptfmt % (Inst[parm][1])
2224                if parm in sigDict:
2225                    sigstr += ptfmt % (sigDict[parm])
2226                else:
2227                    sigstr += 12*' '
2228        print (ptlbls)
2229        print (ptstr)
2230        print (sigstr)
2231
2232    def SetPeaksParms(dataType,Peaks):
2233        peakNames = []
2234        peakVary = []
2235        peakVals = []
2236        if 'C' in dataType:
2237            names = ['pos','int','sig','gam']
2238        elif 'E' in dataType:
2239            names = ['pos','int','sig']
2240        else:   #'T' and 'B'
2241            names = ['pos','int','alp','bet','sig','gam']
2242        for i,peak in enumerate(Peaks):
2243            for j,name in enumerate(names):
2244                peakVals.append(peak[2*j])
2245                parName = name+str(i)
2246                peakNames.append(parName)
2247                if peak[2*j+1]:
2248                    peakVary.append(parName)
2249        return dict(zip(peakNames,peakVals)),peakVary
2250               
2251    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
2252        if 'C' in Inst['Type'][0]:
2253            names = ['pos','int','sig','gam']
2254        elif 'E' in Inst['Type'][0]:
2255            names = ['pos','int','sig']
2256        else:   #'T' & 'B'
2257            names = ['pos','int','alp','bet','sig','gam']
2258        for i,peak in enumerate(Peaks):
2259            pos = parmDict['pos'+str(i)]
2260            if 'difC' in Inst:
2261                dsp = pos/Inst['difC'][1]
2262            for j in range(len(names)):
2263                parName = names[j]+str(i)
2264                if parName in varyList or not peakInstPrmMode:
2265                    peak[2*j] = parmDict[parName]
2266                elif 'alp' in parName:
2267                    if 'T' in Inst['Type'][0]:
2268                        peak[2*j] = G2mth.getTOFalpha(parmDict,dsp)
2269                    else: #'B'
2270                        peak[2*j] = G2mth.getPinkalpha(parmDict,pos)
2271                elif 'bet' in parName:
2272                    if 'T' in Inst['Type'][0]:
2273                        peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
2274                    else:   #'B'
2275                        peak[2*j] = G2mth.getPinkbeta(parmDict,pos)
2276                elif 'sig' in parName:
2277                    if 'T' in Inst['Type'][0]:
2278                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
2279                    elif 'E' in Inst['Type'][0]:
2280                        peak[2*j] = G2mth.getEDsig(parmDict,pos)
2281                    else:   #'C' & 'B'
2282                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
2283                elif 'gam' in parName:
2284                    if 'T' in Inst['Type'][0]:
2285                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
2286                    else:   #'C' & 'B'
2287                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
2288                       
2289    def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW):
2290        print ('Peak coefficients:')
2291        if 'C' in dataType:
2292            names = ['pos','int','sig','gam']
2293        elif 'E' in dataType:
2294            names = ['pos','int','sig']
2295        else:   #'T' & 'B'
2296            names = ['pos','int','alp','bet','sig','gam']           
2297        head = 13*' '
2298        for name in names:
2299            if name in ['alp','bet']:
2300                head += name.center(8)+'esd'.center(8)
2301            else:
2302                head += name.center(10)+'esd'.center(10)
2303        head += 'bins'.center(8)
2304        print (head)
2305        if 'C' in dataType:
2306            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
2307        elif 'E' in dataType:
2308            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f"}
2309        elif 'T' in dataType:
2310            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
2311        else: #'B'
2312            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'alp':"%8.2f",'bet':"%8.4f",'sig':"%10.3f",'gam':"%10.3f"}
2313        for i,peak in enumerate(Peaks):
2314            ptstr =  ':'
2315            for j in range(len(names)):
2316                name = names[j]
2317                parName = name+str(i)
2318                ptstr += ptfmt[name] % (parmDict[parName])
2319                if parName in varyList:
2320                    ptstr += ptfmt[name] % (sigDict[parName])
2321                else:
2322                    if name in ['alp','bet']:
2323                        ptstr += 8*' '
2324                    else:
2325                        ptstr += 10*' '
2326            ptstr += '%9.2f'%(ptsperFW[i])
2327            print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr)
2328               
2329    def devPeakProfile(values,xdata,ydata,fixback, weights,dataType,parmdict,varylist,bakType,dlg):
2330        '''Computes a matrix where each row is the derivative of the calc-obs
2331        values (see :func:`errPeakProfile`) with respect to each parameter
2332        in backVary,insVary,peakVary. Used for peak fitting.
2333        '''
2334        parmdict.update(zip(varylist,values))
2335        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,fixback,varylist,bakType)
2336           
2337    def errPeakProfile(values,xdata,ydata,fixback,weights,dataType,parmdict,varylist,bakType,dlg):
2338        '''Computes a vector with the weighted calc-obs values differences
2339        for peak fitting
2340        '''
2341        parmdict.update(zip(varylist,values))
2342        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,fixback,varylist,bakType)-ydata)
2343        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
2344        if dlg:
2345            dlg.Raise()
2346            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
2347            if not GoOn:
2348                return -M           #abort!!
2349        return M
2350
2351    #---- beginning of DoPeakFit ----------------------------------------------
2352    if controls:
2353        Ftol = controls['min dM/M']
2354    else:
2355        Ftol = 0.0001
2356    if oneCycle:
2357        Ftol = 1.0
2358    x,y,w,yc,yb,yd = data   #these are numpy arrays - remove masks!
2359    if fixback is None:
2360        fixback = np.zeros_like(y)
2361    yc *= 0.                            #set calcd ones to zero
2362    yb *= 0.
2363    yd *= 0.
2364    xBeg = np.searchsorted(x,Limits[0])
2365    xFin = np.searchsorted(x,Limits[1])+1
2366    bakType,bakDict,bakVary = SetBackgroundParms(Background)
2367    dataType,insDict,insVary = SetInstParms(Inst)
2368    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
2369    parmDict = {}
2370    parmDict.update(bakDict)
2371    parmDict.update(insDict)
2372    parmDict.update(peakDict)
2373    parmDict['Pdabc'] = []      #dummy Pdabc
2374    parmDict.update(Inst2)      #put in real one if there
2375    if prevVaryList:
2376        varyList = prevVaryList[:]
2377    else:
2378        varyList = bakVary+insVary+peakVary
2379    fullvaryList = varyList[:]
2380    if not peakInstPrmMode:
2381        for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1','A','B','C',
2382            'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',):
2383            if v in varyList:
2384                raise Exception('Instrumental profile terms cannot be varied '+
2385                                    'after setPeakInstPrmMode(False) is used')
2386    while not noFit:
2387        begin = time.time()
2388        values =  np.array(Dict2Values(parmDict, varyList))
2389        Rvals = {}
2390        badVary = []
2391        result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
2392               args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
2393        ncyc = int(result[2]['nfev']/2)
2394        runtime = time.time()-begin   
2395        chisq = np.sum(result[2]['fvec']**2)
2396        Values2Dict(parmDict, varyList, result[0])
2397        Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactor*w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
2398        Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
2399        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList)))
2400        if ncyc:
2401            G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
2402        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
2403        sig = [0]*len(varyList)
2404        if len(varyList) == 0: break  # if nothing was refined
2405        try:
2406            sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
2407            if np.any(np.isnan(sig)):
2408                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***')
2409            break                   #refinement succeeded - finish up!
2410        except ValueError:          #result[1] is None on singular matrix
2411            G2fil.G2Print ('**** Refinement failed - singular matrix ****')
2412            Ipvt = result[2]['ipvt']
2413            for i,ipvt in enumerate(Ipvt):
2414                if not np.sum(result[2]['fjac'],axis=1)[i]:
2415                    G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
2416                    badVary.append(varyList[ipvt-1])
2417                    del(varyList[ipvt-1])
2418                    break
2419            else: # nothing removed
2420                break
2421    if dlg: dlg.Destroy()
2422    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin],fixback[xBeg:xFin])[0]
2423    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],fixback[xBeg:xFin],varyList,bakType)
2424    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
2425    if noFit:
2426        return
2427    sigDict = dict(zip(varyList,sig))
2428    GetBackgroundParms(parmDict,Background)
2429    if bakVary: BackgroundPrint(Background,sigDict)
2430    GetInstParms(parmDict,Inst,varyList,Peaks)
2431    if insVary: InstPrint(Inst,sigDict)
2432    GetPeaksParms(Inst,parmDict,Peaks,varyList)
2433    binsperFWHM = []
2434    for peak in Peaks:
2435        FWHM = getFWHM(peak[0],Inst)
2436        try:
2437            xpk = x.searchsorted(peak[0])
2438            cw = x[xpk]-x[xpk-1]
2439            binsperFWHM.append(FWHM/cw)
2440        except IndexError:
2441            binsperFWHM.append(0.)
2442    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM)
2443    if len(binsperFWHM):
2444        if min(binsperFWHM) < 1.:
2445            G2fil.G2Print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***')
2446            if 'T' in Inst['Type'][0]:
2447                G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters')
2448            else:
2449                G2fil.G2Print (' Manually increase W in Instrument Parameters')
2450        elif min(binsperFWHM) < 4.:
2451            G2fil.G2Print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***')
2452            G2fil.G2Print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM)))
2453    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
2454   
2455def calcIncident(Iparm,xdata):
2456    'needs a doc string'
2457
2458    def IfunAdv(Iparm,xdata):
2459        Itype = Iparm['Itype']
2460        Icoef = Iparm['Icoeff']
2461        DYI = np.ones((12,xdata.shape[0]))
2462        YI = np.ones_like(xdata)*Icoef[0]
2463       
2464        x = xdata/1000.                 #expressions are in ms
2465        if Itype == 'Exponential':
2466            for i in [1,3,5,7,9]:
2467                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2468                YI += Icoef[i]*Eterm
2469                DYI[i] *= Eterm
2470                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
2471        elif 'Maxwell'in Itype:
2472            Eterm = np.exp(-Icoef[2]/x**2)
2473            DYI[1] = Eterm/x**5
2474            DYI[2] = -Icoef[1]*DYI[1]/x**2
2475            YI += (Icoef[1]*Eterm/x**5)
2476            if 'Exponential' in Itype:
2477                for i in range(3,11,2):
2478                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
2479                    YI += Icoef[i]*Eterm
2480                    DYI[i] *= Eterm
2481                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
2482            else:   #Chebyschev
2483                T = (2./x)-1.
2484                Ccof = np.ones((12,xdata.shape[0]))
2485                Ccof[1] = T
2486                for i in range(2,12):
2487                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
2488                for i in range(1,10):
2489                    YI += Ccof[i]*Icoef[i+2]
2490                    DYI[i+2] =Ccof[i]
2491        return YI,DYI
2492       
2493    Iesd = np.array(Iparm['Iesd'])
2494    Icovar = Iparm['Icovar']
2495    YI,DYI = IfunAdv(Iparm,xdata)
2496    YI = np.where(YI>0,YI,1.)
2497    WYI = np.zeros_like(xdata)
2498    vcov = np.zeros((12,12))
2499    k = 0
2500    for i in range(12):
2501        for j in range(i,12):
2502            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
2503            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
2504            k += 1
2505    M = np.inner(vcov,DYI.T)
2506    WYI = np.sum(M*DYI,axis=0)
2507    WYI = np.where(WYI>0.,WYI,0.)
2508    return YI,WYI
2509
2510#### RMCutilities ################################################################################
2511def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd):
2512    inst = PWDdata['Instrument Parameters'][0]
2513    Xsb = 0.
2514    Ysb = 0.
2515    if 'T' in inst['Type'][1]:
2516        difC = inst['difC'][1]
2517        if useSamBrd[0]:
2518            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2519                Xsb = 1.e-4*difC/Size[1][0]
2520        if useSamBrd[1]:
2521            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2522                Ysb = 1.e-6*difC*Mustrain[1][0]
2523        prms = ['Bank',
2524                'difC','difA','Zero','2-theta',
2525                'alpha','beta-0','beta-1',
2526                'sig-0','sig-1','sig-2',
2527                'Z','X','Y']
2528        fname = Name+'.inst'
2529        fl = open(fname,'w')
2530        fl.write('1\n')
2531        fl.write('%d\n'%int(inst[prms[0]][1]))
2532        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]))
2533        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[5]][1],inst[prms[6]][1],inst[prms[7]][1]))
2534        fl.write('%12.6e%14.6e%14.6e\n'%(inst[prms[8]][1],inst[prms[9]][1],inst[prms[10]][1]))   
2535        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))
2536        fl.write('\n\n\n')
2537        fl.close()
2538    else:
2539        if useSamBrd[0]:
2540            wave = G2mth.getWave(inst)
2541            if 'ellipsoidal' not in Size[0]:    #take the isotropic term only
2542                Xsb = 1.8*wave/(np.pi*Size[1][0])
2543        if useSamBrd[1]:
2544            if 'generalized' not in Mustrain[0]:    #take the isotropic term only
2545                Ysb = 0.0180*Mustrain[1][0]/np.pi
2546        prms = ['Bank',
2547                'Lam','Zero','Polariz.',
2548                'U','V','W',
2549                'X','Y']
2550        fname = Name+'.inst'
2551        fl = open(fname,'w')
2552        fl.write('1\n')
2553        fl.write('%d\n'%int(inst[prms[0]][1]))
2554        fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],-100.*inst[prms[2]][1],inst[prms[3]][1],0))
2555        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1]))
2556        fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0))   
2557        fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0))
2558        fl.write('\n\n\n')
2559        fl.close()
2560    return fname
2561   
2562def MakeBack(PWDdata,Name):
2563    Back = PWDdata['Background'][0]
2564    inst = PWDdata['Instrument Parameters'][0]
2565    if 'chebyschev-1' != Back[0]:
2566        return None
2567    Nback = Back[2]
2568    BackVals = Back[3:]
2569    fname = Name+'.back'
2570    fl = open(fname,'w')
2571    fl.write('%10d\n'%Nback)
2572    for val in BackVals:
2573        if 'T' in inst['Type'][1]:
2574            fl.write('%12.6g\n'%(float(val)))
2575        else:
2576            fl.write('%12.6g\n'%val)
2577    fl.close()
2578    return fname
2579
2580def findDup(Atoms):
2581    Dup = []
2582    Fracs = []
2583    for iat1,at1 in enumerate(Atoms):
2584        if any([at1[0] in dup for dup in Dup]):
2585            continue
2586        else:
2587            Dup.append([at1[0],])
2588            Fracs.append([at1[6],])
2589        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2590            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2591                Dup[-1] += [at2[0],]
2592                Fracs[-1]+= [at2[6],]
2593    return Dup,Fracs
2594
2595def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):   
2596   
2597    Meta = RMCPdict['metadata']
2598    Atseq = RMCPdict['atSeq']
2599    Supercell =  RMCPdict['SuperCell']
2600    generalData = Phase['General']
2601    Dups,Fracs = findDup(Phase['Atoms'])
2602    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2603    ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2604    Sample = PWDdata['Sample Parameters']
2605    Meta['temperature'] = Sample['Temperature']
2606    Meta['pressure'] = Sample['Pressure']
2607    Cell = generalData['Cell'][1:7]
2608    Trans = np.eye(3)*np.array(Supercell)
2609    newPhase = copy.deepcopy(Phase)
2610    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2611    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2612    GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2613    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.
2614    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2615    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2616    Natm = np.count_nonzero(Natm-1)
2617    Atoms = newPhase['Atoms']
2618    reset = False
2619   
2620    if ifSfracs:
2621        Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2622        Natm = np.count_nonzero(Natm-1)
2623        Satoms = []
2624        for i in range(len(Atoms)//Natm):
2625            ind = i*Natm
2626            Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2627        Natoms = []
2628        for satoms in Satoms:
2629            for idup,dup in enumerate(Dups):
2630                ldup = len(dup)
2631                natm = len(satoms)
2632                i = 0
2633                while i < natm:
2634                    if satoms[i][0] in dup:
2635                        atoms = satoms[i:i+ldup]
2636                        try:
2637                            atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2638                            Natoms.append(atom)
2639                        except IndexError:      #what about vacancies?
2640                            if 'Va' not in Atseq:
2641                                reset = True
2642                                Atseq.append('Va')
2643                                RMCPdict['aTypes']['Va'] = 0.0
2644                            atom = atoms[0]
2645                            atom[1] = 'Va'
2646                            Natoms.append(atom)
2647                        i += ldup
2648                    else:
2649                       i += 1
2650    else:
2651        Natoms = Atoms
2652   
2653    NAtype = np.zeros(len(Atseq))
2654    for atom in Natoms:
2655        NAtype[Atseq.index(atom[1])] += 1
2656    NAstr = ['%6d'%i for i in NAtype]
2657    Cell = newPhase['General']['Cell'][1:7]
2658    if os.path.exists(Name+'.his6f'):
2659        os.remove(Name+'.his6f')
2660    if os.path.exists(Name+'.neigh'):
2661        os.remove(Name+'.neigh')
2662    fname = Name+'.rmc6f'
2663    fl = open(fname,'w')
2664    fl.write('(Version 6f format configuration file)\n')
2665    for item in Meta:
2666        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2667    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2668    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2669    fl.write('Number of atoms:                %d\n'%len(Natoms))
2670    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2671    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2672            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2673    A,B = G2lat.cell2AB(Cell,True)
2674    fl.write('Lattice vectors (Ang):\n')   
2675    for i in [0,1,2]:
2676        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2677    fl.write('Atoms (fractional coordinates):\n')
2678    nat = 0
2679    for atm in Atseq:
2680        for iat,atom in enumerate(Natoms):
2681            if atom[1] == atm:
2682                nat += 1
2683                atcode = Atcodes[iat].split(':')
2684                cell = [0,0,0]
2685                if '+' in atcode[1]:
2686                    cell = eval(atcode[1].split('+')[1])
2687                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2688                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2689    fl.close()
2690    return fname,reset
2691
2692def MakeBragg(PWDdata,Name,Phase):
2693    generalData = Phase['General']
2694    Vol = generalData['Cell'][7]
2695    Data = PWDdata['Data']
2696    Inst = PWDdata['Instrument Parameters'][0]
2697    Bank = int(Inst['Bank'][1])
2698    Sample = PWDdata['Sample Parameters']
2699    Scale = Sample['Scale'][0]
2700    Limits = PWDdata['Limits'][1]
2701    Ibeg = np.searchsorted(Data[0],Limits[0])
2702    Ifin = np.searchsorted(Data[0],Limits[1])+1
2703    fname = Name+'.bragg'
2704    fl = open(fname,'w')
2705    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2706    if 'T' in Inst['Type'][0]:
2707        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2708        for i in range(Ibeg,Ifin-1):
2709            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2710    else:
2711        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2712        for i in range(Ibeg,Ifin-1):
2713            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2714    fl.close()
2715    return fname
2716
2717def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2718    Meta = RMCPdict['metadata']
2719    Times = RMCPdict['runTimes']
2720    Atseq = RMCPdict['atSeq']
2721    Natoms = RMCPdict['NoAtoms']
2722    sumatms = np.sum(np.array([Natoms[iatm] for iatm in Natoms]))
2723    Isotope = RMCPdict['Isotope']
2724    Isotopes = RMCPdict['Isotopes']
2725    Atypes = RMCPdict['aTypes']
2726    if 'Va' in Atypes:
2727        Isotope['Va'] = 'Nat. Abund.'
2728        Isotopes['Va'] = {'Nat. Abund.':{'SL':[0.0,0.0]}}
2729    atPairs = RMCPdict['Pairs']
2730    Files = RMCPdict['files']
2731    BraggWt = RMCPdict['histogram'][1]
2732    inst = PWDdata['Instrument Parameters'][0]
2733    try:
2734        refList = PWDdata['Reflection Lists'][Name]['RefList']
2735    except KeyError:
2736        return 'Error - missing reflection list; you must do Refine first'
2737    dMin = refList[-1][4]
2738    gsasType = 'xray2'
2739    if 'T' in inst['Type'][1]:
2740        gsasType = 'gsas3'
2741    elif 'X' in inst['Type'][1]:
2742        XFF = G2elem.GetFFtable(Atseq)
2743        Xfl = open(Name+'.xray','w')
2744        for atm in Atseq:
2745            fa = XFF[atm]['fa']
2746            fb = XFF[atm]['fb']
2747            fc = XFF[atm]['fc']
2748            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2749                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2750        Xfl.close()
2751    lenA = len(Atseq)
2752    Pairs = []
2753    Ncoeff = []
2754    Nblen = [Isotopes[at][Isotope[at]]['SL'][0] for at in Atypes]
2755    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2756        Pairs += pair
2757    for pair in Pairs:
2758        pair = pair.replace(' ','')
2759        at1,at2 = pair.split('-')
2760        if at1 == 'Va' or at2 == 'Va':
2761            ncoef = 0.0
2762        else:
2763            ncoef = Isotopes[at1][Isotope[at1]]['SL'][0]*Natoms[at1]/sumatms
2764            ncoef *= Isotopes[at2][Isotope[at2]]['SL'][0]*Natoms[at2]/sumatms
2765        if at1 != at2:
2766            ncoef *= 2.
2767        Ncoeff += [ncoef,]
2768    pairMin = [atPairs[pair] if pair in atPairs else [0.0,0.,0.] for pair in Pairs ]
2769    maxMoves = [Atypes[atm] if atm in Atypes else 0.0 for atm in Atseq ]
2770    fname = Name+'.dat'
2771    fl = open(fname,'w')
2772    fl.write(' %% Hand edit the following as needed\n')
2773    fl.write('TITLE :: '+Name+'\n')
2774    fl.write('MATERIAL :: '+Meta['material']+'\n')
2775    fl.write('PHASE :: '+Meta['phase']+'\n')
2776    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2777    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2778    if RMCPdict.get('useGPU',False):
2779        fl.write('GPU_ACCELERATOR :: 0\n')
2780    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2781    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2782    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2783    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2784    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2785    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2786    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2787    fl.write('PRINT_PERIOD :: 100\n')
2788    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2789    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2790    fl.write('\n')
2791    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2792    fl.write('\n')
2793    fl.write('FLAGS ::\n')
2794    fl.write('  > NO_MOVEOUT\n')
2795    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2796    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2797    fl.write('\n')
2798    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2799    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2800    fl.write('IGNORE_HISTORY_FILE ::\n')
2801    fl.write('\n')
2802    if 'T' in inst['Type'][1]:
2803        fl.write('NEUTRON_COEFFICIENTS :: '+''.join(['%9.5f'%coeff for coeff in Ncoeff])+'\n')
2804    fl.write('DISTANCE_WINDOW ::\n')
2805    fl.write('  > MNDIST :: %s\n'%minD)
2806    fl.write('  > MXDIST :: %s\n'%maxD)
2807    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2808        fl.write('\n')
2809        fl.write('POTENTIALS ::\n')
2810        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2811        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2812        if len(RMCPdict['Potentials']['Stretch']):
2813            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2814            for bond in RMCPdict['Potentials']['Stretch']:
2815                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2816        if len(RMCPdict['Potentials']['Angles']):
2817            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2818            for angle in RMCPdict['Potentials']['Angles']:
2819                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2820                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2821    if RMCPdict['useBVS']:
2822        fl.write('BVS ::\n')
2823        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2824        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2825        oxid = []
2826        for val in RMCPdict['Oxid']:
2827            if len(val) == 3:
2828                oxid.append(val[0][1:])
2829            else:
2830                oxid.append(val[0][2:])
2831        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2832        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2833        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2834        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))       
2835        fl.write('  > SAVE :: 100000\n')
2836        fl.write('  > UPDATE :: 100000\n')
2837        if len(RMCPdict['Swap']):
2838            fl.write('\n')
2839            fl.write('SWAP_MULTI ::\n')
2840            for swap in RMCPdict['Swap']:
2841                try:
2842                    at1 = Atseq.index(swap[0])
2843                    at2 = Atseq.index(swap[1])
2844                except ValueError:
2845                    break
2846                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2847       
2848    if len(RMCPdict['FxCN']):
2849        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2850        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2851            try:
2852                at1 = Atseq.index(fxcn[0])
2853                at2 = Atseq.index(fxcn[1])
2854            except ValueError:
2855                break
2856            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]))
2857    if len(RMCPdict['AveCN']):
2858        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2859        for iav,avcn in enumerate(RMCPdict['AveCN']):
2860            try:
2861                at1 = Atseq.index(avcn[0])
2862                at2 = Atseq.index(avcn[1])
2863            except ValueError:
2864                break
2865            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]))
2866    for File in Files:
2867        if Files[File][0] and Files[File][0] != 'Select':
2868            if 'Xray' in File and 'F(Q)' in File:
2869                fqdata = open(Files[File][0],'r')
2870                lines = int(fqdata.readline()[:-1])
2871                fqdata.close()
2872            fl.write('\n')
2873            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2874            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2875            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2876            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2877            if 'Xray' not in File:
2878                fl.write('  > START_POINT :: 1\n')
2879                fl.write('  > END_POINT :: 3000\n')
2880                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2881            fl.write('  > CONSTANT_OFFSET 0.000\n')
2882            fl.write('  > NO_FITTED_OFFSET\n')
2883            if RMCPdict['FitScale']:
2884                fl.write('  > FITTED_SCALE\n')
2885            else:
2886                fl.write('  > NO_FITTED_SCALE\n')
2887            if Files[File][3] !='RMC':
2888                fl.write('  > %s\n'%Files[File][3])
2889            if 'reciprocal' in File:
2890                fl.write('  > CONVOLVE ::\n')
2891                if 'Xray' in File:
2892                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2893                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2894                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2895                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2896    fl.write('\n')
2897    fl.write('BRAGG ::\n')
2898    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2899    fl.write('  > RECALCUATE\n')
2900    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2901    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2902    if 'T' in inst['Type'][1]:
2903        fl.write('  > SCATTERING LENGTH :: '+''.join(['%8.4f'%blen for blen in Nblen])+'\n')
2904    fl.write('\n')
2905    fl.write('END  ::\n')
2906    fl.close()
2907    return fname
2908
2909# def FindBonds(Phase,RMCPdict):
2910#     generalData = Phase['General']
2911#     cx,ct,cs,cia = generalData['AtomPtrs']
2912#     atomData = Phase['Atoms']
2913#     Res = 'RMC'
2914#     if 'macro' in generalData['Type']:
2915#         Res = atomData[0][ct-3]
2916#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2917#     Pairs = RMCPdict['Pairs']   #dict!
2918#     BondList = []
2919#     notNames = []
2920#     for FrstName in AtDict:
2921#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2922#         Atyp1 = AtDict[FrstName]
2923#         if 'Va' in Atyp1:
2924#             continue
2925#         for nbr in nbrs:
2926#             Atyp2 = AtDict[nbr[0]]
2927#             if 'Va' in Atyp2:
2928#                 continue
2929#             try:
2930#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2931#             except KeyError:
2932#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2933#             if any(bndData):
2934#                 if bndData[0] <= nbr[1] <= bndData[1]:
2935#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2936#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2937#                     if bondStr not in BondList and revbondStr not in BondList:
2938#                         BondList.append(bondStr)
2939#         notNames.append(FrstName)
2940#     return Res,BondList
2941
2942# def FindAngles(Phase,RMCPdict):
2943#     generalData = Phase['General']
2944#     Cell = generalData['Cell'][1:7]
2945#     Amat = G2lat.cell2AB(Cell)[0]
2946#     cx,ct,cs,cia = generalData['AtomPtrs']
2947#     atomData = Phase['Atoms']
2948#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2949#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2950#     Angles = RMCPdict['Angles']
2951#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2952#     AngleList = []
2953#     for MidName in AtDict:
2954#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2955#         if len(nbrs) < 2: #need 2 neighbors to make an angle
2956#             continue
2957#         Atyp2 = AtDict[MidName]
2958#         for i,nbr1 in enumerate(nbrs):
2959#             Atyp1 = AtDict[nbr1[0]]
2960#             for j,nbr3 in enumerate(nbrs[i+1:]):
2961#                 Atyp3 = AtDict[nbr3[0]]
2962#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2963#                 try:
2964#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2965#                 except KeyError:
2966#                     try:
2967#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2968#                     except KeyError:
2969#                         continue
2970#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2971#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
2972#                 if angData[0] <= calAngle <= angData[1]:
2973#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
2974#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
2975#                     if angStr not in AngleList and revangStr not in AngleList:
2976#                         AngleList.append(angStr)
2977#     return AngleList
2978
2979# def GetSqConvolution(XY,d):
2980
2981#     n = XY.shape[1]
2982#     snew = np.zeros(n)
2983#     dq = np.zeros(n)
2984#     sold = XY[1]
2985#     q = XY[0]
2986#     dq[1:] = np.diff(q)
2987#     dq[0] = dq[1]
2988   
2989#     for j in range(n):
2990#         for i in range(n):
2991#             b = abs(q[i]-q[j])
2992#             t = q[i]+q[j]
2993#             if j == i:
2994#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
2995#             else:
2996#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
2997#         snew[j] /= np.pi*q[j]
2998   
2999#     snew[0] = snew[1]
3000#     return snew
3001
3002# def GetMaxSphere(pdbName):
3003#     try:
3004#         pFil = open(pdbName,'r')
3005#     except FileNotFoundError:
3006#         return None
3007#     while True:
3008#         line = pFil.readline()
3009#         if 'Boundary' in line:
3010#             line = line.split()[3:]
3011#             G = np.array([float(item) for item in line])
3012#             G = np.reshape(G,(3,3))**2
3013#             G = nl.inv(G)
3014#             pFil.close()
3015#             break
3016#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
3017#     return min(dspaces)
3018
3019def findfullrmc():
3020    '''Find where fullrmc is installed. Tries the following:
3021   
3022         1. Returns the Config var 'fullrmc_exec', if defined. No check
3023            is done that the interpreter has fullrmc
3024         2. The current Python interpreter if fullrmc can be imported
3025            and fullrmc is version 5+
3026         3. The path is checked for a fullrmc image as named by Bachir
3027
3028    :returns: the full path to a python executable that is assumed to
3029      have fullrmc installed or None, if it was not found.
3030    '''
3031    is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
3032    if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
3033            GSASIIpath.GetConfigValue('fullrmc_exec')):
3034        return GSASIIpath.GetConfigValue('fullrmc_exec')
3035    try:
3036        import fullrmc
3037        if int(fullrmc.__version__.split('.')[0]) >= 5:
3038            return sys.executable
3039    except:
3040        pass
3041    pathlist = os.environ["PATH"].split(os.pathsep)
3042    for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
3043                  os.path.split(sys.executable)[0]):
3044        if p not in pathlist: pathlist.insert(0,p)
3045    import glob
3046    for p in pathlist:
3047        if sys.platform == "darwin":
3048            lookfor = "fullrmc*macOS*i386-64bit"
3049        elif sys.platform == "win32":
3050            lookfor = "fullrmc*.exe"
3051        else:
3052            lookfor = "fullrmc*"
3053        fl = glob.glob(lookfor)
3054        if len(fl) > 0:
3055            return os.path.abspath(sorted(fl)[0])
3056       
3057def findPDFfit():
3058    '''Find if PDFfit2 is installed (may be local to GSAS-II). Does the following:
3059    :returns: the full path to a python executable or None, if it was not found.
3060    '''
3061    try:
3062        sys.path.append('%s'%GSASIIpath.path2GSAS2) 
3063        from diffpy.pdffit2 import PdfFit
3064        return sys.executable
3065    except:
3066        return None
3067   
3068def GetPDFfitAtomVar(Phase,RMCPdict):
3069    ''' Find dict of independent "@n" variables for PDFfit in atom constraints
3070    '''
3071    General = Phase['General']
3072    Atoms = Phase['Atoms']
3073    cx,ct,cs,cia = General['AtomPtrs']
3074    AtomVar = RMCPdict['AtomVar']
3075    varnames = []
3076    for iat,atom in enumerate(RMCPdict['AtomConstr']):
3077        for it,item in enumerate(atom):
3078            if it > 1 and item:
3079                itnum = item.split('@')[1]
3080                varname = '@%s'%itnum
3081                varnames.append(varname)
3082                if it < 6:
3083                    if varname not in AtomVar:
3084                        AtomVar[varname] = 0.0      #put ISODISTORT mode displ here?
3085                else:
3086                    for i in range(3):
3087                        if varname not in AtomVar:
3088                            AtomVar[varname] = Atoms[iat][cia+i+2]
3089    varnames = set(varnames)
3090    for name in list(AtomVar.keys()):       #clear out unused parameters
3091        if name not in varnames:
3092            del AtomVar[name]
3093   
3094def MakePDFfitAtomsFile(Phase,RMCPdict):
3095    '''Make the PDFfit atoms file
3096    '''
3097    General = Phase['General']
3098    fName = General['Name']+'-PDFfit.stru'
3099    fName = fName.replace(' ','_')
3100    fatm = open(fName,'w')
3101    fatm.write('title  structure of '+General['Name']+'\n')
3102    fatm.write('format pdffit\n')
3103    fatm.write('scale   1.000000\n')    #fixed
3104    # if RMCPdict['shape'] == 'sphere':
3105    #     sharp = '%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0])
3106    # else:
3107    #     sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3108    sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3109    fatm.write('sharp '+sharp)
3110    shape = ''
3111    if RMCPdict['shape'] == 'sphere' and RMCPdict['spdiameter'][0] > 0.:
3112        shape = '   sphere, %10.6f\n'%RMCPdict['spdiameter'][0]
3113    elif RMCPdict['stepcut'] > 0.:
3114        shape = 'stepcut, %10.6f\n'%RMCPdict['stepcut']
3115    if shape:
3116        fatm.write('shape  '+shape)
3117    fatm.write('spcgr   %s\n'%General['SGData']['SpGrp'].replace(' ',''))
3118    cell = General['Cell'][1:7]
3119    fatm.write('cell  %10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n'%(
3120        cell[0],cell[1],cell[2],cell[3],cell[4],cell[5]))
3121    fatm.write('dcell '+5*'  0.000000,'+'  0.000000\n') 
3122    Atoms = Phase['Atoms']
3123    fatm.write('ncell %8d,%8d,%8d,%10d\n'%(1,1,1,len(Atoms)))
3124    fatm.write('atoms\n')
3125    cx,ct,cs,cia = General['AtomPtrs']
3126    for atom in Atoms:
3127        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]))
3128        fatm.write('    '+'%18.8f%18.8f%18.8f%13.4f\n'%(0.,0.,0.,0.))
3129        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(atom[cia+2],atom[cia+3],atom[cia+4]))
3130        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.,))
3131        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(atom[cia+5],atom[cia+6],atom[cia+7]))
3132        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.))
3133    fatm.close()
3134   
3135def MakePDFfitRunFile(Phase,RMCPdict):
3136    '''Make the PDFfit python run file
3137    '''
3138   
3139    def GetCellConstr(SGData):
3140        if SGData['SGLaue'] in ['m3', 'm3m']:
3141            return [1,1,1,0,0,0]
3142        elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3143            return [1,1,2,0,0,0]
3144        elif SGData['SGLaue'] in ['3R','3mR']:
3145            return [1,1,1,2,2,2]
3146        elif SGData['SGLaue'] == 'mmm':
3147            return [1,2,3,0,0,0]
3148        elif SGData['SGLaue'] == '2/m':
3149            if SGData['SGUniq'] == 'a':
3150                return [1,2,3,4,0,0]
3151            elif SGData['SGUniq'] == 'b':
3152                return [1,2,3,0,4,0]
3153            elif SGData['SGUniq'] == 'c':
3154                return [1,2,3,0,0,4]
3155        else:
3156            return [1,2,3,4,5,6]
3157       
3158    General = Phase['General']
3159    Cell = General['Cell'][1:7]
3160    G2path = GSASIIpath.path2GSAS2
3161    rundata = '''
3162#!/usr/bin/env python
3163# -*- coding: utf-8 -*-
3164import sys
3165sys.path.append('%s')
3166from diffpy.pdffit2 import PdfFit
3167pf = PdfFit()
3168'''%G2path
3169    Nd = 0
3170    Np = 0
3171    parms = {}
3172    if 'sequential' in RMCPdict['refinement']:
3173        Np = 3
3174        rundata += '#sequential data here\n'
3175    else:
3176        for file in RMCPdict['files']:
3177            if 'Select' in RMCPdict['files'][file][0]:
3178                continue
3179            if 'Neutron' in file:
3180                Nd += 1
3181                dType = 'Ndata'
3182            else:
3183                Nd += 1
3184                dType = 'Xdata'
3185            rundata += "pf.read_data('%s', '%s', 30.0, %.4f)\n"%(RMCPdict['files'][file][0],dType[0],RMCPdict[dType]['qdamp'][0])
3186            rundata += 'pf.setdata(%d)\n'%Nd
3187            rundata += 'pf.pdfrange(%d, %6.2f, %6.2f)\n'%(Nd,RMCPdict[dType]['Fitrange'][0],RMCPdict[dType]['Fitrange'][1])
3188            for item in ['dscale','qdamp','qbroad']:
3189                if RMCPdict[dType][item][1]:
3190                    Np += 1
3191                    rundata += 'pf.constrain(pf.%s(),"@%d")\n'%(item,Np)
3192                    parms[Np] = [RMCPdict[dType][item][0],item]
3193    fName = General['Name']+'-PDFfit.stru'
3194    fName = fName.replace(' ','_')
3195    Np = 9
3196    rundata += "pf.read_struct('%s')\n"%(fName)
3197    for item in ['delta1','delta2','sratio']:
3198        if RMCPdict[item][1]:
3199            Np += 1
3200            rundata += 'pf.constrain(pf.%s,"@%d")\n'%(item,Np)
3201            parms[Np] = [RMCPdict[item][0],item]
3202    if 'sphere' in RMCPdict['shape'][0] and RMCPdict['spdiameter'][1]:
3203        Np += 1
3204        rundata += 'pf.constrain(pf.spdiameter,"@%d")\n'%Np
3205        parms[Np] = [RMCPdict['spdiameter'][0],'spdiameter']
3206   
3207    if RMCPdict['cellref']:
3208        cellconst = GetCellConstr(RMCPdict['SGData'])
3209        used = []
3210        cellNames = ['a','b','c','alpha','beta','gamma']
3211        for ic in range(6):
3212            if cellconst[ic]:
3213                rundata += 'pf.constrain(pf.lat(%d), "@%d")\n'%(ic+1,Np+cellconst[ic])
3214                if cellconst[ic] not in used:
3215                    parms[Np+cellconst[ic]] = [Cell[ic],cellNames[ic]]
3216                used.append(cellconst[ic])
3217#Atom constraints here -------------------------------------------------------
3218    AtomVar = RMCPdict['AtomVar']
3219    used = []
3220    for iat,atom in enumerate(RMCPdict['AtomConstr']):
3221        for it,item in enumerate(atom):
3222            names = ['pf.x(%d)'%(iat+1),'pf.y(%d)'%(iat+1),'pf.z(%d)'%(iat+1),'pf.occ(%d)'%(iat+1)]
3223            if it > 1 and item:
3224                itnum = item.split('@')[1]
3225                if it < 6:
3226                    rundata += 'pf.constrain(%s,"%s")\n'%(names[it-2],item)
3227                    if itnum not in used:
3228                        parms[itnum] = [AtomVar['@%s'%itnum],names[it-2].split('.')[1]]
3229                        used.append(itnum)
3230                else:
3231                    uijs = ['pf.u11(%d)'%(iat+1),'pf.u22(%d)'%(iat+1),'pf.u33(%d)'%(iat+1)]     
3232                    for i in range(3):
3233                        rundata += 'pf.constrain(%s,"%s")\n'%(uijs[i],item)
3234                        if itnum not in used:
3235                            parms[itnum] = [AtomVar['@%s'%itnum],uijs[i].split('.')[1]]
3236                            used.append(itnum)
3237                           
3238    if 'sequential' in RMCPdict['refinement']:
3239        rundata += '#parameters here\n'
3240    else:       
3241# set parameter values
3242        RMCPdict['Parms'] = parms
3243        for iprm in parms:
3244            rundata += 'pf.setpar(%s,%.6f)\n'%(iprm,parms[iprm][0])
3245                       
3246# Refine & Save results ---------------------------------------------------------------   
3247    rundata += 'pf.refine()\n'
3248    if 'sequential' in RMCPdict['refinement']:
3249        fName = 'Sequential_PDFfit'
3250    else:
3251        fName = General['Name'].replace(' ','_')+'-PDFfit'
3252    Nd = 0   
3253    for file in RMCPdict['files']:
3254        if 'Select' in RMCPdict['files'][file][0]:
3255            continue
3256        Nd += 1
3257        rundata += 'pf.save_pdf(%d, "%s")\n'%(Nd,fName+file[0]+'.fgr')
3258       
3259    rundata += 'pf.save_struct(1, "%s")\n'%(fName+'.rstr')
3260    rundata += 'pf.save_res("%s")\n'%(fName+'.res')
3261 
3262    if 'sequential' in RMCPdict['refinement']:
3263        rfile = open('Seq_PDFfit_template.py','w')
3264    else:
3265        rfile = open(fName+'.py','w')
3266    rfile.writelines(rundata)
3267    rfile.close()
3268   
3269    return fName+'.py'
3270   
3271def UpdatePDFfit(Phase,RMCPdict):
3272    ''' Updates various PDFfit parameters held in GSAS-II
3273    '''
3274   
3275    General = Phase['General']
3276    if RMCPdict['refinement'] == 'normal':
3277        fName = General['Name']+'-PDFfit.rstr'
3278        try:
3279            rstr = open(fName.replace(' ','_'),'r')
3280        except FileNotFoundError:
3281            return [fName,'Not found - PDFfit failed']
3282        lines = rstr.readlines()
3283        rstr.close()
3284        header = [line[:-1].split(' ',1) for line in lines[:7]]
3285        resdict = dict(header)
3286        for item in ['scale','sharp','cell']:
3287            resdict[item] = [float(val) for val in resdict[item].split(',')]
3288        General['Cell'][1:7] = resdict['cell']
3289        for inam,name in enumerate(['delta2','delta1','sratio']):
3290            RMCPdict[name][0] = resdict['sharp'][inam]
3291        if 'shape' in resdict:
3292            if 'sphere' in resdict['shape']:
3293                RMCPdict['spdiameter'][0] = resdict['shape'][-1]
3294            else:
3295                RMCPdict['stepcut'][0] = resdict['shape'][-1]
3296        cx,ct,cs,ci = G2mth.getAtomPtrs(Phase)     
3297        Atoms = Phase['Atoms']
3298        atmBeg = 0
3299        for line in lines:
3300            atmBeg += 1
3301            if 'atoms' in line:
3302                break
3303        for atom in Atoms:
3304            atstr = lines[atmBeg][:-1].split()
3305            Uiistr = lines[atmBeg+2][:-1].split()
3306            Uijstr = lines[atmBeg+4][:-1].split()
3307            atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])]
3308            atom[ci] = 'A'
3309            atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])]
3310            atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])]
3311            atmBeg += 6
3312           
3313    fName = 'Sequential_PDFfit.res'
3314    if RMCPdict['refinement'] == 'normal':
3315        fName = General['Name']+'-PDFfit.res'
3316    try:
3317        res = open(fName.replace(' ','_'),'r')
3318    except FileNotFoundError:
3319        return [fName,'Not found - PDFfit failed']
3320    lines = res.readlines()
3321    res.close()
3322    Ibeg = False
3323    resline = ''
3324    for iline,line in enumerate(lines):
3325        if 'Refinement parameters' in line:
3326            Ibeg = True
3327            continue
3328        if Ibeg:
3329            if '---------' in line:
3330                break
3331            resline += line[:-1]
3332    for iline,line in enumerate(lines):
3333        if 'Rw - ' in line:
3334            Rwp = float(line.split(':')[1])
3335    results = resline.replace('(','').split(')')[:-1]
3336    results = ['@'+result.lstrip() for result in results]
3337    results = [item.split() for item in results]
3338    if RMCPdict['refinement'] == 'normal':
3339        results = dict([[item[0][:-1],float(item[1])] for item in results if item[0][:-1] in RMCPdict['AtomVar']])
3340        RMCPdict['AtomVar'].update(results)
3341        return None
3342    else:   #sequential
3343        newParms = dict([[item[0][1:-1],[float(item[1]),float(item[2])]] for item in results])
3344        return newParms,Rwp
3345       
3346def MakefullrmcRun(pName,Phase,RMCPdict):
3347    '''Creates a script to run fullrmc. Returns the name of the file that was
3348    created.
3349    '''
3350    BondList = {}
3351    for k in RMCPdict['Pairs']:
3352        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
3353            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
3354    AngleList = []
3355    for angle in RMCPdict['Angles']:
3356        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
3357            continue
3358        for i in (0,1,2):
3359            angle[i] = angle[i].strip()
3360        AngleList.append(angle)
3361    # rmin = RMCPdict['min Contact']
3362    cell = Phase['General']['Cell'][1:7]
3363    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
3364    cx,ct,cs,cia = Phase['General']['AtomPtrs']
3365    atomsList = []
3366    for atom in Phase['Atoms']:
3367        el = ''.join([i for i in atom[ct] if i.isalpha()])
3368        atomsList.append([el] + atom[cx:cx+4])
3369    projDir,projName = os.path.split(os.path.abspath(pName))
3370    scrname = pName+'-fullrmc.py'
3371    restart = '%s_restart.pdb'%pName
3372    Files = RMCPdict['files']
3373    rundata = ''
3374    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
3375    rundata += '# created in '+__file__+" v"+filversion.split()[1]
3376    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
3377    rundata += '''
3378# fullrmc imports (all that are potentially useful)
3379import os,glob
3380import time
3381import pickle
3382import numpy as np
3383from fullrmc.Core import Collection
3384from fullrmc.Engine import Engine
3385import fullrmc.Constraints.PairDistributionConstraints as fPDF
3386from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
3387from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
3388from fullrmc.Constraints.BondConstraints import BondConstraint
3389from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
3390from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
3391from fullrmc.Generators.Swaps import SwapPositionsGenerator
3392# utility routines
3393def writeHeader(ENGINE,statFP):
3394    'header for stats file'
3395    statFP.write('generated-steps, total-error, ')
3396    for c in ENGINE.constraints:
3397        statFP.write(c.constraintName)
3398        statFP.write(', ')
3399    statFP.write('\\n')
3400    statFP.flush()
3401   
3402def writeCurrentStatus(ENGINE,statFP,plotF):
3403    'line in stats file & current constraint plots'
3404    statFP.write(str(ENGINE.generated))
3405    statFP.write(', ')
3406    statFP.write(str(ENGINE.totalStandardError))
3407    statFP.write(', ')
3408    for c in ENGINE.constraints:
3409        statFP.write(str(c.standardError))
3410        statFP.write(', ')
3411    statFP.write('\\n')
3412    statFP.flush()
3413    mpl.use('agg')
3414    fp = open(plotF,'wb')
3415    for c in ENGINE.constraints:
3416        p = c.plot(show=False)
3417        p[0].canvas.draw()
3418        image = p[0].canvas.buffer_rgba()
3419        pickle.dump(c.constraintName,fp)
3420        pickle.dump(np.array(image),fp)
3421    fp.close()
3422
3423def calcRmax(ENGINE):
3424    'from Bachir, works for non-othorhombic'
3425    a,b,c = ENGINE.basisVectors
3426    lens = []
3427    ts    = np.linalg.norm(np.cross(a,b))/2
3428    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3429    ts = np.linalg.norm(np.cross(b,c))/2
3430    lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3431    ts = np.linalg.norm(np.cross(a,c))/2
3432    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3433    return min(lens)
3434'''
3435    rundata += '''
3436### When True, erases an existing engine to provide a fresh start
3437FRESH_START = {:}
3438dirName = "{:}"
3439prefix = "{:}"
3440project = prefix + "-fullrmc"
3441time0 = time.time()
3442'''.format(RMCPdict['ReStart'][0],projDir,projName)
3443   
3444    rundata += '# setup structure\n'
3445    rundata += 'cell = ' + str(cell) + '\n'
3446    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3447    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
3448    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3449
3450    rundata += '\n# initialize engine\n'
3451    rundata += '''
3452engineFileName = os.path.join(dirName, project + '.rmc')
3453projectStats = os.path.join(dirName, project + '.stats')
3454projectPlots = os.path.join(dirName, project + '.plots')
3455pdbFile = os.path.join(dirName, project + '_restart.pdb')
3456# check Engine exists if so (and not FRESH_START) load it
3457# otherwise build it
3458ENGINE = Engine(path=None)
3459if not ENGINE.is_engine(engineFileName) or FRESH_START:
3460    ## create structure
3461    ENGINE = Engine(path=engineFileName, freshStart=True)
3462    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
3463                                 atoms      = atomList,
3464                                 unitcellBC = cell,
3465                                 supercell  = supercell)
3466'''   
3467    import atmdata
3468    # rundata += '    # conversion factors (may be needed)\n'
3469    # rundata += '    sumCiBi2 = 0.\n'
3470    # for elem in Phase['General']['AtomTypes']:
3471    #     rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3472    #     rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3473    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3474    # settings that require a new Engine
3475    for File in Files:
3476        filDat = RMCPdict['files'][File]
3477        if not os.path.exists(filDat[0]): continue
3478        sfwt = 'neutronCohb'
3479        if 'Xray' in File:
3480            sfwt = 'atomicNumber'
3481        if 'G(r)' in File:
3482            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3483            if filDat[3] == 0:
3484                #rundata += '''    # read and xform G(r) as defined in RMCProfile
3485    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3486                #rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
3487                #rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3488                rundata += '    # G(r) as defined in RMCProfile\n'
3489                rundata += '    GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3490            elif filDat[3] == 1:
3491                rundata += '    # This is G(r) as defined in PDFFIT\n'
3492                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3493            elif filDat[3] == 2:
3494                rundata += '    # This is g(r)\n'
3495                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3496            else:
3497                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3498            rundata += '    ENGINE.add_constraints([GofR])\n'
3499            rundata += '    GofR.set_limits((None, calcRmax(ENGINE)))\n'
3500        elif '(Q)' in File:
3501            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3502            if filDat[3] == 0:
3503                rundata += '    # F(Q) as defined in RMCProfile\n'
3504                #rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
3505                if filDat[4]:
3506                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3507                rundata += '    SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3508            elif filDat[3] == 1:
3509                rundata += '    # S(Q) as defined in PDFFIT\n'
3510                rundata += '    SOQ[1] -= 1\n'
3511                if filDat[4]:
3512                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3513                rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3514            else:
3515                raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3516            rundata += '    ENGINE.add_constraints([SofQ])\n'
3517        else:
3518            print('What is this?')
3519    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
3520    if BondList:
3521        rundata += '''    B_CONSTRAINT   = BondConstraint()
3522    ENGINE.add_constraints(B_CONSTRAINT)
3523    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3524'''
3525        for pair in BondList:
3526            e1,e2 = pair.split('-')
3527            d1,d2 = BondList[pair]
3528            if d1 == 0: continue
3529            if d2 == 0:
3530                print('Ignoring min distance without maximum')
3531                #rundata += '            ("element","{}","{}",{}),\n'.format(
3532                #                        e1.strip(),e2.strip(),d1)
3533            else:
3534                rundata += '            ("element","{}","{}",{},{}),\n'.format(
3535                                        e1.strip(),e2.strip(),d1,d2)
3536        rundata += '             ])\n'
3537    if AngleList:
3538        rundata += '''    A_CONSTRAINT   = BondsAngleConstraint()
3539    ENGINE.add_constraints(A_CONSTRAINT)
3540    A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3541'''
3542        for item in AngleList:
3543            rundata += ('            '+
3544               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3545        rundata += '             ])\n'
3546    rundata += '''
3547    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3548    ENGINE.save()
3549else:
3550    ENGINE = ENGINE.load(path=engineFileName)
3551
3552ENGINE.set_log_file(os.path.join(dirName,prefix))
3553'''
3554    if RMCPdict['Swaps']:
3555        rundata += '\n#set up for site swaps\n'
3556        rundata += 'aN = ENGINE.allNames\n'
3557        rundata += 'SwapGen = {}\n'
3558        for swap in RMCPdict['Swaps']:
3559            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3560            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3561            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3562        rundata += '    for swaps in SwapGen:\n'
3563        rundata += '        AB = swaps.split("-")\n'
3564        rundata += '        ENGINE.set_groups_as_atoms()\n'
3565        rundata += '        for g in ENGINE.groups:\n'
3566        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
3567        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
3568        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
3569        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
3570        rundata += '            sProb = SwapGen[swaps][2]\n'
3571    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3572    # rundata += 'wtDict = {}\n'
3573    # for File in Files:
3574    #     filDat = RMCPdict['files'][File]
3575    #     if not os.path.exists(filDat[0]): continue
3576    #     if 'Xray' in File:
3577    #         sfwt = 'atomicNumber'
3578    #     else:
3579    #         sfwt = 'neutronCohb'
3580    #     if 'G(r)' in File:
3581    #         typ = 'Pair'
3582    #     elif '(Q)' in File:
3583    #         typ = 'Struct'
3584    #     rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3585    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
3586    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
3587    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3588    rundata += '        c.set_limits((None,calcRmax(ENGINE)))\n'
3589    if RMCPdict['FitScale']:
3590        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3591    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3592    if RMCPdict['FitScale']:
3593        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
3594        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3595    # torsions difficult to implement, must be internal to cell & named with
3596    # fullrmc atom names
3597    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
3598    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
3599    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
3600    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3601    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3602    #     for torsion in RMCPdict['Torsions']:
3603    #         rundata += '    %s\n'%str(tuple(torsion))
3604    #     rundata += '        ]})\n'           
3605    rundata += '''
3606if FRESH_START:
3607    # initialize engine with one step to get starting config energetics
3608    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3609    statFP = open(projectStats,'w')
3610    writeHeader(ENGINE,statFP)
3611    writeCurrentStatus(ENGINE,statFP,projectPlots)
3612else:
3613    statFP = open(projectStats,'a')
3614
3615# setup runs for fullrmc
3616'''
3617    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3618    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3619    rundata += '    ENGINE.set_groups_as_atoms()\n'
3620    rundata += '    expected = ENGINE.generated+steps\n'
3621   
3622    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3623    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3624    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
3625    rundata += 'statFP.close()\n'
3626    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3627    rfile = open(scrname,'w')
3628    rfile.writelines(rundata)
3629    rfile.close()
3630    return scrname
3631   
3632def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3633    bondDist = []
3634    Cell = general['Cell'][1:7]
3635    Supercell =  RMCPdict['SuperCell']
3636    Trans = np.eye(3)*np.array(Supercell)
3637    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3638    Amat,Bmat = G2lat.cell2AB(Cell)
3639    indices = (-1,0,1)
3640    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3641    for bonds in bondList:
3642        Oxyz = np.array(Atoms[bonds[0]][1:])
3643        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
3644        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3645        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3646        for dx in Dx.T:
3647            bondDist.append(np.min(dx))
3648    return np.array(bondDist)
3649   
3650def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3651    bondAngles = []
3652    Cell = general['Cell'][1:7]
3653    Supercell =  RMCPdict['SuperCell']
3654    Trans = np.eye(3)*np.array(Supercell)
3655    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3656    Amat,Bmat = G2lat.cell2AB(Cell)
3657    indices = (-1,0,1)
3658    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3659    for angle in angleList:
3660        Oxyz = np.array(Atoms[angle[0]][1:])
3661        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
3662        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])
3663        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3664        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3665        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3666        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3667        iDAx = np.argmin(DAx,axis=0)
3668        iDBx = np.argmin(DBx,axis=0)
3669        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3670            DAv = DAxV[iA,i]/DAx[iA,i]
3671            DBv = DBxV[iB,i]/DBx[iB,i]
3672            bondAngles.append(npacosd(np.sum(DAv*DBv)))
3673    return np.array(bondAngles)
3674   
3675#### Reflectometry calculations ################################################################################
3676def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3677    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3678   
3679    class RandomDisplacementBounds(object):
3680        """random displacement with bounds"""
3681        def __init__(self, xmin, xmax, stepsize=0.5):
3682            self.xmin = xmin
3683            self.xmax = xmax
3684            self.stepsize = stepsize
3685   
3686        def __call__(self, x):
3687            """take a random step but ensure the new position is within the bounds"""
3688            while True:
3689                # this could be done in a much more clever way, but it will work for example purposes
3690                steps = self.xmax-self.xmin
3691                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
3692                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3693                    break
3694            return xnew
3695   
3696    def GetModelParms():
3697        parmDict = {}
3698        varyList = []
3699        values = []
3700        bounds = []
3701        parmDict['dQ type'] = data['dQ type']
3702        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
3703        for parm in ['Scale','FltBack']:
3704            parmDict[parm] = data[parm][0]
3705            if data[parm][1]:
3706                varyList.append(parm)
3707                values.append(data[parm][0])
3708                bounds.append(Bounds[parm])
3709        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3710        parmDict['nLayers'] = len(parmDict['Layer Seq'])
3711        for ilay,layer in enumerate(data['Layers']):
3712            name = layer['Name']
3713            cid = str(ilay)+';'
3714            parmDict[cid+'Name'] = name
3715            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3716                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3717                if layer.get(parm,[0.,False])[1]:
3718                    varyList.append(cid+parm)
3719                    value = layer[parm][0]
3720                    values.append(value)
3721                    if value:
3722                        bound = [value*Bfac,value/Bfac]
3723                    else:
3724                        bound = [0.,10.]
3725                    bounds.append(bound)
3726            if name not in ['vacuum','unit scatter']:
3727                parmDict[cid+'rho'] = Substances[name]['Scatt density']
3728                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3729        return parmDict,varyList,values,bounds
3730   
3731    def SetModelParms():
3732        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3733        if 'Scale' in varyList:
3734            data['Scale'][0] = parmDict['Scale']
3735            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
3736        G2fil.G2Print (line)
3737        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3738        if 'FltBack' in varyList:
3739            data['FltBack'][0] = parmDict['FltBack']
3740            line += ' esd: %15.3g'%(sigDict['FltBack'])
3741        G2fil.G2Print (line)
3742        for ilay,layer in enumerate(data['Layers']):
3743            name = layer['Name']
3744            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3745            cid = str(ilay)+';'
3746            line = ' '
3747            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3748            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3749            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3750                if parm in layer:
3751                    if parm == 'Rough':
3752                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
3753                    else:
3754                        layer[parm][0] = parmDict[cid+parm]
3755                    line += ' %s: %.3f'%(parm,layer[parm][0])
3756                    if cid+parm in varyList:
3757                        line += ' esd: %.3g'%(sigDict[cid+parm])
3758            G2fil.G2Print (line)
3759            G2fil.G2Print (line2)
3760   
3761    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3762        parmDict.update(zip(varyList,values))
3763        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3764        return M
3765   
3766    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3767        parmDict.update(zip(varyList,values))
3768        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3769        return np.sum(M**2)
3770   
3771    def getREFD(Q,Qsig,parmDict):
3772        Ic = np.ones_like(Q)*parmDict['FltBack']
3773        Scale = parmDict['Scale']
3774        Nlayers = parmDict['nLayers']
3775        Res = parmDict['Res']
3776        depth = np.zeros(Nlayers)
3777        rho = np.zeros(Nlayers)
3778        irho = np.zeros(Nlayers)
3779        sigma = np.zeros(Nlayers)
3780        for ilay,lay in enumerate(parmDict['Layer Seq']):
3781            cid = str(lay)+';'
3782            depth[ilay] = parmDict[cid+'Thick']
3783            sigma[ilay] = parmDict[cid+'Rough']
3784            if parmDict[cid+'Name'] == u'unit scatter':
3785                rho[ilay] = parmDict[cid+'DenMul']
3786                irho[ilay] = parmDict[cid+'iDenMul']
3787            elif 'vacuum' != parmDict[cid+'Name']:
3788                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3789                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3790            if cid+'Mag SLD' in parmDict:
3791                rho[ilay] += parmDict[cid+'Mag SLD']
3792        if parmDict['dQ type'] == 'None':
3793            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3794        elif 'const' in parmDict['dQ type']:
3795            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
3796        else:       #dQ/Q in data
3797            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3798        Ic += AB*Scale
3799        return Ic
3800       
3801    def estimateT0(takestep):
3802        Mmax = -1.e-10
3803        Mmin = 1.e10
3804        for i in range(100):
3805            x0 = takestep(values)
3806            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3807            Mmin = min(M,Mmin)
3808            MMax = max(M,Mmax)
3809        return 1.5*(MMax-Mmin)
3810
3811    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3812    if data.get('2% weight'):
3813        wt = 1./(0.02*Io)**2
3814    Qmin = Limits[1][0]
3815    Qmax = Limits[1][1]
3816    wtFactor = ProfDict['wtFactor']
3817    Bfac = data['Toler']
3818    Ibeg = np.searchsorted(Q,Qmin)
3819    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
3820    Ic[:] = 0
3821    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3822              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3823    parmDict,varyList,values,bounds = GetModelParms()
3824    Msg = 'Failed to converge'
3825    if varyList:
3826        if data['Minimizer'] == 'LMLS': 
3827            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3828                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3829            parmDict.update(zip(varyList,result[0]))
3830            chisq = np.sum(result[2]['fvec']**2)
3831            ncalc = result[2]['nfev']
3832            covM = result[1]
3833            newVals = result[0]
3834        elif data['Minimizer'] == 'Basin Hopping':
3835            xyrng = np.array(bounds).T
3836            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3837            T0 = estimateT0(take_step)
3838            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3839            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3840                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3841                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3842            chisq = result.fun
3843            ncalc = result.nfev
3844            newVals = result.x
3845            covM = []
3846        elif data['Minimizer'] == 'MC/SA Anneal':
3847            xyrng = np.array(bounds).T
3848            result = G2mth.anneal(sumREFD, values, 
3849                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3850                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3851                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3852                ranRange=0.20,autoRan=False,dlg=None)
3853            newVals = result[0]
3854            parmDict.update(zip(varyList,newVals))
3855            chisq = result[1]
3856            ncalc = result[3]
3857            covM = []
3858            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3859        elif data['Minimizer'] == 'L-BFGS-B':
3860            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
3861                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3862            parmDict.update(zip(varyList,result['x']))
3863            chisq = result.fun
3864            ncalc = result.nfev
3865            newVals = result.x
3866            covM = []
3867    else:   #nothing varied
3868        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3869        chisq = np.sum(M**2)
3870        ncalc = 0
3871        covM = []
3872        sig = []
3873        sigDict = {}
3874        result = []
3875    Rvals = {}
3876    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
3877    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
3878    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3879    Ib[Ibeg:Ifin] = parmDict['FltBack']
3880    try:
3881        if not len(varyList):
3882            Msg += ' - nothing refined'
3883            raise ValueError
3884        Nans = np.isnan(newVals)
3885        if np.any(Nans):
3886            name = varyList[Nans.nonzero(True)[0]]
3887            Msg += ' Nan result for '+name+'!'
3888            raise ValueError
3889        Negs = np.less_equal(newVals,0.)
3890        if np.any(Negs):
3891            indx = Negs.nonzero()
3892            name = varyList[indx[0][0]]
3893            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3894                Msg += ' negative coefficient for '+name+'!'
3895                raise ValueError
3896        if len(covM):
3897            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3898            covMatrix = covM*Rvals['GOF']
3899        else:
3900            sig = np.zeros(len(varyList))
3901            covMatrix = []
3902        sigDict = dict(zip(varyList,sig))
3903        G2fil.G2Print (' Results of reflectometry data modelling fit:')
3904        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3905        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3906        SetModelParms()
3907        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3908    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
3909        G2fil.G2Print (Msg)
3910        return False,0,0,0,0,0,0,Msg
3911       
3912def makeSLDprofile(data,Substances):
3913   
3914    sq2 = np.sqrt(2.)
3915    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3916    Nlayers = len(laySeq)
3917    laySeq = np.array(laySeq,dtype=int)
3918    interfaces = np.zeros(Nlayers)
3919    rho = np.zeros(Nlayers)
3920    sigma = np.zeros(Nlayers)
3921    name = data['Layers'][0]['Name']
3922    thick = 0.
3923    for ilay,lay in enumerate(laySeq):
3924        layer = data['Layers'][lay]
3925        name = layer['Name']
3926        if 'Thick' in layer:
3927            thick += layer['Thick'][0]
3928            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3929        if 'Rough' in layer:
3930            sigma[ilay] = max(0.001,layer['Rough'][0])
3931        if name != 'vacuum':
3932            if name == 'unit scatter':
3933                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
3934            else:
3935                rrho = Substances[name]['Scatt density']
3936                irho = Substances[name]['XImag density']
3937                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
3938        if 'Mag SLD' in layer:
3939            rho[ilay] += layer['Mag SLD'][0]
3940    name = data['Layers'][-1]['Name']
3941    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
3942    xr = np.flipud(x)
3943    interfaces[-1] = x[-1]
3944    y = np.ones_like(x)*rho[0]
3945    iBeg = 0
3946    for ilayer in range(Nlayers-1):
3947        delt = rho[ilayer+1]-rho[ilayer]
3948        iPos = np.searchsorted(x,interfaces[ilayer])
3949        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
3950        iBeg = iPos
3951    return x,xr,y   
3952
3953def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3954   
3955    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3956    Qmin = Limits[1][0]
3957    Qmax = Limits[1][1]
3958    iBeg = np.searchsorted(Q,Qmin)
3959    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3960    Ib[:] = data['FltBack'][0]
3961    Ic[:] = 0
3962    Scale = data['Scale'][0]
3963    if data['Layer Seq'] == []:
3964        return
3965    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3966    Nlayers = len(laySeq)
3967    depth = np.zeros(Nlayers)
3968    rho = np.zeros(Nlayers)
3969    irho = np.zeros(Nlayers)
3970    sigma = np.zeros(Nlayers)
3971    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3972        layer = data['Layers'][lay]
3973        name = layer['Name']
3974        if 'Thick' in layer:    #skips first & last layers
3975            depth[ilay] = layer['Thick'][0]
3976        if 'Rough' in layer:    #skips first layer
3977            sigma[ilay] = layer['Rough'][0]
3978        if 'unit scatter' == name:
3979            rho[ilay] = layer['DenMul'][0]
3980            irho[ilay] = layer['iDenMul'][0]
3981        else:
3982            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
3983            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
3984        if 'Mag SLD' in layer:
3985            rho[ilay] += layer['Mag SLD'][0]
3986    if data['dQ type'] == 'None':
3987        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3988    elif 'const' in data['dQ type']:
3989        res = data['Resolution'][0]/(100.*sateln2)
3990        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
3991    else:       #dQ/Q in data
3992        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
3993    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
3994
3995def abeles(kz, depth, rho, irho=0, sigma=0):
3996    """
3997    Optical matrix form of the reflectivity calculation.
3998    O.S. Heavens, Optical Properties of Thin Solid Films
3999   
4000    Reflectometry as a function of kz for a set of slabs.
4001
4002    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
4003        This is :math:`\\tfrac12 Q_z`.       
4004    :param depth:  float[m] (Ang).
4005        thickness of each layer.  The thickness of the incident medium
4006        and substrate are ignored.
4007    :param rho:  float[n,k] (1e-6/Ang^2)
4008        Real scattering length density for each layer for each kz
4009    :param irho:  float[n,k] (1e-6/Ang^2)
4010        Imaginary scattering length density for each layer for each kz
4011        Note: absorption cross section mu = 2 irho/lambda for neutrons
4012    :param sigma: float[m-1] (Ang)
4013        interfacial roughness.  This is the roughness between a layer
4014        and the previous layer. The sigma array should have m-1 entries.
4015
4016    Slabs are ordered with the surface SLD at index 0 and substrate at
4017    index -1, or reversed if kz < 0.
4018    """
4019    def calc(kz, depth, rho, irho, sigma):
4020        if len(kz) == 0: return kz
4021   
4022        # Complex index of refraction is relative to the incident medium.
4023        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
4024        # in place of kz^2, and ignoring rho_o
4025        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
4026        k = kz
4027   
4028        # According to Heavens, the initial matrix should be [ 1 F; F 1],
4029        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
4030        # multiply versus some coding convenience.
4031        B11 = 1
4032        B22 = 1
4033        B21 = 0
4034        B12 = 0
4035        for i in range(0, len(depth)-1):
4036            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
4037            F = (k - k_next) / (k + k_next)
4038            F *= np.exp(-2*k*k_next*sigma[i]**2)
4039            #print "==== layer",i
4040            #print "kz:", kz
4041            #print "k:", k
4042            #print "k_next:",k_next
4043            #print "F:",F
4044            #print "rho:",rho[:,i+1]
4045            #print "irho:",irho[:,i+1]
4046            #print "d:",depth[i],"sigma:",sigma[i]
4047            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
4048            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
4049            M21 = F*M11
4050            M12 = F*M22
4051            C1 = B11*M11 + B21*M12
4052            C2 = B11*M21 + B21*M22
4053            B11 = C1
4054            B21 = C2
4055            C1 = B12*M11 + B22*M12
4056            C2 = B12*M21 + B22*M22
4057            B12 = C1
4058            B22 = C2
4059            k = k_next
4060   
4061        r = B12/B11
4062        return np.absolute(r)**2
4063
4064    if np.isscalar(kz): kz = np.asarray([kz], 'd')
4065
4066    m = len(depth)
4067
4068    # Make everything into arrays
4069    depth = np.asarray(depth,'d')
4070    rho = np.asarray(rho,'d')
4071    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
4072    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
4073
4074    # Repeat rho,irho columns as needed
4075    if len(rho.shape) == 1:
4076        rho = rho[None,:]
4077        irho = irho[None,:]
4078
4079    return calc(kz, depth, rho, irho, sigma)
4080   
4081def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
4082    y = abeles(kz, depth, rho, irho, sigma)
4083    s = dq/2.
4084    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
4085    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
4086    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
4087    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
4088    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
4089    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
4090    y *= 0.137023
4091    return y
4092       
4093def makeRefdFFT(Limits,Profile):
4094    G2fil.G2Print ('make fft')
4095    Q,Io = Profile[:2]
4096    Qmin = Limits[1][0]
4097    Qmax = Limits[1][1]
4098    iBeg = np.searchsorted(Q,Qmin)
4099    iFin = np.searchsorted(Q,Qmax)+1    #include last point
4100    Qf = np.linspace(0.,Q[iFin-1],5000)
4101    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
4102    If = QI(Qf)*Qf**4
4103    R = np.linspace(0.,5000.,5000)
4104    F = fft.rfft(If)
4105    return R,F
4106
4107   
4108#### Stacking fault simulation codes ################################################################################
4109def GetStackParms(Layers):
4110   
4111    Parms = []
4112#cell parms
4113    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
4114        Parms.append('cellA')
4115        Parms.append('cellC')
4116    else:
4117        Parms.append('cellA')
4118        Parms.append('cellB')
4119        Parms.append('cellC')
4120        if Layers['Laue'] != 'mmm':
4121            Parms.append('cellG')
4122#Transition parms
4123    for iY in range(len(Layers['Layers'])):
4124        for iX in range(len(Layers['Layers'])):
4125            Parms.append('TransP;%d;%d'%(iY,iX))
4126            Parms.append('TransX;%d;%d'%(iY,iX))
4127            Parms.append('TransY;%d;%d'%(iY,iX))
4128            Parms.append('TransZ;%d;%d'%(iY,iX))
4129    return Parms
4130
4131def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
4132    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
4133   
4134    :param dict Layers: dict with following items
4135
4136      ::
4137
4138       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
4139       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
4140       'Layers':[],'Stacking':[],'Transitions':[]}
4141       
4142    :param str ctrls: controls string to be written on DIFFaX controls.dif file
4143    :param float scale: scale factor
4144    :param dict background: background parameters
4145    :param list limits: min/max 2-theta to be calculated
4146    :param dict inst: instrument parameters dictionary
4147    :param list profile: powder pattern data
4148   
4149    Note that parameters all updated in place   
4150    '''
4151    import atmdata
4152    path = sys.path
4153    for name in path:
4154        if 'bin' in name:
4155            DIFFaX = name+'/DIFFaX.exe'
4156            G2fil.G2Print (' Execute '+DIFFaX)
4157            break
4158    # make form factor file that DIFFaX wants - atom types are GSASII style
4159    sf = open('data.sfc','w')
4160    sf.write('GSASII special form factor file for DIFFaX\n\n')
4161    atTypes = list(Layers['AtInfo'].keys())
4162    if 'H' not in atTypes:
4163        atTypes.insert(0,'H')
4164    for atType in atTypes:
4165        if atType == 'H': 
4166            blen = -.3741
4167        else:
4168            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
4169        Adat = atmdata.XrayFF[atType]
4170        text = '%4s'%(atType.ljust(4))
4171        for i in range(4):
4172            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
4173        text += '%11.6f%11.6f'%(Adat['fc'],blen)
4174        text += '%3d\n'%(Adat['Z'])
4175        sf.write(text)
4176    sf.close()
4177    #make DIFFaX control.dif file - future use GUI to set some of these flags
4178    cf = open('control.dif','w')
4179    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4180        x0 = profile[0]
4181        iBeg = np.searchsorted(x0,limits[0])
4182        iFin = np.searchsorted(x0,limits[1])+1
4183        if iFin-iBeg > 20000:
4184            iFin = iBeg+20000
4185        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
4186        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4187        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
4188    else:
4189        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4190        inst = {'Type':['XSC','XSC',]}
4191    cf.close()
4192    #make DIFFaX data file
4193    df = open('GSASII-DIFFaX.dat','w')
4194    df.write('INSTRUMENTAL\n')
4195    if 'X' in inst['Type'][0]:
4196        df.write('X-RAY\n')
4197    elif 'N' in inst['Type'][0]:
4198        df.write('NEUTRON\n')
4199    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4200        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
4201        U = ateln2*inst['U'][1]/10000.
4202        V = ateln2*inst['V'][1]/10000.
4203        W = ateln2*inst['W'][1]/10000.
4204        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
4205        HW = np.sqrt(np.mean(HWHM))
4206    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
4207        if 'Mean' in Layers['selInst']:
4208            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
4209        elif 'Gaussian' in Layers['selInst']:
4210            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
4211        else:
4212            df.write('None\n')
4213    else:
4214        df.write('0.10\nNone\n')
4215    df.write('STRUCTURAL\n')
4216    a,b,c = Layers['Cell'][1:4]
4217    gam = Layers['Cell'][6]
4218    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
4219    laue = Layers['Laue']
4220    if laue == '2/m(ab)':
4221        laue = '2/m(1)'
4222    elif laue == '2/m(c)':
4223        laue = '2/m(2)'
4224    if 'unknown' in Layers['Laue']:
4225        df.write('%s %.3f\n'%(laue,Layers['Toler']))
4226    else:   
4227        df.write('%s\n'%(laue))
4228    df.write('%d\n'%(len(Layers['Layers'])))
4229    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
4230        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
4231    layerNames = []
4232    for layer in Layers['Layers']:
4233        layerNames.append(layer['Name'])
4234    for il,layer in enumerate(Layers['Layers']):
4235        if layer['SameAs']:
4236            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
4237            continue
4238        df.write('LAYER %d\n'%(il+1))
4239        if '-1' in layer['Symm']:
4240            df.write('CENTROSYMMETRIC\n')
4241        else:
4242            df.write('NONE\n')
4243        for ia,atom in enumerate(layer['Atoms']):
4244            [name,atype,x,y,z,frac,Uiso] = atom
4245            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
4246                frac /= 2.
4247            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
4248    df.write('STACKING\n')
4249    df.write('%s\n'%(Layers['Stacking'][0]))
4250    if 'recursive' in Layers['Stacking'][0]:
4251        df.write('%s\n'%Layers['Stacking'][1])
4252    else:
4253        if 'list' in Layers['Stacking'][1]:
4254            Slen = len(Layers['Stacking'][2])
4255            iB = 0
4256            iF = 0
4257            while True:
4258                iF += 68
4259                if iF >= Slen:
4260                    break
4261                iF = min(iF,Slen)
4262                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
4263                iB = iF
4264        else:
4265            df.write('%s\n'%Layers['Stacking'][1])   
4266    df.write('TRANSITIONS\n')
4267    for iY in range(len(Layers['Layers'])):
4268        sumPx = 0.
4269        for iX in range(len(Layers['Layers'])):
4270            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
4271            p = round(p,3)
4272            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
4273            sumPx += p
4274        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
4275            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
4276            df.close()
4277            os.remove('data.sfc')
4278            os.remove('control.dif')
4279            os.remove('GSASII-DIFFaX.dat')
4280            return       
4281    df.close()   
4282    time0 = time.time()
4283    try:
4284        subp.call(DIFFaX)
4285    except OSError:
4286        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
4287    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
4288    if os.path.exists('GSASII-DIFFaX.spc'):
4289        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
4290        iFin = iBeg+Xpat.shape[1]
4291        bakType,backDict,backVary = SetBackgroundParms(background)
4292        backDict['Lam1'] = G2mth.getWave(inst)
4293        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
4294        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
4295        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
4296            rv = st.poisson(profile[3][iBeg:iFin])
4297            profile[1][iBeg:iFin] = rv.rvs()
4298            Z = np.ones_like(profile[3][iBeg:iFin])
4299            Z[1::2] *= -1
4300            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
4301            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
4302        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
4303    #cleanup files..
4304        os.remove('GSASII-DIFFaX.spc')
4305    elif os.path.exists('GSASII-DIFFaX.sadp'):
4306        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
4307        Sadp = np.reshape(Sadp,(256,-1))
4308        Layers['Sadp']['Img'] = Sadp
4309        os.remove('GSASII-DIFFaX.sadp')
4310    os.remove('data.sfc')
4311    os.remove('control.dif')
4312    os.remove('GSASII-DIFFaX.dat')
4313   
4314def SetPWDRscan(inst,limits,profile):
4315   
4316    wave = G2mth.getMeanWave(inst)
4317    x0 = profile[0]
4318    iBeg = np.searchsorted(x0,limits[0])
4319    iFin = np.searchsorted(x0,limits[1])
4320    if iFin-iBeg > 20000:
4321        iFin = iBeg+20000
4322    Dx = (x0[iFin]-x0[iBeg])/(