source: trunk/GSASIIpwd.py @ 4588

Last change on this file since 4588 was 4588, checked in by toby, 13 months ago

use G2VarObj for param limits; add more info to seq. ref. done dialog; show Frozen in show LS parameters

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