source: trunk/GSASIIpwd.py @ 5106

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

modifications to handle sequential PDFfit analysis & plot the results.
Modifications to G2strMain & IO suggested by Conrad Gillard

  • 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-03 19:07:59 +0000 (Fri, 03 Dec 2021) $
10# $Author: vondreele $
11# $Revision: 5106 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 5106 2021-12-03 19:07:59Z 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: 5106 $"
38GSASIIpath.SetVersionNumber("$Revision: 5106 $")
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.close()
2470    return fname
2471   
2472def MakeBack(PWDdata,Name):
2473    Back = PWDdata['Background'][0]
2474    inst = PWDdata['Instrument Parameters'][0]
2475    if 'chebyschev-1' != Back[0]:
2476        return None
2477    Nback = Back[2]
2478    BackVals = Back[3:]
2479    fname = Name+'.back'
2480    fl = open(fname,'w')
2481    fl.write('%10d\n'%Nback)
2482    for val in BackVals:
2483        if 'T' in inst['Type'][1]:
2484            fl.write('%12.6g\n'%(float(val)))
2485        else:
2486            fl.write('%12.6g\n'%val)
2487    fl.close()
2488    return fname
2489
2490def findDup(Atoms):
2491    Dup = []
2492    Fracs = []
2493    for iat1,at1 in enumerate(Atoms):
2494        if any([at1[0] in dup for dup in Dup]):
2495            continue
2496        else:
2497            Dup.append([at1[0],])
2498            Fracs.append([at1[6],])
2499        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
2500            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
2501                Dup[-1] += [at2[0],]
2502                Fracs[-1]+= [at2[6],]
2503    return Dup,Fracs
2504
2505def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):   
2506   
2507    Meta = RMCPdict['metadata']
2508    Atseq = RMCPdict['atSeq']
2509    Supercell =  RMCPdict['SuperCell']
2510    generalData = Phase['General']
2511    Dups,Fracs = findDup(Phase['Atoms'])
2512    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
2513    ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
2514    Sample = PWDdata['Sample Parameters']
2515    Meta['temperature'] = Sample['Temperature']
2516    Meta['pressure'] = Sample['Pressure']
2517    Cell = generalData['Cell'][1:7]
2518    Trans = np.eye(3)*np.array(Supercell)
2519    newPhase = copy.deepcopy(Phase)
2520    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
2521    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans)
2522    GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0]
2523    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.
2524    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
2525    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2526    Natm = np.count_nonzero(Natm-1)
2527    Atoms = newPhase['Atoms']
2528    reset = False
2529   
2530    if ifSfracs:
2531        Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
2532        Natm = np.count_nonzero(Natm-1)
2533        Satoms = []
2534        for i in range(len(Atoms)//Natm):
2535            ind = i*Natm
2536            Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
2537        Natoms = []
2538        for satoms in Satoms:
2539            for idup,dup in enumerate(Dups):
2540                ldup = len(dup)
2541                natm = len(satoms)
2542                i = 0
2543                while i < natm:
2544                    if satoms[i][0] in dup:
2545                        atoms = satoms[i:i+ldup]
2546                        try:
2547                            atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
2548                            Natoms.append(atom)
2549                        except IndexError:      #what about vacancies?
2550                            if 'Va' not in Atseq:
2551                                reset = True
2552                                Atseq.append('Va')
2553                                RMCPdict['aTypes']['Va'] = 0.0
2554                            atom = atoms[0]
2555                            atom[1] = 'Va'
2556                            Natoms.append(atom)
2557                        i += ldup
2558                    else:
2559                       i += 1
2560    else:
2561        Natoms = Atoms
2562   
2563    NAtype = np.zeros(len(Atseq))
2564    for atom in Natoms:
2565        NAtype[Atseq.index(atom[1])] += 1
2566    NAstr = ['%6d'%i for i in NAtype]
2567    Cell = newPhase['General']['Cell'][1:7]
2568    if os.path.exists(Name+'.his6f'):
2569        os.remove(Name+'.his6f')
2570    if os.path.exists(Name+'.neigh'):
2571        os.remove(Name+'.neigh')
2572    fname = Name+'.rmc6f'
2573    fl = open(fname,'w')
2574    fl.write('(Version 6f format configuration file)\n')
2575    for item in Meta:
2576        fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item]))
2577    fl.write('Atom types present:                 %s\n'%'    '.join(Atseq))
2578    fl.write('Number of each atom type:       %s\n'%''.join(NAstr))
2579    fl.write('Number of atoms:                %d\n'%len(Natoms))
2580    fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2]))
2581    fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%(
2582            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
2583    A,B = G2lat.cell2AB(Cell,True)
2584    fl.write('Lattice vectors (Ang):\n')   
2585    for i in [0,1,2]:
2586        fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2]))
2587    fl.write('Atoms (fractional coordinates):\n')
2588    nat = 0
2589    for atm in Atseq:
2590        for iat,atom in enumerate(Natoms):
2591            if atom[1] == atm:
2592                nat += 1
2593                atcode = Atcodes[iat].split(':')
2594                cell = [0,0,0]
2595                if '+' in atcode[1]:
2596                    cell = eval(atcode[1].split('+')[1])
2597                fl.write('%6d%4s  [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%(       
2598                        nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2]))
2599    fl.close()
2600    return fname,reset
2601
2602def MakeBragg(PWDdata,Name,Phase):
2603    generalData = Phase['General']
2604    Vol = generalData['Cell'][7]
2605    Data = PWDdata['Data']
2606    Inst = PWDdata['Instrument Parameters'][0]
2607    Bank = int(Inst['Bank'][1])
2608    Sample = PWDdata['Sample Parameters']
2609    Scale = Sample['Scale'][0]
2610    if 'X' in Inst['Type'][1]:
2611        Scale *= 2.
2612    Limits = PWDdata['Limits'][1]
2613    Ibeg = np.searchsorted(Data[0],Limits[0])
2614    Ifin = np.searchsorted(Data[0],Limits[1])+1
2615    fname = Name+'.bragg'
2616    fl = open(fname,'w')
2617    fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol))
2618    if 'T' in Inst['Type'][0]:
2619        fl.write('%12s%12s\n'%('   TOF,ms','  I(obs)'))
2620        for i in range(Ibeg,Ifin-1):
2621            fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i]))
2622    else:
2623        fl.write('%12s%12s\n'%('   2-theta, deg','  I(obs)'))
2624        for i in range(Ibeg,Ifin-1):
2625            fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i]))       
2626    fl.close()
2627    return fname
2628
2629def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict):
2630    Meta = RMCPdict['metadata']
2631    Times = RMCPdict['runTimes']
2632    Atseq = RMCPdict['atSeq']
2633    Natoms = RMCPdict['NoAtoms']
2634    sumatms = np.sum(np.array([Natoms[iatm] for iatm in Natoms]))
2635    Isotope = RMCPdict['Isotope']
2636    Isotopes = RMCPdict['Isotopes']
2637    Atypes = RMCPdict['aTypes']
2638    if 'Va' in Atypes:
2639        Isotope['Va'] = 'Nat. Abund.'
2640        Isotopes['Va'] = {'Nat. Abund.':{'SL':[0.0,0.0]}}
2641    atPairs = RMCPdict['Pairs']
2642    Files = RMCPdict['files']
2643    BraggWt = RMCPdict['histogram'][1]
2644    inst = PWDdata['Instrument Parameters'][0]
2645    try:
2646        refList = PWDdata['Reflection Lists'][Name]['RefList']
2647    except KeyError:
2648        return 'Error - missing reflection list; you must do Refine first'
2649    dMin = refList[-1][4]
2650    gsasType = 'xray2'
2651    if 'T' in inst['Type'][1]:
2652        gsasType = 'gsas3'
2653    elif 'X' in inst['Type'][1]:
2654        XFF = G2elem.GetFFtable(Atseq)
2655        Xfl = open(Name+'.xray','w')
2656        for atm in Atseq:
2657            fa = XFF[atm]['fa']
2658            fb = XFF[atm]['fb']
2659            fc = XFF[atm]['fc']
2660            Xfl.write('%2s  %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%(
2661                    atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc))
2662        Xfl.close()
2663    lenA = len(Atseq)
2664    Pairs = []
2665    Ncoeff = []
2666    Nblen = [Isotopes[at][Isotope[at]]['SL'][0] for at in Atypes]
2667    for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]:
2668        Pairs += pair
2669    for pair in Pairs:
2670        pair = pair.replace(' ','')
2671        at1,at2 = pair.split('-')
2672        if at1 == 'Va' or at2 == 'Va':
2673            ncoef = 0.0
2674        else:
2675            ncoef = Isotopes[at1][Isotope[at1]]['SL'][0]*Natoms[at1]/sumatms
2676            ncoef *= Isotopes[at2][Isotope[at2]]['SL'][0]*Natoms[at2]/sumatms
2677        if at1 != at2:
2678            ncoef *= 2.
2679        Ncoeff += [ncoef,]
2680    pairMin = [atPairs[pair] if pair in atPairs else [0.0,0.,0.] for pair in Pairs ]
2681    maxMoves = [Atypes[atm] if atm in Atypes else 0.0 for atm in Atseq ]
2682    fname = Name+'.dat'
2683    fl = open(fname,'w')
2684    fl.write(' %% Hand edit the following as needed\n')
2685    fl.write('TITLE :: '+Name+'\n')
2686    fl.write('MATERIAL :: '+Meta['material']+'\n')
2687    fl.write('PHASE :: '+Meta['phase']+'\n')
2688    fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n')
2689    fl.write('INVESTIGATOR :: '+Meta['owner']+'\n')
2690    if RMCPdict.get('useGPU',False):
2691        fl.write('GPU_ACCELERATOR :: 0\n')
2692    minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin])
2693    minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin])
2694    maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin])
2695    fl.write('MINIMUM_DISTANCES ::   %s  Angstrom\n'%minHD)
2696    maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves])
2697    fl.write('MAXIMUM_MOVES ::   %s Angstrom\n'%maxMv)
2698    fl.write('R_SPACING ::  0.0200 Angstrom\n')
2699    fl.write('PRINT_PERIOD :: 100\n')
2700    fl.write('TIME_LIMIT ::     %.2f MINUTES\n'%Times[0])
2701    fl.write('SAVE_PERIOD ::    %.2f MINUTES\n'%Times[1])
2702    fl.write('\n')
2703    fl.write('ATOMS :: '+' '.join(Atseq)+'\n')
2704    fl.write('\n')
2705    fl.write('FLAGS ::\n')
2706    fl.write('  > NO_MOVEOUT\n')
2707    fl.write('  > NO_SAVE_CONFIGURATIONS\n')
2708    fl.write('  > NO_RESOLUTION_CONVOLUTION\n')
2709    fl.write('\n')
2710    fl.write('INPUT_CONFIGURATION_FORMAT ::  rmc6f\n')
2711    fl.write('SAVE_CONFIGURATION_FORMAT  ::  rmc6f\n')
2712    fl.write('IGNORE_HISTORY_FILE ::\n')
2713    fl.write('\n')
2714    fl.write('NEUTRON_COEFFICIENTS :: '+''.join(['%9.5f'%coeff for coeff in Ncoeff])+'\n')
2715    fl.write('DISTANCE_WINDOW ::\n')
2716    fl.write('  > MNDIST :: %s\n'%minD)
2717    fl.write('  > MXDIST :: %s\n'%maxD)
2718    if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']):
2719        fl.write('\n')
2720        fl.write('POTENTIALS ::\n')
2721        fl.write('  > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.'])
2722        fl.write('  > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n')
2723        if len(RMCPdict['Potentials']['Stretch']):
2724            fl.write('  > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search'])
2725            for bond in RMCPdict['Potentials']['Stretch']:
2726                fl.write('  > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2]))       
2727        if len(RMCPdict['Potentials']['Angles']):
2728            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
2729            for angle in RMCPdict['Potentials']['Angles']:
2730                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
2731                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
2732    if RMCPdict['useBVS']:
2733        fl.write('BVS ::\n')
2734        fl.write('  > ATOM :: '+' '.join(Atseq)+'\n')
2735        fl.write('  > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))
2736        oxid = []
2737        for val in RMCPdict['Oxid']:
2738            if len(val) == 3:
2739                oxid.append(val[0][1:])
2740            else:
2741                oxid.append(val[0][2:])
2742        fl.write('  > OXID :: %s\n'%' '.join(oxid))
2743        fl.write('  > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']]))
2744        fl.write('  > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']]))
2745        fl.write('  > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']]))       
2746        fl.write('  > SAVE :: 100000\n')
2747        fl.write('  > UPDATE :: 100000\n')
2748        if len(RMCPdict['Swap']):
2749            fl.write('\n')
2750            fl.write('SWAP_MULTI ::\n')
2751            for swap in RMCPdict['Swap']:
2752                try:
2753                    at1 = Atseq.index(swap[0])
2754                    at2 = Atseq.index(swap[1])
2755                except ValueError:
2756                    break
2757                fl.write('  > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2]))
2758       
2759    if len(RMCPdict['FxCN']):
2760        fl.write('FIXED_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['FxCN']))       
2761        for ifx,fxcn in enumerate(RMCPdict['FxCN']):
2762            try:
2763                at1 = Atseq.index(fxcn[0])
2764                at2 = Atseq.index(fxcn[1])
2765            except ValueError:
2766                break
2767            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]))
2768    if len(RMCPdict['AveCN']):
2769        fl.write('AVERAGE_COORDINATION_CONSTRAINTS ::  %d\n'%len(RMCPdict['AveCN']))
2770        for iav,avcn in enumerate(RMCPdict['AveCN']):
2771            try:
2772                at1 = Atseq.index(avcn[0])
2773                at2 = Atseq.index(avcn[1])
2774            except ValueError:
2775                break
2776            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]))
2777    for File in Files:
2778        if Files[File][0] and Files[File][0] != 'Select':
2779            if 'Xray' in File and 'F(Q)' in File:
2780                fqdata = open(Files[File][0],'r')
2781                lines = int(fqdata.readline()[:-1])
2782            fl.write('\n')
2783            fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_'))
2784            fl.write('  > FILENAME :: %s\n'%Files[File][0])
2785            fl.write('  > DATA_TYPE :: %s\n'%Files[File][2])
2786            fl.write('  > FIT_TYPE :: %s\n'%Files[File][2])
2787            if 'Xray' not in File:
2788                fl.write('  > START_POINT :: 1\n')
2789                fl.write('  > END_POINT :: 3000\n')
2790                fl.write('  > WEIGHT :: %.4f\n'%Files[File][1])
2791            fl.write('  > CONSTANT_OFFSET 0.000\n')
2792            fl.write('  > NO_FITTED_OFFSET\n')
2793            if RMCPdict['FitScale']:
2794                fl.write('  > FITTED_SCALE\n')
2795            else:
2796                fl.write('  > NO_FITTED_SCALE\n')
2797            if Files[File][3] !='RMC':
2798                fl.write('  > %s\n'%Files[File][3])
2799            if 'reciprocal' in File:
2800                fl.write('  > CONVOLVE ::\n')
2801                if 'Xray' in File:
2802                    fl.write('  > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines)
2803                    fl.write('  > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1]))
2804                    fl.write('  > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2))
2805                    fl.write('  > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1]))
2806    fl.write('\n')
2807    fl.write('BRAGG ::\n')
2808    fl.write('  > BRAGG_SHAPE :: %s\n'%gsasType)
2809    fl.write('  > RECALCUATE\n')
2810    fl.write('  > DMIN :: %.2f\n'%(dMin-0.02))
2811    fl.write('  > WEIGHT :: %10.3f\n'%BraggWt)
2812    fl.write('  > SCATTERING LENGTH :: '+''.join(['%8.4f'%blen for blen in Nblen])+'\n')
2813    fl.write('\n')
2814    fl.write('END  ::\n')
2815    fl.close()
2816    return fname
2817
2818# def FindBonds(Phase,RMCPdict):
2819#     generalData = Phase['General']
2820#     cx,ct,cs,cia = generalData['AtomPtrs']
2821#     atomData = Phase['Atoms']
2822#     Res = 'RMC'
2823#     if 'macro' in generalData['Type']:
2824#         Res = atomData[0][ct-3]
2825#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2826#     Pairs = RMCPdict['Pairs']   #dict!
2827#     BondList = []
2828#     notNames = []
2829#     for FrstName in AtDict:
2830#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
2831#         Atyp1 = AtDict[FrstName]
2832#         if 'Va' in Atyp1:
2833#             continue
2834#         for nbr in nbrs:
2835#             Atyp2 = AtDict[nbr[0]]
2836#             if 'Va' in Atyp2:
2837#                 continue
2838#             try:
2839#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
2840#             except KeyError:
2841#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
2842#             if any(bndData):
2843#                 if bndData[0] <= nbr[1] <= bndData[1]:
2844#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
2845#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
2846#                     if bondStr not in BondList and revbondStr not in BondList:
2847#                         BondList.append(bondStr)
2848#         notNames.append(FrstName)
2849#     return Res,BondList
2850
2851# def FindAngles(Phase,RMCPdict):
2852#     generalData = Phase['General']
2853#     Cell = generalData['Cell'][1:7]
2854#     Amat = G2lat.cell2AB(Cell)[0]
2855#     cx,ct,cs,cia = generalData['AtomPtrs']
2856#     atomData = Phase['Atoms']
2857#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
2858#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
2859#     Angles = RMCPdict['Angles']
2860#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
2861#     AngleList = []
2862#     for MidName in AtDict:
2863#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
2864#         if len(nbrs) < 2: #need 2 neighbors to make an angle
2865#             continue
2866#         Atyp2 = AtDict[MidName]
2867#         for i,nbr1 in enumerate(nbrs):
2868#             Atyp1 = AtDict[nbr1[0]]
2869#             for j,nbr3 in enumerate(nbrs[i+1:]):
2870#                 Atyp3 = AtDict[nbr3[0]]
2871#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
2872#                 try:
2873#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
2874#                 except KeyError:
2875#                     try:
2876#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
2877#                     except KeyError:
2878#                         continue
2879#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
2880#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
2881#                 if angData[0] <= calAngle <= angData[1]:
2882#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
2883#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
2884#                     if angStr not in AngleList and revangStr not in AngleList:
2885#                         AngleList.append(angStr)
2886#     return AngleList
2887
2888# def GetSqConvolution(XY,d):
2889
2890#     n = XY.shape[1]
2891#     snew = np.zeros(n)
2892#     dq = np.zeros(n)
2893#     sold = XY[1]
2894#     q = XY[0]
2895#     dq[1:] = np.diff(q)
2896#     dq[0] = dq[1]
2897   
2898#     for j in range(n):
2899#         for i in range(n):
2900#             b = abs(q[i]-q[j])
2901#             t = q[i]+q[j]
2902#             if j == i:
2903#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
2904#             else:
2905#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
2906#         snew[j] /= np.pi*q[j]
2907   
2908#     snew[0] = snew[1]
2909#     return snew
2910
2911# def GetMaxSphere(pdbName):
2912#     try:
2913#         pFil = open(pdbName,'r')
2914#     except FileNotFoundError:
2915#         return None
2916#     while True:
2917#         line = pFil.readline()
2918#         if 'Boundary' in line:
2919#             line = line.split()[3:]
2920#             G = np.array([float(item) for item in line])
2921#             G = np.reshape(G,(3,3))**2
2922#             G = nl.inv(G)
2923#             pFil.close()
2924#             break
2925#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
2926#     return min(dspaces)
2927
2928def findfullrmc():
2929    '''Find where fullrmc is installed. Tries the following:
2930   
2931         1. Returns the Config var 'fullrmc_exec', if defined. No check
2932            is done that the interpreter has fullrmc
2933         2. The current Python interpreter if fullrmc can be imported
2934            and fullrmc is version 5+
2935         3. The path is checked for a fullrmc image as named by Bachir
2936
2937    :returns: the full path to a python executable that is assumed to
2938      have fullrmc installed or None, if it was not found.
2939    '''
2940    is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
2941    if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
2942            GSASIIpath.GetConfigValue('fullrmc_exec')):
2943        return GSASIIpath.GetConfigValue('fullrmc_exec')
2944    try:
2945        import fullrmc
2946        if int(fullrmc.__version__.split('.')[0]) >= 5:
2947            return sys.executable
2948    except:
2949        pass
2950    pathlist = os.environ["PATH"].split(os.pathsep)
2951    for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
2952                  os.path.split(sys.executable)[0]):
2953        if p not in pathlist: pathlist.insert(0,p)
2954    import glob
2955    for p in pathlist:
2956        if sys.platform == "darwin":
2957            lookfor = "fullrmc*macOS*i386-64bit"
2958        elif sys.platform == "win32":
2959            lookfor = "fullrmc*.exe"
2960        else:
2961            lookfor = "fullrmc*"
2962        fl = glob.glob(lookfor)
2963        if len(fl) > 0:
2964            return os.path.abspath(sorted(fl)[0])
2965       
2966def findPDFfit():
2967    '''Find if PDFfit2 is installed (may be local to GSAS-II). Does the following:
2968    :returns: the full path to a python executable or None, if it was not found.
2969    '''
2970    try:
2971        sys.path.append('%s'%GSASIIpath.path2GSAS2) 
2972        from diffpy.pdffit2 import PdfFit
2973        return sys.executable
2974    except:
2975        return None
2976   
2977def GetPDFfitAtomVar(Phase,RMCPdict):
2978    ''' Find dict of independent "@n" variables for PDFfit in atom constraints
2979    '''
2980    General = Phase['General']
2981    Atoms = Phase['Atoms']
2982    cx,ct,cs,cia = General['AtomPtrs']
2983    AtomVar = RMCPdict['AtomVar']
2984    varnames = []
2985    for iat,atom in enumerate(RMCPdict['AtomConstr']):
2986        for it,item in enumerate(atom):
2987            if it > 1 and item:
2988                itnum = item.split('@')[1]
2989                varname = '@%s'%itnum
2990                varnames.append(varname)
2991                if it < 6:
2992                    if varname not in AtomVar:
2993                        AtomVar[varname] = 0.0      #put ISODISTORT mode displ here?
2994                else:
2995                    for i in range(3):
2996                        if varname not in AtomVar:
2997                            AtomVar[varname] = Atoms[iat][cia+i+2]
2998    varnames = set(varnames)
2999    for name in list(AtomVar.keys()):       #clear out unused parameters
3000        if name not in varnames:
3001            del AtomVar[name]
3002   
3003def MakePDFfitAtomsFile(Phase,RMCPdict):
3004    '''Make the PDFfit atoms file
3005    '''
3006    General = Phase['General']
3007    fName = General['Name']+'-PDFfit.stru'
3008    fName = fName.replace(' ','_')
3009    fatm = open(fName,'w')
3010    fatm.write('title  structure of '+General['Name']+'\n')
3011    fatm.write('format pdffit\n')
3012    fatm.write('scale   1.000000\n')    #fixed
3013    # if RMCPdict['shape'] == 'sphere':
3014    #     sharp = '%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0])
3015    # else:
3016    #     sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3017    sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut'])
3018    fatm.write('sharp '+sharp)
3019    shape = ''
3020    if RMCPdict['shape'] == 'sphere' and RMCPdict['spdiameter'][0] > 0.:
3021        shape = '   sphere, %10.6f\n'%RMCPdict['spdiameter'][0]
3022    elif RMCPdict['stepcut'] > 0.:
3023        shape = 'stepcut, %10.6f\n'%RMCPdict['stepcut']
3024    if shape:
3025        fatm.write('shape  '+shape)
3026    fatm.write('spcgr   %s\n'%General['SGData']['SpGrp'].replace(' ',''))
3027    cell = General['Cell'][1:7]
3028    fatm.write('cell  %10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n'%(
3029        cell[0],cell[1],cell[2],cell[3],cell[4],cell[5]))
3030    fatm.write('dcell '+5*'  0.000000,'+'  0.000000\n') 
3031    Atoms = Phase['Atoms']
3032    fatm.write('ncell %8d,%8d,%8d,%10d\n'%(1,1,1,len(Atoms)))
3033    fatm.write('atoms\n')
3034    cx,ct,cs,cia = General['AtomPtrs']
3035    for atom in Atoms:
3036        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]))
3037        fatm.write('    '+'%18.8f%18.8f%18.8f%13.4f\n'%(0.,0.,0.,0.))
3038        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(atom[cia+2],atom[cia+3],atom[cia+4]))
3039        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.,))
3040        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(atom[cia+5],atom[cia+6],atom[cia+7]))
3041        fatm.write('    '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.))
3042    fatm.close()
3043   
3044def MakePDFfitRunFile(Phase,RMCPdict):
3045    '''Make the PDFfit python run file
3046    '''
3047   
3048    def GetCellConstr(SGData):
3049        if SGData['SGLaue'] in ['m3', 'm3m']:
3050            return [1,1,1,0,0,0]
3051        elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']:
3052            return [1,1,2,0,0,0]
3053        elif SGData['SGLaue'] in ['3R','3mR']:
3054            return [1,1,1,2,2,2]
3055        elif SGData['SGLaue'] == 'mmm':
3056            return [1,2,3,0,0,0]
3057        elif SGData['SGLaue'] == '2/m':
3058            if SGData['SGUniq'] == 'a':
3059                return [1,2,3,4,0,0]
3060            elif SGData['SGUniq'] == 'b':
3061                return [1,2,3,0,4,0]
3062            elif SGData['SGUniq'] == 'c':
3063                return [1,2,3,0,0,4]
3064        else:
3065            return [1,2,3,4,5,6]
3066       
3067    General = Phase['General']
3068    Cell = General['Cell'][1:7]
3069    G2path = GSASIIpath.path2GSAS2
3070    rundata = '''
3071#!/usr/bin/env python
3072# -*- coding: utf-8 -*-
3073import sys
3074sys.path.append('%s')
3075from diffpy.pdffit2 import PdfFit
3076pf = PdfFit()
3077'''%G2path
3078    Nd = 0
3079    Np = 0
3080    parms = {}
3081    if 'sequential' in RMCPdict['refinement']:
3082        Np = 3
3083        rundata += '#sequential data here\n'
3084    else:
3085        for file in RMCPdict['files']:
3086            if 'Select' in RMCPdict['files'][file][0]:
3087                continue
3088            if 'Neutron' in file:
3089                Nd += 1
3090                dType = 'Ndata'
3091            else:
3092                Nd += 1
3093                dType = 'Xdata'
3094            rundata += "pf.read_data('%s', '%s', 30.0, %.4f)\n"%(RMCPdict['files'][file][0],dType[0],RMCPdict[dType]['qdamp'][0])
3095            rundata += 'pf.setdata(%d)\n'%Nd
3096            rundata += 'pf.pdfrange(%d, %6.2f, %6.2f)\n'%(Nd,RMCPdict[dType]['Fitrange'][0],RMCPdict[dType]['Fitrange'][1])
3097            for item in ['dscale','qdamp','qbroad']:
3098                if RMCPdict[dType][item][1]:
3099                    Np += 1
3100                    rundata += 'pf.constrain(pf.%s(),"@%d")\n'%(item,Np)
3101                    parms[Np] = [RMCPdict[dType][item][0],item]
3102    fName = General['Name']+'-PDFfit.stru'
3103    fName = fName.replace(' ','_')
3104    Np = 9
3105    rundata += "pf.read_struct('%s')\n"%(fName)
3106    for item in ['delta1','delta2','sratio']:
3107        if RMCPdict[item][1]:
3108            Np += 1
3109            rundata += 'pf.constrain(pf.%s,"@%d")\n'%(item,Np)
3110            parms[Np] = [RMCPdict[item][0],item]
3111    if 'sphere' in RMCPdict['shape'][0] and RMCPdict['spdiameter'][1]:
3112        Np += 1
3113        rundata += 'pf.constrain(pf.spdiameter,"@%d")\n'%Np
3114        parms[Np] = [RMCPdict['spdiameter'][0],'spdiameter']
3115   
3116    if RMCPdict['cellref']:
3117        cellconst = GetCellConstr(RMCPdict['SGData'])
3118        used = []
3119        cellNames = ['a','b','c','alpha','beta','gamma']
3120        for ic in range(6):
3121            if cellconst[ic]:
3122                rundata += 'pf.constrain(pf.lat(%d), "@%d")\n'%(ic+1,Np+cellconst[ic])
3123                if cellconst[ic] not in used:
3124                    parms[Np+cellconst[ic]] = [Cell[ic],cellNames[ic]]
3125                used.append(cellconst[ic])
3126#Atom constraints here -------------------------------------------------------
3127    AtomVar = RMCPdict['AtomVar']
3128    used = []
3129    for iat,atom in enumerate(RMCPdict['AtomConstr']):
3130        for it,item in enumerate(atom):
3131            names = ['pf.x(%d)'%(iat+1),'pf.y(%d)'%(iat+1),'pf.z(%d)'%(iat+1),'pf.occ(%d)'%(iat+1)]
3132            if it > 1 and item:
3133                itnum = item.split('@')[1]
3134                if it < 6:
3135                    rundata += 'pf.constrain(%s,"%s")\n'%(names[it-2],item)
3136                    if itnum not in used:
3137                        parms[itnum] = [AtomVar['@%s'%itnum],names[it-2].split('.')[1]]
3138                        used.append(itnum)
3139                else:
3140                    uijs = ['pf.u11(%d)'%(iat+1),'pf.u22(%d)'%(iat+1),'pf.u33(%d)'%(iat+1)]     
3141                    for i in range(3):
3142                        rundata += 'pf.constrain(%s,"%s")\n'%(uijs[i],item)
3143                        if itnum not in used:
3144                            parms[itnum] = [AtomVar['@%s'%itnum],uijs[i].split('.')[1]]
3145                            used.append(itnum)
3146                           
3147    if 'sequential' in RMCPdict['refinement']:
3148        rundata += '#parameters here\n'
3149    else:       
3150# set parameter values
3151        RMCPdict['Parms'] = parms
3152        for iprm in parms:
3153            rundata += 'pf.setpar(%s,%.6f)\n'%(iprm,parms[iprm][0])
3154                       
3155# Refine & Save results ---------------------------------------------------------------   
3156    rundata += 'pf.refine()\n'
3157    if 'sequential' in RMCPdict['refinement']:
3158        fName = 'Sequential_PDFfit'
3159    else:
3160        fName = General['Name'].replace(' ','_')+'-PDFfit'
3161    Nd = 0   
3162    for file in RMCPdict['files']:
3163        if 'Select' in RMCPdict['files'][file][0]:
3164            continue
3165        Nd += 1
3166        rundata += 'pf.save_pdf(%d, "%s")\n'%(Nd,fName+file[0]+'.fgr')
3167       
3168    rundata += 'pf.save_struct(1, "%s")\n'%(fName+'.rstr')
3169    rundata += 'pf.save_res("%s")\n'%(fName+'.res')
3170 
3171    if 'sequential' in RMCPdict['refinement']:
3172        rfile = open('Seq_PDFfit_template.py','w')
3173    else:
3174        rfile = open(fName+'.py','w')
3175    rfile.writelines(rundata)
3176    rfile.close()
3177   
3178    return fName+'.py'
3179   
3180def UpdatePDFfit(Phase,RMCPdict):
3181    ''' Updates various PDFfit parameters held in GSAS-II
3182    '''
3183   
3184    General = Phase['General']
3185    if RMCPdict['refinement'] == 'normal':
3186        fName = General['Name']+'-PDFfit.rstr'
3187        try:
3188            rstr = open(fName.replace(' ','_'),'r')
3189        except FileNotFoundError:
3190            return [fName,'Not found - PDFfit failed']
3191        lines = rstr.readlines()
3192        rstr.close()
3193        header = [line[:-1].split(' ',1) for line in lines[:7]]
3194        resdict = dict(header)
3195        for item in ['scale','sharp','cell']:
3196            resdict[item] = [float(val) for val in resdict[item].split(',')]
3197        General['Cell'][1:7] = resdict['cell']
3198        for inam,name in enumerate(['delta2','delta1','sratio']):
3199            RMCPdict[name][0] = resdict['sharp'][inam]
3200        if 'shape' in resdict:
3201            if 'sphere' in resdict['shape']:
3202                RMCPdict['spdiameter'][0] = resdict['shape'][-1]
3203            else:
3204                RMCPdict['stepcut'][0] = resdict['shape'][-1]
3205        cx,ct,cs,ci = G2mth.getAtomPtrs(Phase)     
3206        Atoms = Phase['Atoms']
3207        atmBeg = 0
3208        for line in lines:
3209            atmBeg += 1
3210            if 'atoms' in line:
3211                break
3212        for atom in Atoms:
3213            atstr = lines[atmBeg][:-1].split()
3214            Uiistr = lines[atmBeg+2][:-1].split()
3215            Uijstr = lines[atmBeg+4][:-1].split()
3216            atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])]
3217            atom[ci] = 'A'
3218            atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])]
3219            atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])]
3220            atmBeg += 6
3221           
3222    fName = 'Sequential_PDFfit.res'
3223    if RMCPdict['refinement'] == 'normal':
3224        fName = General['Name']+'-PDFfit.res'
3225    try:
3226        res = open(fName.replace(' ','_'),'r')
3227    except FileNotFoundError:
3228        return [fName,'Not found - PDFfit failed']
3229    lines = res.readlines()
3230    res.close()
3231    Ibeg = False
3232    resline = ''
3233    for iline,line in enumerate(lines):
3234        if 'Refinement parameters' in line:
3235            Ibeg = True
3236            continue
3237        if Ibeg:
3238            if '---------' in line:
3239                break
3240            resline += line[:-1]
3241    for iline,line in enumerate(lines):
3242        if 'Rw - ' in line:
3243            Rwp = float(line.split(':')[1])
3244    results = resline.replace('(','').split(')')[:-1]
3245    results = ['@'+result.lstrip() for result in results]
3246    results = [item.split() for item in results]
3247    if RMCPdict['refinement'] == 'normal':
3248        results = dict([[item[0][:-1],float(item[1])] for item in results if item[0][:-1] in RMCPdict['AtomVar']])
3249        RMCPdict['AtomVar'].update(results)
3250        return None
3251    else:   #sequential
3252        newParms = dict([[item[0][1:-1],[float(item[1]),float(item[2])]] for item in results])
3253        return newParms,Rwp
3254       
3255def MakefullrmcRun(pName,Phase,RMCPdict):
3256    '''Creates a script to run fullrmc. Returns the name of the file that was
3257    created.
3258    '''
3259    BondList = {}
3260    for k in RMCPdict['Pairs']:
3261        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
3262            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
3263    AngleList = []
3264    for angle in RMCPdict['Angles']:
3265        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
3266            continue
3267        for i in (0,1,2):
3268            angle[i] = angle[i].strip()
3269        AngleList.append(angle)
3270    # rmin = RMCPdict['min Contact']
3271    cell = Phase['General']['Cell'][1:7]
3272    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
3273    cx,ct,cs,cia = Phase['General']['AtomPtrs']
3274    atomsList = []
3275    for atom in Phase['Atoms']:
3276        el = ''.join([i for i in atom[ct] if i.isalpha()])
3277        atomsList.append([el] + atom[cx:cx+4])
3278    projDir,projName = os.path.split(os.path.abspath(pName))
3279    scrname = pName+'-fullrmc.py'
3280    restart = '%s_restart.pdb'%pName
3281    Files = RMCPdict['files']
3282    rundata = ''
3283    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
3284    rundata += '# created in '+__file__+" v"+filversion.split()[1]
3285    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
3286    rundata += '''
3287# fullrmc imports (all that are potentially useful)
3288import os,glob
3289import time
3290import pickle
3291import numpy as np
3292from fullrmc.Core import Collection
3293from fullrmc.Engine import Engine
3294import fullrmc.Constraints.PairDistributionConstraints as fPDF
3295from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
3296from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
3297from fullrmc.Constraints.BondConstraints import BondConstraint
3298from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
3299from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
3300from fullrmc.Generators.Swaps import SwapPositionsGenerator
3301# utility routines
3302def writeHeader(ENGINE,statFP):
3303    'header for stats file'
3304    statFP.write('generated-steps, total-error, ')
3305    for c in ENGINE.constraints:
3306        statFP.write(c.constraintName)
3307        statFP.write(', ')
3308    statFP.write('\\n')
3309    statFP.flush()
3310   
3311def writeCurrentStatus(ENGINE,statFP,plotF):
3312    'line in stats file & current constraint plots'
3313    statFP.write(str(ENGINE.generated))
3314    statFP.write(', ')
3315    statFP.write(str(ENGINE.totalStandardError))
3316    statFP.write(', ')
3317    for c in ENGINE.constraints:
3318        statFP.write(str(c.standardError))
3319        statFP.write(', ')
3320    statFP.write('\\n')
3321    statFP.flush()
3322    mpl.use('agg')
3323    fp = open(plotF,'wb')
3324    for c in ENGINE.constraints:
3325        p = c.plot(show=False)
3326        p[0].canvas.draw()
3327        image = p[0].canvas.buffer_rgba()
3328        pickle.dump(c.constraintName,fp)
3329        pickle.dump(np.array(image),fp)
3330    fp.close()
3331
3332def calcRmax(ENGINE):
3333    'from Bachir, works for non-othorhombic'
3334    a,b,c = ENGINE.basisVectors
3335    lens = []
3336    ts    = np.linalg.norm(np.cross(a,b))/2
3337    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3338    ts = np.linalg.norm(np.cross(b,c))/2
3339    lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3340    ts = np.linalg.norm(np.cross(a,c))/2
3341    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3342    return min(lens)
3343'''
3344    rundata += '''
3345### When True, erases an existing engine to provide a fresh start
3346FRESH_START = {:}
3347dirName = "{:}"
3348prefix = "{:}"
3349project = prefix + "-fullrmc"
3350time0 = time.time()
3351'''.format(RMCPdict['ReStart'][0],projDir,projName)
3352   
3353    rundata += '# setup structure\n'
3354    rundata += 'cell = ' + str(cell) + '\n'
3355    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3356    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
3357    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3358
3359    rundata += '\n# initialize engine\n'
3360    rundata += '''
3361engineFileName = os.path.join(dirName, project + '.rmc')
3362projectStats = os.path.join(dirName, project + '.stats')
3363projectPlots = os.path.join(dirName, project + '.plots')
3364pdbFile = os.path.join(dirName, project + '_restart.pdb')
3365# check Engine exists if so (and not FRESH_START) load it
3366# otherwise build it
3367ENGINE = Engine(path=None)
3368if not ENGINE.is_engine(engineFileName) or FRESH_START:
3369    ## create structure
3370    ENGINE = Engine(path=engineFileName, freshStart=True)
3371    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
3372                                 atoms      = atomList,
3373                                 unitcellBC = cell,
3374                                 supercell  = supercell)
3375'''   
3376    import atmdata
3377    # rundata += '    # conversion factors (may be needed)\n'
3378    # rundata += '    sumCiBi2 = 0.\n'
3379    # for elem in Phase['General']['AtomTypes']:
3380    #     rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3381    #     rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3382    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3383    # settings that require a new Engine
3384    for File in Files:
3385        filDat = RMCPdict['files'][File]
3386        if not os.path.exists(filDat[0]): continue
3387        sfwt = 'neutronCohb'
3388        if 'Xray' in File:
3389            sfwt = 'atomicNumber'
3390        if 'G(r)' in File:
3391            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3392            if filDat[3] == 0:
3393                #rundata += '''    # read and xform G(r) as defined in RMCProfile
3394    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3395                #rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
3396                #rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3397                rundata += '    # G(r) as defined in RMCProfile\n'
3398                rundata += '    GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3399            elif filDat[3] == 1:
3400                rundata += '    # This is G(r) as defined in PDFFIT\n'
3401                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3402            elif filDat[3] == 2:
3403                rundata += '    # This is g(r)\n'
3404                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3405            else:
3406                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3407            rundata += '    ENGINE.add_constraints([GofR])\n'
3408            rundata += '    GofR.set_limits((None, calcRmax(ENGINE)))\n'
3409        elif '(Q)' in File:
3410            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3411            if filDat[3] == 0:
3412                rundata += '    # F(Q) as defined in RMCProfile\n'
3413                #rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
3414                if filDat[4]:
3415                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3416                rundata += '    SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3417            elif filDat[3] == 1:
3418                rundata += '    # S(Q) as defined in PDFFIT\n'
3419                rundata += '    SOQ[1] -= 1\n'
3420                if filDat[4]:
3421                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3422                rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3423            else:
3424                raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3425            rundata += '    ENGINE.add_constraints([SofQ])\n'
3426        else:
3427            print('What is this?')
3428    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
3429    if BondList:
3430        rundata += '''    B_CONSTRAINT   = BondConstraint()
3431    ENGINE.add_constraints(B_CONSTRAINT)
3432    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3433'''
3434        for pair in BondList:
3435            e1,e2 = pair.split('-')
3436            d1,d2 = BondList[pair]
3437            if d1 == 0: continue
3438            if d2 == 0:
3439                print('Ignoring min distance without maximum')
3440                #rundata += '            ("element","{}","{}",{}),\n'.format(
3441                #                        e1.strip(),e2.strip(),d1)
3442            else:
3443                rundata += '            ("element","{}","{}",{},{}),\n'.format(
3444                                        e1.strip(),e2.strip(),d1,d2)
3445        rundata += '             ])\n'
3446    if AngleList:
3447        rundata += '''    A_CONSTRAINT   = BondsAngleConstraint()
3448    ENGINE.add_constraints(A_CONSTRAINT)
3449    A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3450'''
3451        for item in AngleList:
3452            rundata += ('            '+
3453               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3454        rundata += '             ])\n'
3455    rundata += '''
3456    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3457    ENGINE.save()
3458else:
3459    ENGINE = ENGINE.load(path=engineFileName)
3460
3461ENGINE.set_log_file(os.path.join(dirName,prefix))
3462'''
3463    if RMCPdict['Swaps']:
3464        rundata += '\n#set up for site swaps\n'
3465        rundata += 'aN = ENGINE.allNames\n'
3466        rundata += 'SwapGen = {}\n'
3467        for swap in RMCPdict['Swaps']:
3468            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3469            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3470            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3471        rundata += '    for swaps in SwapGen:\n'
3472        rundata += '        AB = swaps.split("-")\n'
3473        rundata += '        ENGINE.set_groups_as_atoms()\n'
3474        rundata += '        for g in ENGINE.groups:\n'
3475        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
3476        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
3477        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
3478        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
3479        rundata += '            sProb = SwapGen[swaps][2]\n'
3480    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3481    # rundata += 'wtDict = {}\n'
3482    # for File in Files:
3483    #     filDat = RMCPdict['files'][File]
3484    #     if not os.path.exists(filDat[0]): continue
3485    #     if 'Xray' in File:
3486    #         sfwt = 'atomicNumber'
3487    #     else:
3488    #         sfwt = 'neutronCohb'
3489    #     if 'G(r)' in File:
3490    #         typ = 'Pair'
3491    #     elif '(Q)' in File:
3492    #         typ = 'Struct'
3493    #     rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3494    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
3495    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
3496    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3497    rundata += '        c.set_limits((None,calcRmax(ENGINE)))\n'
3498    if RMCPdict['FitScale']:
3499        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3500    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3501    if RMCPdict['FitScale']:
3502        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
3503        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3504    # torsions difficult to implement, must be internal to cell & named with
3505    # fullrmc atom names
3506    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
3507    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
3508    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
3509    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3510    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3511    #     for torsion in RMCPdict['Torsions']:
3512    #         rundata += '    %s\n'%str(tuple(torsion))
3513    #     rundata += '        ]})\n'           
3514    rundata += '''
3515if FRESH_START:
3516    # initialize engine with one step to get starting config energetics
3517    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3518    statFP = open(projectStats,'w')
3519    writeHeader(ENGINE,statFP)
3520    writeCurrentStatus(ENGINE,statFP,projectPlots)
3521else:
3522    statFP = open(projectStats,'a')
3523
3524# setup runs for fullrmc
3525'''
3526    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3527    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3528    rundata += '    ENGINE.set_groups_as_atoms()\n'
3529    rundata += '    expected = ENGINE.generated+steps\n'
3530   
3531    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3532    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3533    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
3534    rundata += 'statFP.close()\n'
3535    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3536    rfile = open(scrname,'w')
3537    rfile.writelines(rundata)
3538    rfile.close()
3539    return scrname
3540   
3541def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3542    bondDist = []
3543    Cell = general['Cell'][1:7]
3544    Supercell =  RMCPdict['SuperCell']
3545    Trans = np.eye(3)*np.array(Supercell)
3546    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3547    Amat,Bmat = G2lat.cell2AB(Cell)
3548    indices = (-1,0,1)
3549    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3550    for bonds in bondList:
3551        Oxyz = np.array(Atoms[bonds[0]][1:])
3552        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
3553        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3554        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3555        for dx in Dx.T:
3556            bondDist.append(np.min(dx))
3557    return np.array(bondDist)
3558   
3559def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3560    bondAngles = []
3561    Cell = general['Cell'][1:7]
3562    Supercell =  RMCPdict['SuperCell']
3563    Trans = np.eye(3)*np.array(Supercell)
3564    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3565    Amat,Bmat = G2lat.cell2AB(Cell)
3566    indices = (-1,0,1)
3567    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3568    for angle in angleList:
3569        Oxyz = np.array(Atoms[angle[0]][1:])
3570        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
3571        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])       
3572        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3573        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3574        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3575        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3576        iDAx = np.argmin(DAx,axis=0)
3577        iDBx = np.argmin(DBx,axis=0)
3578        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3579            DAv = DAxV[iA,i]/DAx[iA,i]
3580            DBv = DBxV[iB,i]/DBx[iB,i]
3581            bondAngles.append(npacosd(np.sum(DAv*DBv)))
3582    return np.array(bondAngles)
3583   
3584   
3585#### Reflectometry calculations ################################################################################
3586def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3587    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3588   
3589    class RandomDisplacementBounds(object):
3590        """random displacement with bounds"""
3591        def __init__(self, xmin, xmax, stepsize=0.5):
3592            self.xmin = xmin
3593            self.xmax = xmax
3594            self.stepsize = stepsize
3595   
3596        def __call__(self, x):
3597            """take a random step but ensure the new position is within the bounds"""
3598            while True:
3599                # this could be done in a much more clever way, but it will work for example purposes
3600                steps = self.xmax-self.xmin
3601                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
3602                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3603                    break
3604            return xnew
3605   
3606    def GetModelParms():
3607        parmDict = {}
3608        varyList = []
3609        values = []
3610        bounds = []
3611        parmDict['dQ type'] = data['dQ type']
3612        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
3613        for parm in ['Scale','FltBack']:
3614            parmDict[parm] = data[parm][0]
3615            if data[parm][1]:
3616                varyList.append(parm)
3617                values.append(data[parm][0])
3618                bounds.append(Bounds[parm])
3619        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3620        parmDict['nLayers'] = len(parmDict['Layer Seq'])
3621        for ilay,layer in enumerate(data['Layers']):
3622            name = layer['Name']
3623            cid = str(ilay)+';'
3624            parmDict[cid+'Name'] = name
3625            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3626                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3627                if layer.get(parm,[0.,False])[1]:
3628                    varyList.append(cid+parm)
3629                    value = layer[parm][0]
3630                    values.append(value)
3631                    if value:
3632                        bound = [value*Bfac,value/Bfac]
3633                    else:
3634                        bound = [0.,10.]
3635                    bounds.append(bound)
3636            if name not in ['vacuum','unit scatter']:
3637                parmDict[cid+'rho'] = Substances[name]['Scatt density']
3638                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3639        return parmDict,varyList,values,bounds
3640   
3641    def SetModelParms():
3642        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3643        if 'Scale' in varyList:
3644            data['Scale'][0] = parmDict['Scale']
3645            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
3646        G2fil.G2Print (line)
3647        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3648        if 'FltBack' in varyList:
3649            data['FltBack'][0] = parmDict['FltBack']
3650            line += ' esd: %15.3g'%(sigDict['FltBack'])
3651        G2fil.G2Print (line)
3652        for ilay,layer in enumerate(data['Layers']):
3653            name = layer['Name']
3654            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3655            cid = str(ilay)+';'
3656            line = ' '
3657            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3658            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3659            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3660                if parm in layer:
3661                    if parm == 'Rough':
3662                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
3663                    else:
3664                        layer[parm][0] = parmDict[cid+parm]
3665                    line += ' %s: %.3f'%(parm,layer[parm][0])
3666                    if cid+parm in varyList:
3667                        line += ' esd: %.3g'%(sigDict[cid+parm])
3668            G2fil.G2Print (line)
3669            G2fil.G2Print (line2)
3670   
3671    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3672        parmDict.update(zip(varyList,values))
3673        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3674        return M
3675   
3676    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3677        parmDict.update(zip(varyList,values))
3678        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3679        return np.sum(M**2)
3680   
3681    def getREFD(Q,Qsig,parmDict):
3682        Ic = np.ones_like(Q)*parmDict['FltBack']
3683        Scale = parmDict['Scale']
3684        Nlayers = parmDict['nLayers']
3685        Res = parmDict['Res']
3686        depth = np.zeros(Nlayers)
3687        rho = np.zeros(Nlayers)
3688        irho = np.zeros(Nlayers)
3689        sigma = np.zeros(Nlayers)
3690        for ilay,lay in enumerate(parmDict['Layer Seq']):
3691            cid = str(lay)+';'
3692            depth[ilay] = parmDict[cid+'Thick']
3693            sigma[ilay] = parmDict[cid+'Rough']
3694            if parmDict[cid+'Name'] == u'unit scatter':
3695                rho[ilay] = parmDict[cid+'DenMul']
3696                irho[ilay] = parmDict[cid+'iDenMul']
3697            elif 'vacuum' != parmDict[cid+'Name']:
3698                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3699                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3700            if cid+'Mag SLD' in parmDict:
3701                rho[ilay] += parmDict[cid+'Mag SLD']
3702        if parmDict['dQ type'] == 'None':
3703            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3704        elif 'const' in parmDict['dQ type']:
3705            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
3706        else:       #dQ/Q in data
3707            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3708        Ic += AB*Scale
3709        return Ic
3710       
3711    def estimateT0(takestep):
3712        Mmax = -1.e-10
3713        Mmin = 1.e10
3714        for i in range(100):
3715            x0 = takestep(values)
3716            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3717            Mmin = min(M,Mmin)
3718            MMax = max(M,Mmax)
3719        return 1.5*(MMax-Mmin)
3720
3721    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3722    if data.get('2% weight'):
3723        wt = 1./(0.02*Io)**2
3724    Qmin = Limits[1][0]
3725    Qmax = Limits[1][1]
3726    wtFactor = ProfDict['wtFactor']
3727    Bfac = data['Toler']
3728    Ibeg = np.searchsorted(Q,Qmin)
3729    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
3730    Ic[:] = 0
3731    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3732              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3733    parmDict,varyList,values,bounds = GetModelParms()
3734    Msg = 'Failed to converge'
3735    if varyList:
3736        if data['Minimizer'] == 'LMLS': 
3737            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3738                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3739            parmDict.update(zip(varyList,result[0]))
3740            chisq = np.sum(result[2]['fvec']**2)
3741            ncalc = result[2]['nfev']
3742            covM = result[1]
3743            newVals = result[0]
3744        elif data['Minimizer'] == 'Basin Hopping':
3745            xyrng = np.array(bounds).T
3746            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3747            T0 = estimateT0(take_step)
3748            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3749            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3750                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3751                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3752            chisq = result.fun
3753            ncalc = result.nfev
3754            newVals = result.x
3755            covM = []
3756        elif data['Minimizer'] == 'MC/SA Anneal':
3757            xyrng = np.array(bounds).T
3758            result = G2mth.anneal(sumREFD, values, 
3759                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3760                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3761                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3762                ranRange=0.20,autoRan=False,dlg=None)
3763            newVals = result[0]
3764            parmDict.update(zip(varyList,newVals))
3765            chisq = result[1]
3766            ncalc = result[3]
3767            covM = []
3768            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3769        elif data['Minimizer'] == 'L-BFGS-B':
3770            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
3771                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3772            parmDict.update(zip(varyList,result['x']))
3773            chisq = result.fun
3774            ncalc = result.nfev
3775            newVals = result.x
3776            covM = []
3777    else:   #nothing varied
3778        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3779        chisq = np.sum(M**2)
3780        ncalc = 0
3781        covM = []
3782        sig = []
3783        sigDict = {}
3784        result = []
3785    Rvals = {}
3786    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
3787    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
3788    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3789    Ib[Ibeg:Ifin] = parmDict['FltBack']
3790    try:
3791        if not len(varyList):
3792            Msg += ' - nothing refined'
3793            raise ValueError
3794        Nans = np.isnan(newVals)
3795        if np.any(Nans):
3796            name = varyList[Nans.nonzero(True)[0]]
3797            Msg += ' Nan result for '+name+'!'
3798            raise ValueError
3799        Negs = np.less_equal(newVals,0.)
3800        if np.any(Negs):
3801            indx = Negs.nonzero()
3802            name = varyList[indx[0][0]]
3803            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3804                Msg += ' negative coefficient for '+name+'!'
3805                raise ValueError
3806        if len(covM):
3807            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3808            covMatrix = covM*Rvals['GOF']
3809        else:
3810            sig = np.zeros(len(varyList))
3811            covMatrix = []
3812        sigDict = dict(zip(varyList,sig))
3813        G2fil.G2Print (' Results of reflectometry data modelling fit:')
3814        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3815        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3816        SetModelParms()
3817        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3818    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
3819        G2fil.G2Print (Msg)
3820        return False,0,0,0,0,0,0,Msg
3821       
3822def makeSLDprofile(data,Substances):
3823   
3824    sq2 = np.sqrt(2.)
3825    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3826    Nlayers = len(laySeq)
3827    laySeq = np.array(laySeq,dtype=int)
3828    interfaces = np.zeros(Nlayers)
3829    rho = np.zeros(Nlayers)
3830    sigma = np.zeros(Nlayers)
3831    name = data['Layers'][0]['Name']
3832    thick = 0.
3833    for ilay,lay in enumerate(laySeq):
3834        layer = data['Layers'][lay]
3835        name = layer['Name']
3836        if 'Thick' in layer:
3837            thick += layer['Thick'][0]
3838            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3839        if 'Rough' in layer:
3840            sigma[ilay] = max(0.001,layer['Rough'][0])
3841        if name != 'vacuum':
3842            if name == 'unit scatter':
3843                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
3844            else:
3845                rrho = Substances[name]['Scatt density']
3846                irho = Substances[name]['XImag density']
3847                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
3848        if 'Mag SLD' in layer:
3849            rho[ilay] += layer['Mag SLD'][0]
3850    name = data['Layers'][-1]['Name']
3851    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
3852    xr = np.flipud(x)
3853    interfaces[-1] = x[-1]
3854    y = np.ones_like(x)*rho[0]
3855    iBeg = 0
3856    for ilayer in range(Nlayers-1):
3857        delt = rho[ilayer+1]-rho[ilayer]
3858        iPos = np.searchsorted(x,interfaces[ilayer])
3859        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
3860        iBeg = iPos
3861    return x,xr,y   
3862
3863def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3864   
3865    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3866    Qmin = Limits[1][0]
3867    Qmax = Limits[1][1]
3868    iBeg = np.searchsorted(Q,Qmin)
3869    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3870    Ib[:] = data['FltBack'][0]
3871    Ic[:] = 0
3872    Scale = data['Scale'][0]
3873    if data['Layer Seq'] == []:
3874        return
3875    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3876    Nlayers = len(laySeq)
3877    depth = np.zeros(Nlayers)
3878    rho = np.zeros(Nlayers)
3879    irho = np.zeros(Nlayers)
3880    sigma = np.zeros(Nlayers)
3881    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3882        layer = data['Layers'][lay]
3883        name = layer['Name']
3884        if 'Thick' in layer:    #skips first & last layers
3885            depth[ilay] = layer['Thick'][0]
3886        if 'Rough' in layer:    #skips first layer
3887            sigma[ilay] = layer['Rough'][0]
3888        if 'unit scatter' == name:
3889            rho[ilay] = layer['DenMul'][0]
3890            irho[ilay] = layer['iDenMul'][0]
3891        else:
3892            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
3893            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
3894        if 'Mag SLD' in layer:
3895            rho[ilay] += layer['Mag SLD'][0]
3896    if data['dQ type'] == 'None':
3897        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3898    elif 'const' in data['dQ type']:
3899        res = data['Resolution'][0]/(100.*sateln2)
3900        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
3901    else:       #dQ/Q in data
3902        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
3903    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
3904
3905def abeles(kz, depth, rho, irho=0, sigma=0):
3906    """
3907    Optical matrix form of the reflectivity calculation.
3908    O.S. Heavens, Optical Properties of Thin Solid Films
3909   
3910    Reflectometry as a function of kz for a set of slabs.
3911
3912    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
3913        This is :math:`\\tfrac12 Q_z`.       
3914    :param depth:  float[m] (Ang).
3915        thickness of each layer.  The thickness of the incident medium
3916        and substrate are ignored.
3917    :param rho:  float[n,k] (1e-6/Ang^2)
3918        Real scattering length density for each layer for each kz
3919    :param irho:  float[n,k] (1e-6/Ang^2)
3920        Imaginary scattering length density for each layer for each kz
3921        Note: absorption cross section mu = 2 irho/lambda for neutrons
3922    :param sigma: float[m-1] (Ang)
3923        interfacial roughness.  This is the roughness between a layer
3924        and the previous layer. The sigma array should have m-1 entries.
3925
3926    Slabs are ordered with the surface SLD at index 0 and substrate at
3927    index -1, or reversed if kz < 0.
3928    """
3929    def calc(kz, depth, rho, irho, sigma):
3930        if len(kz) == 0: return kz
3931   
3932        # Complex index of refraction is relative to the incident medium.
3933        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
3934        # in place of kz^2, and ignoring rho_o
3935        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
3936        k = kz
3937   
3938        # According to Heavens, the initial matrix should be [ 1 F; F 1],
3939        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
3940        # multiply versus some coding convenience.
3941        B11 = 1
3942        B22 = 1
3943        B21 = 0
3944        B12 = 0
3945        for i in range(0, len(depth)-1):
3946            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
3947            F = (k - k_next) / (k + k_next)
3948            F *= np.exp(-2*k*k_next*sigma[i]**2)
3949            #print "==== layer",i
3950            #print "kz:", kz
3951            #print "k:", k
3952            #print "k_next:",k_next
3953            #print "F:",F
3954            #print "rho:",rho[:,i+1]
3955            #print "irho:",irho[:,i+1]
3956            #print "d:",depth[i],"sigma:",sigma[i]
3957            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
3958            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
3959            M21 = F*M11
3960            M12 = F*M22
3961            C1 = B11*M11 + B21*M12
3962            C2 = B11*M21 + B21*M22
3963            B11 = C1
3964            B21 = C2
3965            C1 = B12*M11 + B22*M12
3966            C2 = B12*M21 + B22*M22
3967            B12 = C1
3968            B22 = C2
3969            k = k_next
3970   
3971        r = B12/B11
3972        return np.absolute(r)**2
3973
3974    if np.isscalar(kz): kz = np.asarray([kz], 'd')
3975
3976    m = len(depth)
3977
3978    # Make everything into arrays
3979    depth = np.asarray(depth,'d')
3980    rho = np.asarray(rho,'d')
3981    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
3982    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
3983
3984    # Repeat rho,irho columns as needed
3985    if len(rho.shape) == 1:
3986        rho = rho[None,:]
3987        irho = irho[None,:]
3988
3989    return calc(kz, depth, rho, irho, sigma)
3990   
3991def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
3992    y = abeles(kz, depth, rho, irho, sigma)
3993    s = dq/2.
3994    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
3995    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
3996    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
3997    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
3998    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
3999    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
4000    y *= 0.137023
4001    return y
4002       
4003def makeRefdFFT(Limits,Profile):
4004    G2fil.G2Print ('make fft')
4005    Q,Io = Profile[:2]
4006    Qmin = Limits[1][0]
4007    Qmax = Limits[1][1]
4008    iBeg = np.searchsorted(Q,Qmin)
4009    iFin = np.searchsorted(Q,Qmax)+1    #include last point
4010    Qf = np.linspace(0.,Q[iFin-1],5000)
4011    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
4012    If = QI(Qf)*Qf**4
4013    R = np.linspace(0.,5000.,5000)
4014    F = fft.rfft(If)
4015    return R,F
4016
4017   
4018#### Stacking fault simulation codes ################################################################################
4019def GetStackParms(Layers):
4020   
4021    Parms = []
4022#cell parms
4023    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
4024        Parms.append('cellA')
4025        Parms.append('cellC')
4026    else:
4027        Parms.append('cellA')
4028        Parms.append('cellB')
4029        Parms.append('cellC')
4030        if Layers['Laue'] != 'mmm':
4031            Parms.append('cellG')
4032#Transition parms
4033    for iY in range(len(Layers['Layers'])):
4034        for iX in range(len(Layers['Layers'])):
4035            Parms.append('TransP;%d;%d'%(iY,iX))
4036            Parms.append('TransX;%d;%d'%(iY,iX))
4037            Parms.append('TransY;%d;%d'%(iY,iX))
4038            Parms.append('TransZ;%d;%d'%(iY,iX))
4039    return Parms
4040
4041def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
4042    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
4043   
4044    :param dict Layers: dict with following items
4045
4046      ::
4047
4048       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
4049       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
4050       'Layers':[],'Stacking':[],'Transitions':[]}
4051       
4052    :param str ctrls: controls string to be written on DIFFaX controls.dif file
4053    :param float scale: scale factor
4054    :param dict background: background parameters
4055    :param list limits: min/max 2-theta to be calculated
4056    :param dict inst: instrument parameters dictionary
4057    :param list profile: powder pattern data
4058   
4059    Note that parameters all updated in place   
4060    '''
4061    import atmdata
4062    path = sys.path
4063    for name in path:
4064        if 'bin' in name:
4065            DIFFaX = name+'/DIFFaX.exe'
4066            G2fil.G2Print (' Execute '+DIFFaX)
4067            break
4068    # make form factor file that DIFFaX wants - atom types are GSASII style
4069    sf = open('data.sfc','w')
4070    sf.write('GSASII special form factor file for DIFFaX\n\n')
4071    atTypes = list(Layers['AtInfo'].keys())
4072    if 'H' not in atTypes:
4073        atTypes.insert(0,'H')
4074    for atType in atTypes:
4075        if atType == 'H': 
4076            blen = -.3741
4077        else:
4078            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
4079        Adat = atmdata.XrayFF[atType]
4080        text = '%4s'%(atType.ljust(4))
4081        for i in range(4):
4082            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
4083        text += '%11.6f%11.6f'%(Adat['fc'],blen)
4084        text += '%3d\n'%(Adat['Z'])
4085        sf.write(text)
4086    sf.close()
4087    #make DIFFaX control.dif file - future use GUI to set some of these flags
4088    cf = open('control.dif','w')
4089    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4090        x0 = profile[0]
4091        iBeg = np.searchsorted(x0,limits[0])
4092        iFin = np.searchsorted(x0,limits[1])+1
4093        if iFin-iBeg > 20000:
4094            iFin = iBeg+20000
4095        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
4096        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4097        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
4098    else:
4099        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4100        inst = {'Type':['XSC','XSC',]}
4101    cf.close()
4102    #make DIFFaX data file
4103    df = open('GSASII-DIFFaX.dat','w')
4104    df.write('INSTRUMENTAL\n')
4105    if 'X' in inst['Type'][0]:
4106        df.write('X-RAY\n')
4107    elif 'N' in inst['Type'][0]:
4108        df.write('NEUTRON\n')
4109    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4110        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
4111        U = ateln2*inst['U'][1]/10000.
4112        V = ateln2*inst['V'][1]/10000.
4113        W = ateln2*inst['W'][1]/10000.
4114        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
4115        HW = np.sqrt(np.mean(HWHM))
4116    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
4117        if 'Mean' in Layers['selInst']:
4118            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
4119        elif 'Gaussian' in Layers['selInst']:
4120            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
4121        else:
4122            df.write('None\n')
4123    else:
4124        df.write('0.10\nNone\n')
4125    df.write('STRUCTURAL\n')
4126    a,b,c = Layers['Cell'][1:4]
4127    gam = Layers['Cell'][6]
4128    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
4129    laue = Layers['Laue']
4130    if laue == '2/m(ab)':
4131        laue = '2/m(1)'
4132    elif laue == '2/m(c)':
4133        laue = '2/m(2)'
4134    if 'unknown' in Layers['Laue']:
4135        df.write('%s %.3f\n'%(laue,Layers['Toler']))
4136    else:   
4137        df.write('%s\n'%(laue))
4138    df.write('%d\n'%(len(Layers['Layers'])))
4139    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
4140        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
4141    layerNames = []
4142    for layer in Layers['Layers']:
4143        layerNames.append(layer['Name'])
4144    for il,layer in enumerate(Layers['Layers']):
4145        if layer['SameAs']:
4146            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
4147            continue
4148        df.write('LAYER %d\n'%(il+1))
4149        if '-1' in layer['Symm']:
4150            df.write('CENTROSYMMETRIC\n')
4151        else:
4152            df.write('NONE\n')
4153        for ia,atom in enumerate(layer['Atoms']):
4154            [name,atype,x,y,z,frac,Uiso] = atom
4155            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
4156                frac /= 2.
4157            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
4158    df.write('STACKING\n')
4159    df.write('%s\n'%(Layers['Stacking'][0]))
4160    if 'recursive' in Layers['Stacking'][0]:
4161        df.write('%s\n'%Layers['Stacking'][1])
4162    else:
4163        if 'list' in Layers['Stacking'][1]:
4164            Slen = len(Layers['Stacking'][2])
4165            iB = 0
4166            iF = 0
4167            while True:
4168                iF += 68
4169                if iF >= Slen:
4170                    break
4171                iF = min(iF,Slen)
4172                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
4173                iB = iF
4174        else:
4175            df.write('%s\n'%Layers['Stacking'][1])   
4176    df.write('TRANSITIONS\n')
4177    for iY in range(len(Layers['Layers'])):
4178        sumPx = 0.
4179        for iX in range(len(Layers['Layers'])):
4180            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
4181            p = round(p,3)
4182            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
4183            sumPx += p
4184        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
4185            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
4186            df.close()
4187            os.remove('data.sfc')
4188            os.remove('control.dif')
4189            os.remove('GSASII-DIFFaX.dat')
4190            return       
4191    df.close()   
4192    time0 = time.time()
4193    try:
4194        subp.call(DIFFaX)
4195    except OSError:
4196        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
4197    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
4198    if os.path.exists('GSASII-DIFFaX.spc'):
4199        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
4200        iFin = iBeg+Xpat.shape[1]
4201        bakType,backDict,backVary = SetBackgroundParms(background)
4202        backDict['Lam1'] = G2mth.getWave(inst)
4203        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
4204        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
4205        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
4206            rv = st.poisson(profile[3][iBeg:iFin])
4207            profile[1][iBeg:iFin] = rv.rvs()
4208            Z = np.ones_like(profile[3][iBeg:iFin])
4209            Z[1::2] *= -1
4210            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
4211            profile[2][iBeg:iFin] = np.where(profile[1][iBeg:iFin]>0.,1./profile[1][iBeg:iFin],1.0)
4212        profile[5][iBeg:iFin] = profile[1][iBeg:iFin]-profile[3][iBeg:iFin]
4213    #cleanup files..
4214        os.remove('GSASII-DIFFaX.spc')
4215    elif os.path.exists('GSASII-DIFFaX.sadp'):
4216        Sadp = np.fromfile('GSASII-DIFFaX.sadp','>u2')
4217        Sadp = np.reshape(Sadp,(256,-1))
4218        Layers['Sadp']['Img'] = Sadp
4219        os.remove('GSASII-DIFFaX.sadp')
4220    os.remove('data.sfc')
4221    os.remove('control.dif')
4222    os.remove('GSASII-DIFFaX.dat')
4223   
4224def SetPWDRscan(inst,limits,profile):
4225   
4226    wave = G2mth.getMeanWave(inst)
4227    x0 = profile[0]
4228    iBeg = np.searchsorted(x0,limits[0])
4229    iFin = np.searchsorted(x0,limits[1])
4230    if iFin-iBeg > 20000:
4231        iFin = iBeg+20000
4232    Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
4233    pyx.pygetinst(wave,x0[iBeg],x0[iFin],Dx)
4234    return iFin-iBeg
4235       
4236def SetStackingSF(Layers,debug):
4237# Load scattering factors into DIFFaX arrays
4238    import atmdata
4239    atTypes = Layers['AtInfo'].keys()
4240    aTypes = []
4241    for atype in atTypes:
4242        aTypes.append('%4s'%(atype.ljust(4)))
4243    SFdat = []
4244    for atType in atTypes:
4245        Adat = atmdata.XrayFF[atType]
4246        SF = np.zeros(9)
4247        SF[:8:2] = Adat['fa']
4248        SF[1:8:2] = Adat['fb']
4249        SF[8] = Adat['fc']
4250        SFdat.append(SF)
4251    SFdat = np.array(SFdat)
4252    pyx.pyloadscf(len(atTypes),aTypes,SFdat.T,debug)
4253   
4254def SetStackingClay(Layers,Type):
4255# Controls
4256    rand.seed()
4257    ranSeed = rand.randint(1,2**16-1)
4258    try:   
4259        laueId = ['-1','2/m(ab)','2/m(c)','mmm','-3','-3m','4/m','4/mmm',
4260            '6/m','6/mmm'].index(Layers['Laue'])+1
4261    except ValueError:  #for 'unknown'
4262        laueId = -1
4263    if 'SADP' in Type:
4264        planeId = ['h0l','0kl','hhl','h-hl'].index(Layers['Sadp']['Plane'])+1
4265        lmax = int(Layers['Sadp']['Lmax'])
4266    else:
4267        planeId = 0
4268        lmax = 0
4269# Sequences
4270    StkType = ['recursive','explicit'].index(Layers['Stacking'][0])
4271    try:
4272        StkParm = ['infinite','random','list'].index(Layers['Stacking'][1])
4273    except ValueError:
4274        StkParm = -1
4275    if StkParm == 2:    #list
4276        StkSeq = [int(val) for val in Layers['Stacking'][2].split()]
4277        Nstk = len(StkSeq)
4278    else:
4279        Nstk = 1
4280        StkSeq = [0,]
4281    if StkParm == -1:
4282        StkParm = int(Layers['Stacking'][1])
4283    Wdth = Layers['Width'][0]
4284    mult = 1
4285    controls = [laueId,planeId,lmax,mult,StkType,StkParm,ranSeed]
4286    LaueSym = Layers['Laue'].ljust(12)
4287    pyx.pygetclay(controls,LaueSym,Wdth,Nstk,StkSeq)
4288    return laueId,controls
4289   
4290def SetCellAtoms(Layers):
4291    Cell = Layers['Cell'][1:4]+Layers['Cell'][6:7]
4292# atoms in layers
4293    atTypes = list(Layers['AtInfo'].keys())
4294    AtomXOU = []
4295    AtomTp = []
4296    LayerSymm = []
4297    LayerNum = []
4298    layerNames = []
4299    Natm = 0
4300    Nuniq = 0
4301    for layer in Layers['Layers']:
4302        layerNames.append(layer['Name'])
4303    for il,layer in enumerate(Layers['Layers']):
4304        if layer['SameAs']:
4305            LayerNum.append(layerNames.index(layer['SameAs'])+1)
4306            continue
4307        else:
4308            LayerNum.append(il+1)
4309            Nuniq += 1
4310        if '-1' in layer['Symm']:
4311            LayerSymm.append(1)
4312        else:
4313            LayerSymm.append(0)
4314        for ia,atom in enumerate(layer['Atoms']):
4315            [name,atype,x,y,z,frac,Uiso] = atom
4316            Natm += 1
4317            AtomTp.append('%4s'%(atype.ljust(4)))
4318            Ta = atTypes.index(atype)+1
4319            AtomXOU.append([float(Nuniq),float(ia+1),float(Ta),x,y,z,frac,Uiso*78.9568])
4320    AtomXOU = np.array(AtomXOU)
4321    Nlayers = len(layerNames)
4322    pyx.pycellayer(Cell,Natm,AtomTp,AtomXOU.T,Nuniq,LayerSymm,Nlayers,LayerNum)
4323    return Nlayers
4324   
4325def SetStackingTrans(Layers,Nlayers):
4326# Transitions
4327    TransX = []
4328    TransP = []
4329    for Ytrans in Layers['Transitions']:
4330        TransP.append([trans[0] for trans in Ytrans])   #get just the numbers
4331        TransX.append([trans[1:4] for trans in Ytrans])   #get just the numbers
4332    TransP = np.array(TransP,dtype='float').T
4333    TransX = np.array(TransX,dtype='float')
4334#    GSASIIpath.IPyBreak()
4335    pyx.pygettrans(Nlayers,TransP,TransX)
4336   
4337def CalcStackingPWDR