source: trunk/GSASIIpwd.py @ 5152

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

fix RMC GUI options for PDFfit
fix background Bragg peak fits for EDX data

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 197.3 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-21 20:55:16 +0000 (Fri, 21 Jan 2022) $
10# $Author: vondreele $
11# $Revision: 5152 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 5152 2022-01-21 20:55:16Z 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: 5152 $"
38GSASIIpath.SetVersionNumber("$Revision: 5152 $")
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 = 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['@%s'%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['@%s'%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(%s,%.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 UpdatePDFfit(Phase,RMCPdict):
3291    ''' Updates various PDFfit parameters held in GSAS-II
3292    '''
3293   
3294    General = Phase['General']
3295    if RMCPdict['refinement'] == 'normal':
3296        fName = General['Name']+'-PDFfit.rstr'
3297        try:
3298            rstr = open(fName.replace(' ','_'),'r')
3299        except FileNotFoundError:
3300            return [fName,'Not found - PDFfit failed']
3301        lines = rstr.readlines()
3302        rstr.close()
3303        header = [line[:-1].split(' ',1) for line in lines[:7]]
3304        resdict = dict(header)
3305        for item in ['scale','sharp','cell']:
3306            resdict[item] = [float(val) for val in resdict[item].split(',')]
3307        General['Cell'][1:7] = resdict['cell']
3308        for inam,name in enumerate(['delta2','delta1','sratio']):
3309            RMCPdict[name][0] = resdict['sharp'][inam]
3310        if 'shape' in resdict:
3311            if 'sphere' in resdict['shape']:
3312                RMCPdict['spdiameter'][0] = resdict['shape'][-1]
3313            else:
3314                RMCPdict['stepcut'][0] = resdict['shape'][-1]
3315        cx,ct,cs,ci = G2mth.getAtomPtrs(Phase)     
3316        Atoms = Phase['Atoms']
3317        atmBeg = 0
3318        for line in lines:
3319            atmBeg += 1
3320            if 'atoms' in line:
3321                break
3322        for atom in Atoms:
3323            atstr = lines[atmBeg][:-1].split()
3324            Uiistr = lines[atmBeg+2][:-1].split()
3325            Uijstr = lines[atmBeg+4][:-1].split()
3326            atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])]
3327            atom[ci] = 'A'
3328            atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])]
3329            atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])]
3330            atmBeg += 6
3331           
3332    fName = 'Sequential_PDFfit.res'
3333    if RMCPdict['refinement'] == 'normal':
3334        fName = General['Name']+'-PDFfit.res'
3335    try:
3336        res = open(fName.replace(' ','_'),'r')
3337    except FileNotFoundError:
3338        return [fName,'Not found - PDFfit failed']
3339    lines = res.readlines()
3340    res.close()
3341    Ibeg = False
3342    resline = ''
3343    for iline,line in enumerate(lines):
3344        if 'Refinement parameters' in line:
3345            Ibeg = True
3346            continue
3347        if Ibeg:
3348            if '---------' in line:
3349                break
3350            resline += line[:-1]
3351    for iline,line in enumerate(lines):
3352        if 'Rw - ' in line:
3353            if 'nan' in line:
3354                Rwp = 100.0
3355            else:
3356                Rwp = float(line.split(':')[1])
3357    results = resline.replace('(','').split(')')[:-1]
3358    results = ['@'+result.lstrip() for result in results]
3359    results = [item.split() for item in results]
3360    if RMCPdict['refinement'] == 'normal':
3361        results = dict([[item[0][:-1],float(item[1])] for item in results if item[0][:-1] in RMCPdict['AtomVar']])
3362        RMCPdict['AtomVar'].update(results)
3363        return None
3364    else:   #sequential
3365        newParms = dict([[item[0][1:-1],[float(item[1]),float(item[2])]] for item in results])
3366        return newParms,Rwp
3367       
3368def MakefullrmcRun(pName,Phase,RMCPdict):
3369    '''Creates a script to run fullrmc. Returns the name of the file that was
3370    created.
3371    '''
3372    BondList = {}
3373    for k in RMCPdict['Pairs']:
3374        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
3375            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
3376    AngleList = []
3377    for angle in RMCPdict['Angles']:
3378        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
3379            continue
3380        for i in (0,1,2):
3381            angle[i] = angle[i].strip()
3382        AngleList.append(angle)
3383    # rmin = RMCPdict['min Contact']
3384    cell = Phase['General']['Cell'][1:7]
3385    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
3386    cx,ct,cs,cia = Phase['General']['AtomPtrs']
3387    atomsList = []
3388    for atom in Phase['Atoms']:
3389        el = ''.join([i for i in atom[ct] if i.isalpha()])
3390        atomsList.append([el] + atom[cx:cx+4])
3391    projDir,projName = os.path.split(os.path.abspath(pName))
3392    scrname = pName+'-fullrmc.py'
3393    restart = '%s_restart.pdb'%pName
3394    Files = RMCPdict['files']
3395    rundata = ''
3396    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
3397    rundata += '# created in '+__file__+" v"+filversion.split()[1]
3398    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
3399    rundata += '''
3400# fullrmc imports (all that are potentially useful)
3401import os,glob
3402import time
3403import pickle
3404import numpy as np
3405from fullrmc.Core import Collection
3406from fullrmc.Engine import Engine
3407import fullrmc.Constraints.PairDistributionConstraints as fPDF
3408from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
3409from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
3410from fullrmc.Constraints.BondConstraints import BondConstraint
3411from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
3412from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
3413from fullrmc.Generators.Swaps import SwapPositionsGenerator
3414# utility routines
3415def writeHeader(ENGINE,statFP):
3416    'header for stats file'
3417    statFP.write('generated-steps, total-error, ')
3418    for c in ENGINE.constraints:
3419        statFP.write(c.constraintName)
3420        statFP.write(', ')
3421    statFP.write('\\n')
3422    statFP.flush()
3423   
3424def writeCurrentStatus(ENGINE,statFP,plotF):
3425    'line in stats file & current constraint plots'
3426    statFP.write(str(ENGINE.generated))
3427    statFP.write(', ')
3428    statFP.write(str(ENGINE.totalStandardError))
3429    statFP.write(', ')
3430    for c in ENGINE.constraints:
3431        statFP.write(str(c.standardError))
3432        statFP.write(', ')
3433    statFP.write('\\n')
3434    statFP.flush()
3435    mpl.use('agg')
3436    fp = open(plotF,'wb')
3437    for c in ENGINE.constraints:
3438        p = c.plot(show=False)
3439        p[0].canvas.draw()
3440        image = p[0].canvas.buffer_rgba()
3441        pickle.dump(c.constraintName,fp)
3442        pickle.dump(np.array(image),fp)
3443    fp.close()
3444
3445def calcRmax(ENGINE):
3446    'from Bachir, works for non-othorhombic'
3447    a,b,c = ENGINE.basisVectors
3448    lens = []
3449    ts    = np.linalg.norm(np.cross(a,b))/2
3450    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] )
3451    ts = np.linalg.norm(np.cross(b,c))/2
3452    lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] )
3453    ts = np.linalg.norm(np.cross(a,c))/2
3454    lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] )
3455    return min(lens)
3456'''
3457    rundata += '''
3458### When True, erases an existing engine to provide a fresh start
3459FRESH_START = {:}
3460dirName = "{:}"
3461prefix = "{:}"
3462project = prefix + "-fullrmc"
3463time0 = time.time()
3464'''.format(RMCPdict['ReStart'][0],projDir,projName)
3465   
3466    rundata += '# setup structure\n'
3467    rundata += 'cell = ' + str(cell) + '\n'
3468    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
3469    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
3470    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
3471
3472    rundata += '\n# initialize engine\n'
3473    rundata += '''
3474engineFileName = os.path.join(dirName, project + '.rmc')
3475projectStats = os.path.join(dirName, project + '.stats')
3476projectPlots = os.path.join(dirName, project + '.plots')
3477pdbFile = os.path.join(dirName, project + '_restart.pdb')
3478# check Engine exists if so (and not FRESH_START) load it
3479# otherwise build it
3480ENGINE = Engine(path=None)
3481if not ENGINE.is_engine(engineFileName) or FRESH_START:
3482    ## create structure
3483    ENGINE = Engine(path=engineFileName, freshStart=True)
3484    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
3485                                 atoms      = atomList,
3486                                 unitcellBC = cell,
3487                                 supercell  = supercell)
3488'''   
3489    import atmdata
3490    # rundata += '    # conversion factors (may be needed)\n'
3491    # rundata += '    sumCiBi2 = 0.\n'
3492    # for elem in Phase['General']['AtomTypes']:
3493    #     rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
3494    #     rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
3495    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
3496    # settings that require a new Engine
3497    for File in Files:
3498        filDat = RMCPdict['files'][File]
3499        if not os.path.exists(filDat[0]): continue
3500        sfwt = 'neutronCohb'
3501        if 'Xray' in File:
3502            sfwt = 'atomicNumber'
3503        if 'G(r)' in File:
3504            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3505            if filDat[3] == 0:
3506                #rundata += '''    # read and xform G(r) as defined in RMCProfile
3507    # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
3508                #rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
3509                #rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3510                rundata += '    # G(r) as defined in RMCProfile\n'
3511                rundata += '    GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3512            elif filDat[3] == 1:
3513                rundata += '    # This is G(r) as defined in PDFFIT\n'
3514                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3515            elif filDat[3] == 2:
3516                rundata += '    # This is g(r)\n'
3517                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
3518            else:
3519                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
3520            rundata += '    ENGINE.add_constraints([GofR])\n'
3521            rundata += '    GofR.set_limits((None, calcRmax(ENGINE)))\n'
3522        elif '(Q)' in File:
3523            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
3524            if filDat[3] == 0:
3525                rundata += '    # F(Q) as defined in RMCProfile\n'
3526                #rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
3527                if filDat[4]:
3528                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3529                rundata += '    SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3530            elif filDat[3] == 1:
3531                rundata += '    # S(Q) as defined in PDFFIT\n'
3532                rundata += '    SOQ[1] -= 1\n'
3533                if filDat[4]:
3534                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n'
3535                rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
3536            else:
3537                raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
3538            rundata += '    ENGINE.add_constraints([SofQ])\n'
3539        else:
3540            print('What is this?')
3541    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
3542    if BondList:
3543        rundata += '''    B_CONSTRAINT   = BondConstraint()
3544    ENGINE.add_constraints(B_CONSTRAINT)
3545    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
3546'''
3547        for pair in BondList:
3548            e1,e2 = pair.split('-')
3549            d1,d2 = BondList[pair]
3550            if d1 == 0: continue
3551            if d2 == 0:
3552                print('Ignoring min distance without maximum')
3553                #rundata += '            ("element","{}","{}",{}),\n'.format(
3554                #                        e1.strip(),e2.strip(),d1)
3555            else:
3556                rundata += '            ("element","{}","{}",{},{}),\n'.format(
3557                                        e1.strip(),e2.strip(),d1,d2)
3558        rundata += '             ])\n'
3559    if AngleList:
3560        rundata += '''    A_CONSTRAINT   = BondsAngleConstraint()
3561    ENGINE.add_constraints(A_CONSTRAINT)
3562    A_CONSTRAINT.create_supercell_angles(anglesDefinition=[
3563'''
3564        for item in AngleList:
3565            rundata += ('            '+
3566               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
3567        rundata += '             ])\n'
3568    rundata += '''
3569    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
3570    ENGINE.save()
3571else:
3572    ENGINE = ENGINE.load(path=engineFileName)
3573
3574ENGINE.set_log_file(os.path.join(dirName,prefix))
3575'''
3576    if RMCPdict['Swaps']:
3577        rundata += '\n#set up for site swaps\n'
3578        rundata += 'aN = ENGINE.allNames\n'
3579        rundata += 'SwapGen = {}\n'
3580        for swap in RMCPdict['Swaps']:
3581            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
3582            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
3583            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
3584        rundata += '    for swaps in SwapGen:\n'
3585        rundata += '        AB = swaps.split("-")\n'
3586        rundata += '        ENGINE.set_groups_as_atoms()\n'
3587        rundata += '        for g in ENGINE.groups:\n'
3588        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
3589        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
3590        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
3591        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
3592        rundata += '            sProb = SwapGen[swaps][2]\n'
3593    rundata += '\n# set weights -- do this now so values can be changed without a restart\n'
3594    # rundata += 'wtDict = {}\n'
3595    # for File in Files:
3596    #     filDat = RMCPdict['files'][File]
3597    #     if not os.path.exists(filDat[0]): continue
3598    #     if 'Xray' in File:
3599    #         sfwt = 'atomicNumber'
3600    #     else:
3601    #         sfwt = 'neutronCohb'
3602    #     if 'G(r)' in File:
3603    #         typ = 'Pair'
3604    #     elif '(Q)' in File:
3605    #         typ = 'Struct'
3606    #     rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1])
3607    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
3608    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
3609    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
3610    rundata += '        c.set_limits((None,calcRmax(ENGINE)))\n'
3611    if RMCPdict['FitScale']:
3612        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3613    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
3614    if RMCPdict['FitScale']:
3615        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
3616        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
3617    # torsions difficult to implement, must be internal to cell & named with
3618    # fullrmc atom names
3619    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
3620    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
3621    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
3622    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
3623    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
3624    #     for torsion in RMCPdict['Torsions']:
3625    #         rundata += '    %s\n'%str(tuple(torsion))
3626    #     rundata += '        ]})\n'           
3627    rundata += '''
3628if FRESH_START:
3629    # initialize engine with one step to get starting config energetics
3630    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
3631    statFP = open(projectStats,'w')
3632    writeHeader(ENGINE,statFP)
3633    writeCurrentStatus(ENGINE,statFP,projectPlots)
3634else:
3635    statFP = open(projectStats,'a')
3636
3637# setup runs for fullrmc
3638'''
3639    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
3640    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
3641    rundata += '    ENGINE.set_groups_as_atoms()\n'
3642    rundata += '    expected = ENGINE.generated+steps\n'
3643   
3644    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
3645    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
3646    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
3647    rundata += 'statFP.close()\n'
3648    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
3649    rfile = open(scrname,'w')
3650    rfile.writelines(rundata)
3651    rfile.close()
3652    return scrname
3653   
3654def GetRMCBonds(general,RMCPdict,Atoms,bondList):
3655    bondDist = []
3656    Cell = general['Cell'][1:7]
3657    Supercell =  RMCPdict['SuperCell']
3658    Trans = np.eye(3)*np.array(Supercell)
3659    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3660    Amat,Bmat = G2lat.cell2AB(Cell)
3661    indices = (-1,0,1)
3662    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3663    for bonds in bondList:
3664        Oxyz = np.array(Atoms[bonds[0]][1:])
3665        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
3666        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
3667        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
3668        for dx in Dx.T:
3669            bondDist.append(np.min(dx))
3670    return np.array(bondDist)
3671   
3672def GetRMCAngles(general,RMCPdict,Atoms,angleList):
3673    bondAngles = []
3674    Cell = general['Cell'][1:7]
3675    Supercell =  RMCPdict['SuperCell']
3676    Trans = np.eye(3)*np.array(Supercell)
3677    Cell = G2lat.TransformCell(Cell,Trans)[:6]
3678    Amat,Bmat = G2lat.cell2AB(Cell)
3679    indices = (-1,0,1)
3680    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3681    for angle in angleList:
3682        Oxyz = np.array(Atoms[angle[0]][1:])
3683        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
3684        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])
3685        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
3686        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
3687        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
3688        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
3689        iDAx = np.argmin(DAx,axis=0)
3690        iDBx = np.argmin(DBx,axis=0)
3691        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
3692            DAv = DAxV[iA,i]/DAx[iA,i]
3693            DBv = DBxV[iB,i]/DBx[iB,i]
3694            bondAngles.append(npacosd(np.sum(DAv*DBv)))
3695    return np.array(bondAngles)
3696   
3697#### Reflectometry calculations ################################################################################
3698def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
3699    G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0]))
3700   
3701    class RandomDisplacementBounds(object):
3702        """random displacement with bounds"""
3703        def __init__(self, xmin, xmax, stepsize=0.5):
3704            self.xmin = xmin
3705            self.xmax = xmax
3706            self.stepsize = stepsize
3707   
3708        def __call__(self, x):
3709            """take a random step but ensure the new position is within the bounds"""
3710            while True:
3711                # this could be done in a much more clever way, but it will work for example purposes
3712                steps = self.xmax-self.xmin
3713                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
3714                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
3715                    break
3716            return xnew
3717   
3718    def GetModelParms():
3719        parmDict = {}
3720        varyList = []
3721        values = []
3722        bounds = []
3723        parmDict['dQ type'] = data['dQ type']
3724        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
3725        for parm in ['Scale','FltBack']:
3726            parmDict[parm] = data[parm][0]
3727            if data[parm][1]:
3728                varyList.append(parm)
3729                values.append(data[parm][0])
3730                bounds.append(Bounds[parm])
3731        parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int)
3732        parmDict['nLayers'] = len(parmDict['Layer Seq'])
3733        for ilay,layer in enumerate(data['Layers']):
3734            name = layer['Name']
3735            cid = str(ilay)+';'
3736            parmDict[cid+'Name'] = name
3737            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3738                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
3739                if layer.get(parm,[0.,False])[1]:
3740                    varyList.append(cid+parm)
3741                    value = layer[parm][0]
3742                    values.append(value)
3743                    if value:
3744                        bound = [value*Bfac,value/Bfac]
3745                    else:
3746                        bound = [0.,10.]
3747                    bounds.append(bound)
3748            if name not in ['vacuum','unit scatter']:
3749                parmDict[cid+'rho'] = Substances[name]['Scatt density']
3750                parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
3751        return parmDict,varyList,values,bounds
3752   
3753    def SetModelParms():
3754        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
3755        if 'Scale' in varyList:
3756            data['Scale'][0] = parmDict['Scale']
3757            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
3758        G2fil.G2Print (line)
3759        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
3760        if 'FltBack' in varyList:
3761            data['FltBack'][0] = parmDict['FltBack']
3762            line += ' esd: %15.3g'%(sigDict['FltBack'])
3763        G2fil.G2Print (line)
3764        for ilay,layer in enumerate(data['Layers']):
3765            name = layer['Name']
3766            G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name))
3767            cid = str(ilay)+';'
3768            line = ' '
3769            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
3770            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul'])
3771            for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']:
3772                if parm in layer:
3773                    if parm == 'Rough':
3774                        layer[parm][0] = abs(parmDict[cid+parm])    #make positive
3775                    else:
3776                        layer[parm][0] = parmDict[cid+parm]
3777                    line += ' %s: %.3f'%(parm,layer[parm][0])
3778                    if cid+parm in varyList:
3779                        line += ' esd: %.3g'%(sigDict[cid+parm])
3780            G2fil.G2Print (line)
3781            G2fil.G2Print (line2)
3782   
3783    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3784        parmDict.update(zip(varyList,values))
3785        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3786        return M
3787   
3788    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
3789        parmDict.update(zip(varyList,values))
3790        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
3791        return np.sum(M**2)
3792   
3793    def getREFD(Q,Qsig,parmDict):
3794        Ic = np.ones_like(Q)*parmDict['FltBack']
3795        Scale = parmDict['Scale']
3796        Nlayers = parmDict['nLayers']
3797        Res = parmDict['Res']
3798        depth = np.zeros(Nlayers)
3799        rho = np.zeros(Nlayers)
3800        irho = np.zeros(Nlayers)
3801        sigma = np.zeros(Nlayers)
3802        for ilay,lay in enumerate(parmDict['Layer Seq']):
3803            cid = str(lay)+';'
3804            depth[ilay] = parmDict[cid+'Thick']
3805            sigma[ilay] = parmDict[cid+'Rough']
3806            if parmDict[cid+'Name'] == u'unit scatter':
3807                rho[ilay] = parmDict[cid+'DenMul']
3808                irho[ilay] = parmDict[cid+'iDenMul']
3809            elif 'vacuum' != parmDict[cid+'Name']:
3810                rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
3811                irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
3812            if cid+'Mag SLD' in parmDict:
3813                rho[ilay] += parmDict[cid+'Mag SLD']
3814        if parmDict['dQ type'] == 'None':
3815            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
3816        elif 'const' in parmDict['dQ type']:
3817            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
3818        else:       #dQ/Q in data
3819            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
3820        Ic += AB*Scale
3821        return Ic
3822       
3823    def estimateT0(takestep):
3824        Mmax = -1.e-10
3825        Mmin = 1.e10
3826        for i in range(100):
3827            x0 = takestep(values)
3828            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3829            Mmin = min(M,Mmin)
3830            MMax = max(M,Mmax)
3831        return 1.5*(MMax-Mmin)
3832
3833    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3834    if data.get('2% weight'):
3835        wt = 1./(0.02*Io)**2
3836    Qmin = Limits[1][0]
3837    Qmax = Limits[1][1]
3838    wtFactor = ProfDict['wtFactor']
3839    Bfac = data['Toler']
3840    Ibeg = np.searchsorted(Q,Qmin)
3841    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
3842    Ic[:] = 0
3843    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
3844              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
3845    parmDict,varyList,values,bounds = GetModelParms()
3846    Msg = 'Failed to converge'
3847    if varyList:
3848        if data['Minimizer'] == 'LMLS': 
3849            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6,
3850                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3851            parmDict.update(zip(varyList,result[0]))
3852            chisq = np.sum(result[2]['fvec']**2)
3853            ncalc = result[2]['nfev']
3854            covM = result[1]
3855            newVals = result[0]
3856        elif data['Minimizer'] == 'Basin Hopping':
3857            xyrng = np.array(bounds).T
3858            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
3859            T0 = estimateT0(take_step)
3860            G2fil.G2Print (' Estimated temperature: %.3g'%(T0))
3861            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
3862                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
3863                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
3864            chisq = result.fun
3865            ncalc = result.nfev
3866            newVals = result.x
3867            covM = []
3868        elif data['Minimizer'] == 'MC/SA Anneal':
3869            xyrng = np.array(bounds).T
3870            result = G2mth.anneal(sumREFD, values, 
3871                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
3872                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
3873                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
3874                ranRange=0.20,autoRan=False,dlg=None)
3875            newVals = result[0]
3876            parmDict.update(zip(varyList,newVals))
3877            chisq = result[1]
3878            ncalc = result[3]
3879            covM = []
3880            G2fil.G2Print (' MC/SA final temperature: %.4g'%(result[2]))
3881        elif data['Minimizer'] == 'L-BFGS-B':
3882            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
3883                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
3884            parmDict.update(zip(varyList,result['x']))
3885            chisq = result.fun
3886            ncalc = result.nfev
3887            newVals = result.x
3888            covM = []
3889    else:   #nothing varied
3890        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
3891        chisq = np.sum(M**2)
3892        ncalc = 0
3893        covM = []
3894        sig = []
3895        sigDict = {}
3896        result = []
3897    Rvals = {}
3898    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
3899    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
3900    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
3901    Ib[Ibeg:Ifin] = parmDict['FltBack']
3902    try:
3903        if not len(varyList):
3904            Msg += ' - nothing refined'
3905            raise ValueError
3906        Nans = np.isnan(newVals)
3907        if np.any(Nans):
3908            name = varyList[Nans.nonzero(True)[0]]
3909            Msg += ' Nan result for '+name+'!'
3910            raise ValueError
3911        Negs = np.less_equal(newVals,0.)
3912        if np.any(Negs):
3913            indx = Negs.nonzero()
3914            name = varyList[indx[0][0]]
3915            if name != 'FltBack' and name.split(';')[1] in ['Thick',]:
3916                Msg += ' negative coefficient for '+name+'!'
3917                raise ValueError
3918        if len(covM):
3919            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
3920            covMatrix = covM*Rvals['GOF']
3921        else:
3922            sig = np.zeros(len(varyList))
3923            covMatrix = []
3924        sigDict = dict(zip(varyList,sig))
3925        G2fil.G2Print (' Results of reflectometry data modelling fit:')
3926        G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
3927        G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
3928        SetModelParms()
3929        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
3930    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
3931        G2fil.G2Print (Msg)
3932        return False,0,0,0,0,0,0,Msg
3933       
3934def makeSLDprofile(data,Substances):
3935   
3936    sq2 = np.sqrt(2.)
3937    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3938    Nlayers = len(laySeq)
3939    laySeq = np.array(laySeq,dtype=int)
3940    interfaces = np.zeros(Nlayers)
3941    rho = np.zeros(Nlayers)
3942    sigma = np.zeros(Nlayers)
3943    name = data['Layers'][0]['Name']
3944    thick = 0.
3945    for ilay,lay in enumerate(laySeq):
3946        layer = data['Layers'][lay]
3947        name = layer['Name']
3948        if 'Thick' in layer:
3949            thick += layer['Thick'][0]
3950            interfaces[ilay] = layer['Thick'][0]+interfaces[ilay-1]
3951        if 'Rough' in layer:
3952            sigma[ilay] = max(0.001,layer['Rough'][0])
3953        if name != 'vacuum':
3954            if name == 'unit scatter':
3955                rho[ilay] = np.sqrt(layer['DenMul'][0]**2+layer['iDenMul'][0]**2)
3956            else:
3957                rrho = Substances[name]['Scatt density']
3958                irho = Substances[name]['XImag density']
3959                rho[ilay] = np.sqrt(rrho**2+irho**2)*layer['DenMul'][0]
3960        if 'Mag SLD' in layer:
3961            rho[ilay] += layer['Mag SLD'][0]
3962    name = data['Layers'][-1]['Name']
3963    x = np.linspace(-0.15*thick,1.15*thick,1000,endpoint=True)
3964    xr = np.flipud(x)
3965    interfaces[-1] = x[-1]
3966    y = np.ones_like(x)*rho[0]
3967    iBeg = 0
3968    for ilayer in range(Nlayers-1):
3969        delt = rho[ilayer+1]-rho[ilayer]
3970        iPos = np.searchsorted(x,interfaces[ilayer])
3971        y[iBeg:] += (delt/2.)*sp.erfc((interfaces[ilayer]-x[iBeg:])/(sq2*sigma[ilayer+1]))
3972        iBeg = iPos
3973    return x,xr,y   
3974
3975def REFDModelFxn(Profile,Inst,Limits,Substances,data):
3976   
3977    Q,Io,wt,Ic,Ib,Qsig = Profile[:6]
3978    Qmin = Limits[1][0]
3979    Qmax = Limits[1][1]
3980    iBeg = np.searchsorted(Q,Qmin)
3981    iFin = np.searchsorted(Q,Qmax)+1    #include last point
3982    Ib[:] = data['FltBack'][0]
3983    Ic[:] = 0
3984    Scale = data['Scale'][0]
3985    if data['Layer Seq'] == []:
3986        return
3987    laySeq = ['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),]
3988    Nlayers = len(laySeq)
3989    depth = np.zeros(Nlayers)
3990    rho = np.zeros(Nlayers)
3991    irho = np.zeros(Nlayers)
3992    sigma = np.zeros(Nlayers)
3993    for ilay,lay in enumerate(np.array(laySeq,dtype=int)):
3994        layer = data['Layers'][lay]
3995        name = layer['Name']
3996        if 'Thick' in layer:    #skips first & last layers
3997            depth[ilay] = layer['Thick'][0]
3998        if 'Rough' in layer:    #skips first layer
3999            sigma[ilay] = layer['Rough'][0]
4000        if 'unit scatter' == name:
4001            rho[ilay] = layer['DenMul'][0]
4002            irho[ilay] = layer['iDenMul'][0]
4003        else:
4004            rho[ilay] = Substances[name]['Scatt density']*layer['DenMul'][0]
4005            irho[ilay] = Substances[name].get('XImag density',0.)*layer['DenMul'][0]
4006        if 'Mag SLD' in layer:
4007            rho[ilay] += layer['Mag SLD'][0]
4008    if data['dQ type'] == 'None':
4009        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
4010    elif 'const' in data['dQ type']:
4011        res = data['Resolution'][0]/(100.*sateln2)
4012        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
4013    else:       #dQ/Q in data
4014        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
4015    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
4016
4017def abeles(kz, depth, rho, irho=0, sigma=0):
4018    """
4019    Optical matrix form of the reflectivity calculation.
4020    O.S. Heavens, Optical Properties of Thin Solid Films
4021   
4022    Reflectometry as a function of kz for a set of slabs.
4023
4024    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
4025        This is :math:`\\tfrac12 Q_z`.       
4026    :param depth:  float[m] (Ang).
4027        thickness of each layer.  The thickness of the incident medium
4028        and substrate are ignored.
4029    :param rho:  float[n,k] (1e-6/Ang^2)
4030        Real scattering length density for each layer for each kz
4031    :param irho:  float[n,k] (1e-6/Ang^2)
4032        Imaginary scattering length density for each layer for each kz
4033        Note: absorption cross section mu = 2 irho/lambda for neutrons
4034    :param sigma: float[m-1] (Ang)
4035        interfacial roughness.  This is the roughness between a layer
4036        and the previous layer. The sigma array should have m-1 entries.
4037
4038    Slabs are ordered with the surface SLD at index 0 and substrate at
4039    index -1, or reversed if kz < 0.
4040    """
4041    def calc(kz, depth, rho, irho, sigma):
4042        if len(kz) == 0: return kz
4043   
4044        # Complex index of refraction is relative to the incident medium.
4045        # We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o
4046        # in place of kz^2, and ignoring rho_o
4047        kz_sq = kz**2 + 4e-6*np.pi*rho[:,0]
4048        k = kz
4049   
4050        # According to Heavens, the initial matrix should be [ 1 F; F 1],
4051        # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix
4052        # multiply versus some coding convenience.
4053        B11 = 1
4054        B22 = 1
4055        B21 = 0
4056        B12 = 0
4057        for i in range(0, len(depth)-1):
4058            k_next = np.sqrt(kz_sq - 4e-6*np.pi*(rho[:,i+1] + 1j*irho[:,i+1]))
4059            F = (k - k_next) / (k + k_next)
4060            F *= np.exp(-2*k*k_next*sigma[i]**2)
4061            #print "==== layer",i
4062            #print "kz:", kz
4063            #print "k:", k
4064            #print "k_next:",k_next
4065            #print "F:",F
4066            #print "rho:",rho[:,i+1]
4067            #print "irho:",irho[:,i+1]
4068            #print "d:",depth[i],"sigma:",sigma[i]
4069            M11 = np.exp(1j*k*depth[i]) if i>0 else 1
4070            M22 = np.exp(-1j*k*depth[i]) if i>0 else 1
4071            M21 = F*M11
4072            M12 = F*M22
4073            C1 = B11*M11 + B21*M12
4074            C2 = B11*M21 + B21*M22
4075            B11 = C1
4076            B21 = C2
4077            C1 = B12*M11 + B22*M12
4078            C2 = B12*M21 + B22*M22
4079            B12 = C1
4080            B22 = C2
4081            k = k_next
4082   
4083        r = B12/B11
4084        return np.absolute(r)**2
4085
4086    if np.isscalar(kz): kz = np.asarray([kz], 'd')
4087
4088    m = len(depth)
4089
4090    # Make everything into arrays
4091    depth = np.asarray(depth,'d')
4092    rho = np.asarray(rho,'d')
4093    irho = irho*np.ones_like(rho) if np.isscalar(irho) else np.asarray(irho,'d')
4094    sigma = sigma*np.ones(m-1,'d') if np.isscalar(sigma) else np.asarray(sigma,'d')
4095
4096    # Repeat rho,irho columns as needed
4097    if len(rho.shape) == 1:
4098        rho = rho[None,:]
4099        irho = irho[None,:]
4100
4101    return calc(kz, depth, rho, irho, sigma)
4102   
4103def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
4104    y = abeles(kz, depth, rho, irho, sigma)
4105    s = dq/2.
4106    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
4107    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma)) 
4108    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma)) 
4109    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
4110    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
4111    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
4112    y *= 0.137023
4113    return y
4114       
4115def makeRefdFFT(Limits,Profile):
4116    G2fil.G2Print ('make fft')
4117    Q,Io = Profile[:2]
4118    Qmin = Limits[1][0]
4119    Qmax = Limits[1][1]
4120    iBeg = np.searchsorted(Q,Qmin)
4121    iFin = np.searchsorted(Q,Qmax)+1    #include last point
4122    Qf = np.linspace(0.,Q[iFin-1],5000)
4123    QI = si.interp1d(Q[iBeg:iFin],Io[iBeg:iFin],bounds_error=False,fill_value=0.0)
4124    If = QI(Qf)*Qf**4
4125    R = np.linspace(0.,5000.,5000)
4126    F = fft.rfft(If)
4127    return R,F
4128
4129   
4130#### Stacking fault simulation codes ################################################################################
4131def GetStackParms(Layers):
4132   
4133    Parms = []
4134#cell parms
4135    if Layers['Laue'] in ['-3','-3m','4/m','4/mmm','6/m','6/mmm']:
4136        Parms.append('cellA')
4137        Parms.append('cellC')
4138    else:
4139        Parms.append('cellA')
4140        Parms.append('cellB')
4141        Parms.append('cellC')
4142        if Layers['Laue'] != 'mmm':
4143            Parms.append('cellG')
4144#Transition parms
4145    for iY in range(len(Layers['Layers'])):
4146        for iX in range(len(Layers['Layers'])):
4147            Parms.append('TransP;%d;%d'%(iY,iX))
4148            Parms.append('TransX;%d;%d'%(iY,iX))
4149            Parms.append('TransY;%d;%d'%(iY,iX))
4150            Parms.append('TransZ;%d;%d'%(iY,iX))
4151    return Parms
4152
4153def StackSim(Layers,ctrls,scale=0.,background={},limits=[],inst={},profile=[]):
4154    '''Simulate powder or selected area diffraction pattern from stacking faults using DIFFaX
4155   
4156    :param dict Layers: dict with following items
4157
4158      ::
4159
4160       {'Laue':'-1','Cell':[False,1.,1.,1.,90.,90.,90,1.],
4161       'Width':[[10.,10.],[False,False]],'Toler':0.01,'AtInfo':{},
4162       'Layers':[],'Stacking':[],'Transitions':[]}
4163       
4164    :param str ctrls: controls string to be written on DIFFaX controls.dif file
4165    :param float scale: scale factor
4166    :param dict background: background parameters
4167    :param list limits: min/max 2-theta to be calculated
4168    :param dict inst: instrument parameters dictionary
4169    :param list profile: powder pattern data
4170   
4171    Note that parameters all updated in place   
4172    '''
4173    import atmdata
4174    path = sys.path
4175    for name in path:
4176        if 'bin' in name:
4177            DIFFaX = name+'/DIFFaX.exe'
4178            G2fil.G2Print (' Execute '+DIFFaX)
4179            break
4180    # make form factor file that DIFFaX wants - atom types are GSASII style
4181    sf = open('data.sfc','w')
4182    sf.write('GSASII special form factor file for DIFFaX\n\n')
4183    atTypes = list(Layers['AtInfo'].keys())
4184    if 'H' not in atTypes:
4185        atTypes.insert(0,'H')
4186    for atType in atTypes:
4187        if atType == 'H': 
4188            blen = -.3741
4189        else:
4190            blen = Layers['AtInfo'][atType]['Isotopes']['Nat. Abund.']['SL'][0]
4191        Adat = atmdata.XrayFF[atType]
4192        text = '%4s'%(atType.ljust(4))
4193        for i in range(4):
4194            text += '%11.6f%11.6f'%(Adat['fa'][i],Adat['fb'][i])
4195        text += '%11.6f%11.6f'%(Adat['fc'],blen)
4196        text += '%3d\n'%(Adat['Z'])
4197        sf.write(text)
4198    sf.close()
4199    #make DIFFaX control.dif file - future use GUI to set some of these flags
4200    cf = open('control.dif','w')
4201    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4202        x0 = profile[0]
4203        iBeg = np.searchsorted(x0,limits[0])
4204        iFin = np.searchsorted(x0,limits[1])+1
4205        if iFin-iBeg > 20000:
4206            iFin = iBeg+20000
4207        Dx = (x0[iFin]-x0[iBeg])/(iFin-iBeg)
4208        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4209        cf.write('%.6f %.6f %.6f\n1\n1\nend\n'%(x0[iBeg],x0[iFin],Dx))
4210    else:
4211        cf.write('GSASII-DIFFaX.dat\n'+ctrls)
4212        inst = {'Type':['XSC','XSC',]}
4213    cf.close()
4214    #make DIFFaX data file
4215    df = open('GSASII-DIFFaX.dat','w')
4216    df.write('INSTRUMENTAL\n')
4217    if 'X' in inst['Type'][0]:
4218        df.write('X-RAY\n')
4219    elif 'N' in inst['Type'][0]:
4220        df.write('NEUTRON\n')
4221    if ctrls == '0\n0\n3\n' or ctrls == '0\n1\n3\n': 
4222        df.write('%.4f\n'%(G2mth.getMeanWave(inst)))
4223        U = ateln2*inst['U'][1]/10000.
4224        V = ateln2*inst['V'][1]/10000.
4225        W = ateln2*inst['W'][1]/10000.
4226        HWHM = U*nptand(x0[iBeg:iFin]/2.)**2+V*nptand(x0[iBeg:iFin]/2.)+W
4227        HW = np.sqrt(np.mean(HWHM))
4228    #    df.write('PSEUDO-VOIGT 0.015 -0.0036 0.009 0.605 TRIM\n')
4229        if 'Mean' in Layers['selInst']:
4230            df.write('GAUSSIAN %.6f TRIM\n'%(HW))     #fast option - might not really matter
4231        elif 'Gaussian' in Layers['selInst']:
4232            df.write('GAUSSIAN %.6f %.6f %.6f TRIM\n'%(U,V,W))    #slow - make a GUI option?
4233        else:
4234            df.write('None\n')
4235    else:
4236        df.write('0.10\nNone\n')
4237    df.write('STRUCTURAL\n')
4238    a,b,c = Layers['Cell'][1:4]
4239    gam = Layers['Cell'][6]
4240    df.write('%.4f %.4f %.4f %.3f\n'%(a,b,c,gam))
4241    laue = Layers['Laue']
4242    if laue == '2/m(ab)':
4243        laue = '2/m(1)'
4244    elif laue == '2/m(c)':
4245        laue = '2/m(2)'
4246    if 'unknown' in Layers['Laue']:
4247        df.write('%s %.3f\n'%(laue,Layers['Toler']))
4248    else:   
4249        df.write('%s\n'%(laue))
4250    df.write('%d\n'%(len(Layers['Layers'])))
4251    if Layers['Width'][0][0] < 1. or Layers['Width'][0][1] < 1.:
4252        df.write('%.1f %.1f\n'%(Layers['Width'][0][0]*10000.,Layers['Width'][0][0]*10000.))    #mum to A
4253    layerNames = []
4254    for layer in Layers['Layers']:
4255        layerNames.append(layer['Name'])
4256    for il,layer in enumerate(Layers['Layers']):
4257        if layer['SameAs']:
4258            df.write('LAYER %d = %d\n'%(il+1,layerNames.index(layer['SameAs'])+1))
4259            continue
4260        df.write('LAYER %d\n'%(il+1))
4261        if '-1' in layer['Symm']:
4262            df.write('CENTROSYMMETRIC\n')
4263        else:
4264            df.write('NONE\n')
4265        for ia,atom in enumerate(layer['Atoms']):
4266            [name,atype,x,y,z,frac,Uiso] = atom
4267            if '-1' in layer['Symm'] and [x,y,z] == [0.,0.,0.]:
4268                frac /= 2.
4269            df.write('%4s %3d %.5f %.5f %.5f %.4f %.2f\n'%(atype.ljust(6),ia,x,y,z,78.9568*Uiso,frac))
4270    df.write('STACKING\n')
4271    df.write('%s\n'%(Layers['Stacking'][0]))
4272    if 'recursive' in Layers['Stacking'][0]:
4273        df.write('%s\n'%Layers['Stacking'][1])
4274    else:
4275        if 'list' in Layers['Stacking'][1]:
4276            Slen = len(Layers['Stacking'][2])
4277            iB = 0
4278            iF = 0
4279            while True:
4280                iF += 68
4281                if iF >= Slen:
4282                    break
4283                iF = min(iF,Slen)
4284                df.write('%s\n'%(Layers['Stacking'][2][iB:iF]))
4285                iB = iF
4286        else:
4287            df.write('%s\n'%Layers['Stacking'][1])   
4288    df.write('TRANSITIONS\n')
4289    for iY in range(len(Layers['Layers'])):
4290        sumPx = 0.
4291        for iX in range(len(Layers['Layers'])):
4292            p,dx,dy,dz = Layers['Transitions'][iY][iX][:4]
4293            p = round(p,3)
4294            df.write('%.3f %.5f %.5f %.5f\n'%(p,dx,dy,dz))
4295            sumPx += p
4296        if sumPx != 1.0:    #this has to be picky since DIFFaX is.
4297            G2fil.G2Print ('ERROR - Layer probabilities sum to %.3f DIFFaX will insist it = 1.0'%sumPx)
4298            df.close()
4299            os.remove('data.sfc')
4300            os.remove('control.dif')
4301            os.remove('GSASII-DIFFaX.dat')
4302            return       
4303    df.close()   
4304    time0 = time.time()
4305    try:
4306        subp.call(DIFFaX)
4307    except OSError:
4308        G2fil.G2Print('DIFFax.exe is not available for this platform',mode='warn')
4309    G2fil.G2Print (' DIFFaX time = %.2fs'%(time.time()-time0))
4310    if os.path.exists('GSASII-DIFFaX.spc'):
4311        Xpat = np.loadtxt('GSASII-DIFFaX.spc').T
4312        iFin = iBeg+Xpat.shape[1]
4313        bakType,backDict,backVary = SetBackgroundParms(background)
4314        backDict['Lam1'] = G2mth.getWave(inst)
4315        profile[4][iBeg:iFin] = getBackground('',backDict,bakType,inst['Type'][0],profile[0][iBeg:iFin])[0]   
4316        profile[3][iBeg:iFin] = Xpat[-1]*scale+profile[4][iBeg:iFin]
4317        if not np.any(profile[1]):                   #fill dummy data x,y,w,yc,yb,yd
4318            rv = st.poisson(profile[3][iBeg:iFin])
4319            profile[1][iBeg:iFin] = rv.rvs()
4320            Z = np.ones_like(profile[3][iBeg:iFin])
4321            Z[1::2] *= -1
4322            profile[1][iBeg:iFin] = profile[3][iBeg:iFin]+np.abs(profile[1][iBeg:iFin]-profile[3][iBeg:iFin])*Z
4323            profile[2][iBeg:iFin] = np.where