source: trunk/GSASIIpwd.py @ 5281

Last change on this file since 5281 was 5281, checked in by toby, 7 months ago

start putting paths into PDFfit script

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