Changeset 4951 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Jun 11, 2021 8:42:20 AM (2 years ago)
Author:
toby
Message:

working fullrmc implementation. Much more to do (see TODO in phsGUI), but good start. Needs external fullrmc Python from Bachir until 5.0 is released.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r4919 r4951  
    576576                data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg))
    577577        else:
    578             G2fil.G2Print('  end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f}) *** {} ***\n'.format(
     578            G2fil.G2Print('  end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} RMS:{:.4f}) *** {} ***\n'.format(
    579579                data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg))
    580580    return res
    … …  
    28802880#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
    28812881#     return min(dspaces)
    2882    
     2882
     2883def findfullrmc():
     2884    '''Find where fullrmc is installed. Tries the following:
     2885   
     2886         1. Returns the Config var 'fullrmc_exec', if defined. No check
     2887            is done that the interpreter has fullrmc
     2888         2. The current Python interpreter if fullrmc can be imported
     2889            and fullrmc is version 5+
     2890         3. The path is checked for a fullrmc image as named by Bachir
     2891
     2892    :returns: the full path to a python executable that is assumed to
     2893      have fullrmc installed or None, if it was not found.
     2894    '''
     2895    is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK)
     2896    if GSASIIpath.GetConfigValue('fullrmc_exec') is not None and is_exe(
     2897            GSASIIpath.GetConfigValue('fullrmc_exec')):
     2898        return GSASIIpath.GetConfigValue('fullrmc_exec')
     2899    try:
     2900        import fullrmc
     2901        if int(fullrmc.__version__.split('.')[0]) >= 5:
     2902            return sys.executable
     2903    except:
     2904        pass
     2905    pathlist = os.environ["PATH"].split(os.pathsep)
     2906    for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(),
     2907                  os.path.split(sys.executable)[0]):
     2908        if p not in pathlist: pathlist.insert(0,p)
     2909    import glob
     2910    for p in pathlist:
     2911        if sys.platform == "darwin":
     2912            lookfor = "fullrmc*macOS*i386-64bit"
     2913        elif sys.platform == "win32":
     2914            lookfor = "fullrmc*.exe"
     2915        else:
     2916            lookfor = "fullrmc*"
     2917        fl = glob.glob(lookfor)
     2918        if len(fl) > 0:
     2919            return os.path.abspath(sorted(fl)[0])
     2920       
    28832921def MakefullrmcRun(pName,Phase,RMCPdict):
     2922    '''Creates a script to run fullrmc. Returns the name of the file that was
     2923    created.
     2924    '''
    28842925    BondList = {}
    28852926    for k in RMCPdict['Pairs']:
    … …  
    28932934            angle[i] = angle[i].strip()
    28942935        AngleList.append(angle)
    2895     rmin = RMCPdict['min Contact']
     2936    # rmin = RMCPdict['min Contact']
    28962937    cell = Phase['General']['Cell'][1:7]
    28972938    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
    … …  
    29012942        el = ''.join([i for i in atom[ct] if i.isalpha()])
    29022943        atomsList.append([el] + atom[cx:cx+4])
    2903     rname = pName+'-fullrmc.py'
     2944    projDir,projName = os.path.split(pName)
     2945    scrname = pName+'-fullrmc.py'
    29042946    restart = '%s_restart.pdb'%pName
    29052947    Files = RMCPdict['files']
    29062948    rundata = ''
    2907     rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
     2949    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%scrname
    29082950    rundata += '# created in '+__file__+" v"+filversion.split()[1]
    29092951    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
    … …  
    29122954import os,glob
    29132955import time
    2914 #import matplotlib.pyplot as plt
     2956import pickle
    29152957import numpy as np
    29162958from fullrmc.Core import Collection
    … …  
    29232965from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
    29242966from fullrmc.Generators.Swaps import SwapPositionsGenerator
    2925 
    2926 ### When True, erases an existing enging to provide a fresh start
    2927 FRESH_START = {}
     2967def writeHeader(ENGINE,statFP):
     2968    statFP.write('generated-steps, total-error, ')
     2969    for c in ENGINE.constraints:
     2970        statFP.write(c.constraintName)
     2971        statFP.write(', ')
     2972    statFP.write('\\n')
     2973    statFP.flush()
     2974   
     2975def writeCurrentStatus(ENGINE,statFP,plotF):
     2976    statFP.write(str(ENGINE.generated))
     2977    statFP.write(', ')
     2978    statFP.write(str(ENGINE.totalStandardError))
     2979    statFP.write(', ')
     2980    for c in ENGINE.constraints:
     2981        statFP.write(str(c.standardError))
     2982        statFP.write(', ')
     2983    statFP.write('\\n')
     2984    statFP.flush()
     2985    mpl.use('agg')
     2986    fp = open(plotF,'wb')
     2987    for c in ENGINE.constraints:
     2988        p = c.plot(show=False)
     2989        p[0].canvas.draw()
     2990        image = p[0].canvas.buffer_rgba()
     2991        pickle.dump(c.constraintName,fp)
     2992        pickle.dump(np.array(image),fp)
     2993    fp.close()
     2994'''
     2995    rundata += '''
     2996### When True, erases an existing engine to provide a fresh start
     2997FRESH_START = {:}
     2998dirName = "{:}"
     2999prefix = "{:}"
     3000project = prefix + "-fullrmc"
    29283001time0 = time.time()
    2929 '''.format(RMCPdict['ReStart'][0])
     3002'''.format(RMCPdict['ReStart'][0],projDir,projName)
    29303003   
    29313004    rundata += '# setup structure\n'
    … …  
    29363009
    29373010    rundata += '\n# initialize engine\n'
    2938     rundata += 'engineFileName = "%s.rmc"\n'%pName
    2939 
    2940     rundata += '''\n# check Engine exists if so (and not FRESH_START) load it
     3011    rundata += '''
     3012engineFileName = os.path.join(dirName, project + '.rmc')
     3013projectStats = os.path.join(dirName, project + '.stats')
     3014projectPlots = os.path.join(dirName, project + '.plots')
     3015pdbFile = os.path.join(dirName, project + '_restart.pdb')
     3016# check Engine exists if so (and not FRESH_START) load it
    29413017# otherwise build it
    29423018ENGINE = Engine(path=None)
    … …  
    29653041            sfwt = 'atomicNumber'
    29663042        if 'G(r)' in File:
    2967             rundata += '    GR = np.loadtxt("%s").T\n'%filDat[0]
     3043            rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
    29683044            if filDat[3] == 0:
    29693045                rundata += '''    # read and xform G(r) as defined in RMCProfile
    … …  
    29823058            rundata += '    GofR.set_limits((None, rmax))\n'
    29833059        elif '(Q)' in File:
    2984             rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
     3060            rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
    29853061            if filDat[3] == 0:
    29863062                rundata += '    # Read & xform F(Q) as defined in RMCProfile to S(Q)-1\n'
    … …  
    30153091               '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item))
    30163092        rundata += '             ])\n'
    3017     rundata += '    for f in glob.glob("{}_*.log"): os.remove(f)\n'.format(pName)
    30183093    rundata += '''
     3094    for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f)
    30193095    ENGINE.save()
    30203096else:
    30213097    ENGINE = ENGINE.load(path=engineFileName)
     3098    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
     3099
     3100ENGINE.set_log_file(os.path.join(dirName,prefix))
    30223101'''
    3023     rundata += 'ENGINE.set_log_file("{}")\n'.format(pName)
    30243102    if RMCPdict['Swaps']:
    30253103        rundata += '\n#set up for site swaps\n'
    … …  
    30553133    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
    30563134    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
    3057     rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
     3135    # rundata += '        c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n'
    30583136    rundata += '        c.set_limits((None,rmax))\n'
    30593137    if RMCPdict['FitScale']:
    30603138        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    3061     rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
    3062     rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
     3139    # rundata += '        c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n'
    30633140    if RMCPdict['FitScale']:
     3141        rundata += '    elif type(c) is ReducedStructureFactorConstraint:\n'
    30643142        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    30653143    # torsions difficult to implement, must be internal to cell & named with
    … …  
    30733151    #         rundata += '    %s\n'%str(tuple(torsion))
    30743152    #     rundata += '        ]})\n'           
    3075     rundata += '\n# setup runs for fullrmc\n'
    3076 
    3077     rundata += 'steps = 10000\n'
     3153    rundata += '''
     3154if FRESH_START:
     3155    # initialize engine with one step to get starting config energetics
     3156    ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1)
     3157    statFP = open(projectStats,'w')
     3158    writeHeader(ENGINE,statFP)
     3159    writeCurrentStatus(ENGINE,statFP,projectPlots)
     3160else:
     3161    statFP = open(projectStats,'a')
     3162
     3163# setup runs for fullrmc
     3164'''
     3165    rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle'])
    30783166    rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles'])
    30793167    rundata += '    ENGINE.set_groups_as_atoms()\n'
    30803168    rundata += '    expected = ENGINE.generated+steps\n'
    30813169   
    3082     rundata += '    ENGINE.run(restartPdb="%s",numberOfSteps=steps, saveFrequency=1000)\n'%restart
     3170    rundata += '    ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n'
     3171    rundata += '    writeCurrentStatus(ENGINE,statFP,projectPlots)\n'
    30833172    rundata += '    if ENGINE.generated != expected: break # run was stopped\n'
     3173    rundata += 'statFP.close()\n'
    30843174    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
    3085     rfile = open(rname,'w')
     3175    rfile = open(scrname,'w')
    30863176    rfile.writelines(rundata)
    30873177    rfile.close()
    3088     return rname
     3178    return scrname
    30893179   
    30903180def GetRMCBonds(general,RMCPdict,Atoms,bondList):
    … …  
    43244414    G2fil.G2Print ("OK")
    43254415    plotter.StartEventLoop()
     4416   
     4417#    GSASIIpath.SetBinaryPath(True,False)
     4418#    print('found',findfullrmc())
Note: See TracChangeset for help on using the changeset viewer.