Changeset 3990 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
May 22, 2019 2:04:11 PM (2 years ago)
Author:
vondreele
Message:

MEM Dysnomia now works but can be uncontrolled. No feedback into GSAS-II yet

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r3959 r3990  
    28102810    Layers['Sadp']['Img'] = Sapd
    28112811    print (' GETSAD time = %.2fs'%(time.time()-time0))
    2812 #    GSASIIpath.IPyBreak()
    2813    
     2812   
     2813###############################################################################
     2814#### Maximum Entropy Method - Dysnomia
     2815###############################################################################
     2816   
     2817def makePRFfile(data,MEMtype):
     2818    ''' makes Dysnomia .prf control file from Dysnomia GUI controls
     2819   
     2820    ;param dict data: GSAS-II phase data
     2821    :param int MEMtype: 1 for neutron data with negative scattering lengths
     2822                        0 otherwise
     2823    :returns str: name of Dysnomia control file
     2824    '''
     2825
     2826    generalData = data['General']
     2827    pName = generalData['Name'].replace(' ','_')
     2828    DysData = data['Dysnomia']
     2829    prfName = pName+'.prf'
     2830    prf = open(prfName,'w')
     2831    prf.write('$PREFERENCES\n')
     2832    prf.write(pName+'.mem\n') #or .fos?
     2833    prf.write(pName+'.out\n')
     2834    prf.write(pName+'.pgrid\n')
     2835    prf.write(pName+'.fba\n')
     2836    prf.write(pName+'_eps.raw\n')
     2837    prf.write('%d\n'%MEMtype)
     2838    if DysData['DenStart'] == 'uniform':
     2839        prf.write('0\n')
     2840    else:
     2841        prf.write('1\n')
     2842    if DysData['Optimize'] == 'ZSPA':
     2843        prf.write('0\n')
     2844    else:
     2845        prf.write('1\n')
     2846    prf.write('1\n')
     2847    if DysData['Lagrange'][0] == 'user':
     2848        prf.write('0\n')
     2849    else:
     2850        prf.write('1\n')
     2851    prf.write('%.4f %d\n'%(DysData['Lagrange'][1],DysData['wt pwr']))
     2852    prf.write('%.3f\n'%DysData['Lagrange'][2])
     2853    prf.write('%.2f\n'%DysData['E_factor'])
     2854    prf.write('1\n')
     2855    prf.write('0\n')
     2856    prf.write('%d\n'%DysData['Ncyc'])
     2857    prf.write('1\n')
     2858    prf.write('1 0 0 0 0 0 0 0\n')
     2859    if DysData['prior'] == 'uniform':
     2860        prf.write('0\n')
     2861    else:
     2862        prf.write('1\n')
     2863    prf.close()
     2864    return prfName
     2865
     2866def makeMEMfile(data,reflData,MEMtype):
     2867    ''' make Dysnomia .mem file of reflection data, etc.
     2868    ;param dict data: GSAS-II phase data
     2869    :param list reflData: GSAS-II reflection data
     2870    :param int MEMtype: 1 for neutron data with negative scattering lengths
     2871                        0 otherwise
     2872    '''
     2873   
     2874    DysData = data['Dysnomia']
     2875    generalData = data['General']
     2876    pName = generalData['Name'].replace(' ','_')
     2877    memName = pName+'.mem'
     2878    Map = generalData['Map']
     2879    Type = Map['Type']
     2880    UseList = Map['RefList']
     2881    mem = open(memName,'w')
     2882    mem.write('%s\n'%(generalData['Name']+' from '+UseList[0]))
     2883    a,b,c,alp,bet,gam = generalData['Cell'][1:7]
     2884    mem.write('%10.5f%10.5f%10.5f%10.5f,%10.5f%10.5f\n'%(a,b,c,alp,bet,gam))
     2885    mem.write('      0.0000000      0.0000000     -1    0    0    0     P\n')   #dummy PO stuff
     2886    SGSym = generalData['SGData']['SpGrp']
     2887    try:
     2888        SGId = G2spc.spgbyNum.index(SGSym)
     2889    except IndexError:
     2890        return False
     2891    org = 1
     2892    if SGSym in G2spc.spg2origins:
     2893        org = 2
     2894    mapsize = Map['rho'].shape
     2895    sumZ = 0.
     2896    sumpos = 0.
     2897    sumneg = 0.
     2898    mem.write('%5d%5d%5d%5d%5d\n'%(SGId,org,mapsize[0],mapsize[1],mapsize[2]))
     2899    for atm in generalData['NoAtoms']:
     2900        Nat = generalData['NoAtoms'][atm]
     2901        AtInfo = G2elem.GetAtomInfo(atm)
     2902        sumZ += Nat*AtInfo['Z']
     2903        isotope = generalData['Isotope'][atm]
     2904        blen = generalData['Isotopes'][atm][isotope]['SL'][0]
     2905        if blen < 0.:
     2906            sumneg += blen*Nat
     2907        else:
     2908            sumpos += blen*Nat
     2909    if 'X' in Type:
     2910        mem.write('%10.2f  0.001\n'%sumZ)
     2911    elif 'N' in Type and MEMtype:
     2912        mem.write('%10.3f%10.3f 0.001\n'%(sumpos,sumneg))
     2913    else:
     2914        mem.write('%10.3 0.001\n'%sumpos)
     2915       
     2916    refs = []
     2917    prevpos = 0.
     2918    for ref in reflData:
     2919        if 'T' in Type:
     2920            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,x,x,x,prfo = ref[:16]
     2921            FWHM = getgamFW(gam,sig)
     2922        else:
     2923            h,k,l,mult,dsp,pos,sig,gam,Fobs,Fcalc,phase,x,prfo = ref[:13]
     2924            g = gam/100.    #centideg -> deg
     2925            s = np.sqrt(max(sig,0.0001))/100.   #var -> sig in deg
     2926            FWHM = getgamFW(g,s)
     2927        delt = pos-prevpos
     2928        refs.append([h,k,l,mult,pos,FWHM,Fobs,phase,delt])
     2929        prevpos = pos
     2930           
     2931    ovlp = DysData['overlap']
     2932    refs1 = []
     2933    refs2 = []
     2934    nref2 = 0
     2935    iref = 0
     2936    Nref = len(refs)
     2937    start = False
     2938    while iref < Nref-1:
     2939        if refs[iref+1][-1] < ovlp*refs[iref][5]:
     2940            if refs[iref][-1] > ovlp*refs[iref][5]:
     2941                refs2.append([])
     2942                start = True
     2943            if nref2 == len(refs2):
     2944                refs2.append([])
     2945            refs2[nref2].append(refs[iref])
     2946        else:
     2947            if start:
     2948                refs2[nref2].append(refs[iref])
     2949                start = False
     2950                nref2 += 1
     2951            else:
     2952                refs1.append(refs[iref])
     2953        iref += 1
     2954    if start:
     2955        refs2[nref2].append(refs[iref])
     2956    else:
     2957        refs1.append(refs[iref])
     2958   
     2959    mem.write('%5d\n'%len(refs1))
     2960    for ref in refs1:
     2961        h,k,l = ref[:3]
     2962        Fobs = np.sqrt(ref[6])
     2963        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)))
     2964    mem.write('%5d\n'%len(refs2))
     2965    for ref2 in refs2:
     2966        if not len(ref2):
     2967            break
     2968        mem.write('%5d\n'%len(ref2))
     2969        Gsum = 0.
     2970        Msum = 0
     2971        for ref in ref2:
     2972            Gsum += ref[6]*ref[3]
     2973            Msum += ref[3]
     2974        G = np.sqrt(Gsum/Msum)
     2975        h,k,l = ref2[0][:3]
     2976        mem.write('%5d%5d%5d%10.3f%10.3f%5d\n'%(h,k,l,G,max(0.01*G,0.1),ref2[0][3]))
     2977        for ref in ref2[1:]:
     2978            h,k,l,m = ref[:4]
     2979            mem.write('%5d%5d%5d%5d\n'%(h,k,l,m))
     2980    mem.write('0\n')
     2981    mem.close()
     2982    return True
    28142983#testing data
    28152984NeedTestData = True
Note: See TracChangeset for help on using the changeset viewer.