source: trunk/GSASIIpwd.py

Last change on this file was 5478, checked in by vondreele, 2 weeks ago

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