source: trunk/GSASIIpwd.py @ 4466

Last change on this file since 4466 was 4466, checked in by vondreele, 17 months ago

change Pawley refinement menu items for setting/unsetting refine flags - now more obvious
tune up monoclinic indexing - finer volume step & more tries at each step - improves jadarite indexing

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