source: trunk/GSASIIpwd.py @ 5161

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

cleanup constraint display - remove excessive wx.RIGHT,wx.LEFT, etc. - now cleaner display
fix to ShowScrolledColText? - should be txtSizer.AddGrowableCol?(0); (1) caused crash
fix issues with seq PDFfit not copy forward, etc.

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