source: trunk/GSASIIpwd.py @ 5107

Last change on this file since 5107 was 5107, checked in by vondreele, 8 months ago

fix instrument parameter file prep for x-ray data in RMCProfile

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