Changeset 4415 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
May 7, 2020 12:17:26 PM (5 years ago)
Author:
vondreele
Message:

change export_PDB to handle fullrmc pdb files
implement Wenqian's mod for the geometric image correction
TransformPhase? & FillUnitCell? now have option to not force atoms to new unit cell (default=True)
refactor FillAtomLookup?
change FindAllNeighbors? to optionally give short output
remove result of double click on Uiso column heading
many additions & changes for fulrmc GUI, etc.
comment out obsolete MPLsubplots from G2plot - no longer needed (pre mpl 2.2)
modify PDB importer to handle fullrmc pdb files.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified trunk/GSASIIpwd.py

    r4405 r4415  
    1818import time
    1919import os
     20import os.path
    2021import subprocess as subp
    2122import copy
    … …  
    24622463    return fname
    24632464
     2465def FindBonds(Phase,RMCPdict):
     2466    generalData = Phase['General']
     2467    cx,ct,cs,cia = generalData['AtomPtrs']
     2468    atomData = Phase['Atoms']
     2469    Res = 'RMC'
     2470    if 'macro' in generalData['Type']:
     2471        Res = atomData[0][ct-3]
     2472    AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
     2473    Pairs = RMCPdict['Pairs']   #dict!
     2474    BondList = []
     2475    notNames = []
     2476    for FrstName in AtDict:
     2477        nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
     2478        Atyp1 = AtDict[FrstName]
     2479        for nbr in nbrs:
     2480            Atyp2 = AtDict[nbr[0]]
     2481            try:
     2482                bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
     2483            except KeyError:
     2484                bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
     2485            if any(bndData):
     2486                if bndData[0] <= nbr[1] <= bndData[1]:
     2487                    bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
     2488                    revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
     2489                    if bondStr not in BondList and revbondStr not in BondList:
     2490                        BondList.append(bondStr)
     2491        notNames.append(FrstName)
     2492    return Res,BondList
     2493
     2494def FindAngles(Phase,RMCPdict):
     2495    generalData = Phase['General']
     2496    Cell = generalData['Cell'][1:7]
     2497    Amat = G2lat.cell2AB(Cell)[0]
     2498    cx,ct,cs,cia = generalData['AtomPtrs']
     2499    atomData = Phase['Atoms']
     2500    AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
     2501    AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
     2502    Angles = RMCPdict['Angles']
     2503    AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
     2504    AngleList = []
     2505    for MidName in AtDict:
     2506        nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
     2507        if len(nbrs) < 2: #need 2 neighbors to make an angle
     2508            continue
     2509        Atyp2 = AtDict[MidName]
     2510        for i,nbr1 in enumerate(nbrs):
     2511            Atyp1 = AtDict[nbr1[0]]
     2512            for j,nbr3 in enumerate(nbrs[i+1:]):
     2513                Atyp3 = AtDict[nbr3[0]]
     2514                IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
     2515                try:
     2516                    angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
     2517                except KeyError:
     2518                    try:
     2519                        angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
     2520                    except KeyError:
     2521                        continue
     2522                XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
     2523                calAngle = G2mth.getRestAngle(XYZ,Amat)
     2524                if angData[0] < calAngle < angData[1]:
     2525                    AngleList.append(str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n')
     2526    return AngleList
     2527
     2528def GetSqConvolution(XY,d):
     2529
     2530    n = XY.shape[1]
     2531    snew = np.zeros(n)
     2532    dq = np.zeros(n)
     2533    sold = XY[1]
     2534    q = XY[0]
     2535    dq[1:] = np.diff(q)
     2536    dq[0] = dq[1]
     2537   
     2538    for j in range(n):
     2539        for i in range(n):
     2540            b = abs(q[i]-q[j])
     2541            t = q[i]+q[j]
     2542            if j == i:
     2543                snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
     2544            else:
     2545                snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
     2546        snew[j] /= np.pi*q[j]
     2547   
     2548    snew[0] = snew[1]
     2549    return snew
     2550
     2551def GetMaxSphere(pdbName):
     2552    try:
     2553        pFil = open(pdbName,'r')
     2554    except FileNotFoundError:
     2555        return None
     2556    while True:
     2557        line = pFil.readline()
     2558        if 'Boundary' in line:
     2559            line = line.split()[3:]
     2560            G = np.array([float(item) for item in line])
     2561            G = np.reshape(G,(3,3))**2
     2562            G = nl.inv(G)
     2563            pFil.close()
     2564            break
     2565    dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
     2566    return min(dspaces)
     2567   
    24642568def MakefullrmcRun(pName,Phase,RMCPdict):
     2569    Res,BondList = FindBonds(Phase,RMCPdict)
     2570    AngleList = FindAngles(Phase,RMCPdict)
     2571    rmin = RMCPdict['min Contact']
     2572    rmax = GetMaxSphere(RMCPdict['atomPDB'])
     2573    if rmax is None:
     2574        return None
    24652575    rname = pName+'-run.py'
     2576    restart = '%s_restart.pdb'%pName
    24662577    Files = RMCPdict['files']
     2578    wtDict = {}
    24672579    rundata = ''
    24682580    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
    … …  
    24702582# fullrmc imports (all that are potentially useful)
    24712583import numpy as np
     2584from fullrmc.sincConvolution import sincConvolution
    24722585from fullrmc.Globals import LOGGER
    24732586from fullrmc.Engine import Engine
    24742587from fullrmc.Constraints.PairDistributionConstraints import PairDistributionConstraint
    24752588from fullrmc.Constraints.PairCorrelationConstraints import PairCorrelationConstraint
    2476 from fullrmc.Constraints.StructureFactorConstraints import StructureFactorConstraint
     2589from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint
    24772590from fullrmc.Constraints.DistanceConstraints import InterMolecularDistanceConstraint
    24782591from fullrmc.Constraints.BondConstraints import BondConstraint
    24792592from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
    2480 from fullrmc.Constraints.ImproperAngleConstraints import ImproperAngleConstraint
     2593from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
    24812594from fullrmc.Core.MoveGenerator import MoveGeneratorCollector
    24822595from fullrmc.Core.GroupSelector import RecursiveGroupSelector
    … …  
    24882601from fullrmc.Core.Collection import get_principal_axis
    24892602# engine setup\n'''
    2490     rundata += 'LOGGER.set_log_file_basename(%s)\n'%pName
     2603    rundata += 'LOGGER.set_log_file_basename("%s")\n'%pName
    24912604    rundata += 'engineFileName = "%s.rmc"\n'%pName
    24922605    rundata += 'ENGINE = Engine(path=None)\n'
    24932606    rundata += 'if not ENGINE.is_engine(engineFileName):\n'
    2494 # create engine
     2607    rundata += '# create engine & set atomic (pdb) model\n'
    24952608    rundata += '    ENGINE = Engine(path=engineFileName)\n'
    2496     rundata += '    ENGINE.set_pdb(%s)\n'%RMCPdict['atomPDB']
    2497 ## create experimental constraints
     2609    rundata += '    ENGINE.set_pdb("%s")\n'%RMCPdict['atomPDB']
     2610    rundata += '# create experimental constraints must restart to change these\n'
    24982611    for File in Files:
    24992612        filDat = RMCPdict['files'][File]
    … …  
    25032616                sfwt = 'xrays'
    25042617            if 'G(r)' in File:
    2505                 rundata += '    GofR = PairDistributionConstraint(experimentalData=%s, weighting="%s")\n'%(filDat[0],sfwt)
    2506                 rundata += '    GofR.set_variance_squared(%f)\n'%filDat[1]
     2618                rundata += '    RGR = np.loadtxt("%s").T\n'%filDat[0]
     2619                if filDat[3]:
     2620                    rundata += '    RGR[1] *= RGR[0]\n'
     2621                rundata += '    GofR = PairDistributionConstraint(experimentalData=RGR.T, weighting="%s")\n'%sfwt
    25072622                rundata += '    ENGINE.add_constraints([GofR])\n'
     2623                wtDict['Pair'] = filDat[1]
    25082624            else:
    2509                 rundata += '    FofQ = StructureFactorConstraint(experimentalData=%s, weighting="%s")\n'%(filDat[0],sfwt)
    2510                 rundata += '    FofQ.set_variance_squared(%f)\n'%filDat[1]
     2625                rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
     2626                if filDat[3]:
     2627                    rundata += '    SOQ[1] = sincConvolution(SOQ,%.3f)\n'%rmax
     2628                rundata += '    FofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
    25112629                rundata += '    ENGINE.add_constraints([FofQ])\n'
     2630                wtDict['Struct'] = filDat[1]
     2631    rundata += '    ENGINE.add_constraints(InterMolecularDistanceConstraint())\n'
     2632    if len(BondList):
     2633        rundata += '    B_CONSTRAINT   = BondConstraint()\n'
     2634        rundata += '    ENGINE.add_constraints(B_CONSTRAINT)\n'
     2635    if len(AngleList):
     2636        rundata += '    A_CONSTRAINT   = BondsAngleConstraint()\n'
     2637        rundata += '    ENGINE.add_constraints(A_CONSTRAINT)\n'
     2638    if len(RMCPdict['Torsions']):
     2639        rundata += '    T_CONSTRAINT   = DihedralAngleConstraint()\n'
     2640        rundata += '    ENGINE.add_constraints(T_CONSTRAINT)\n'
    25122641    rundata += '    ENGINE.save()\n'
    2513        
    2514        
    2515    
    25162642    rundata += 'else:\n'
    25172643    rundata += '    ENGINE = ENGINE.load(path=engineFileName)\n'
    2518 
     2644    rundata += '#fill & change constraints - can be done without restart\n' 
     2645    rundata += 'Constraints = ENGINE.constraints\n'
     2646    rundata += 'for constraint in Constraints:\n'
     2647    rundata += '    strcons = str(type(constraint))\n'
     2648    if len(BondList):
     2649        rundata += '    if "BondConstraint" in strcons:\n'
     2650        rundata += '        constraint.set_variance_squared(%f)\n'%RMCPdict['Bond Weight']
     2651        rundata += '        constraint.create_bonds_by_definition(bondsDefinition={"%s":[\n'%Res
     2652        for bond in BondList:
     2653            rundata += '        %s'%bond
     2654        rundata += '        ]})\n'
     2655    if len(AngleList):
     2656        rundata += '    elif "BondsAngleConstraint" in strcons:\n'
     2657        rundata += '        constraint.set_variance_squared(%f)\n'%RMCPdict['Angle Weight']
     2658        rundata += '        constraint.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
     2659        for angle in AngleList:
     2660            rundata += '        %s'%angle
     2661        rundata += '        ]})\n'
     2662    if len(RMCPdict['Torsions']):
     2663        rundata += '    elif "DihedralAngleConstraint" in strcons:\n'
     2664        rundata += '        constraint.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
     2665        rundata += '        constraint.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
     2666        for torsion in RMCPdict['Torsions']:
     2667            rundata += '    %s\n'%str(tuple(torsion))
     2668        rundata += '        ]})\n'
     2669    rundata += '    elif "InterMolecular" in strcons:\n'
     2670    rundata += '        constraint.set_default_distance(%f)\n'%RMCPdict['min Contact']
     2671    rundata += '    elif "PairDistribution" in strcons:\n'
     2672    rundata += '        constraint.set_variance_squared(%f)\n'%wtDict['Pair']
     2673    rundata += '        constraint.set_limits((%.3f,%.3f))\n'%(rmin,rmax)
     2674    if RMCPdict['FitScale']:
     2675        rundata += '        constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
     2676    rundata += '    elif "StructureFactor" in strcons:\n'
     2677    rundata += '        constraint.set_variance_squared(%f)\n'%wtDict['Struct']
     2678    if RMCPdict['FitScale']:
     2679        rundata += '        constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
     2680    rundata += 'ENGINE.set_groups_as_atoms()\n'
     2681    if len(RMCPdict['Swaps']):
     2682        print(RMCPdict['Swaps'])
     2683
     2684    # allElements = ENGINE.allElements
     2685    # niSwapList = [[idx] for idx in range(len(allElements)) if allElements[idx]=='ni']
     2686    # tiSwapList = [[idx] for idx in range(len(allElements)) if allElements[idx]=='ti']
     2687    # # create swap generator
     2688    # toNiSG = SwapPositionsGenerator(swapList=niSwapList)
     2689    # toTiSG = SwapPositionsGenerator(swapList=tiSwapList)
     2690       
     2691
     2692    rundata += 'ENGINE.run(restartPdb="%s",numberOfSteps=%d, saveFrequency=1000)\n'%(restart,RMCPdict['Cycles'])
     2693    rundata += 'ENGINE.save()\n'
     2694    rundata += '#setup runs for fullrmc\n'
     2695   
     2696   
     2697   
     2698   
     2699   
     2700       
     2701   
    25192702
    25202703    rfile = open(rname,'w')
    … …  
    25232706   
    25242707    return rname
    2525    
     2708
    25262709
    25272710def MakefullrmcPDB(Name,Phase,RMCPdict):
    … …  
    25342717    newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
    25352718    newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
    2536     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False)
     2719    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=False)
    25372720    Atoms = newPhase['Atoms']
    25382721    XYZ = np.array([atom[3:6] for atom in Atoms]).T
    … …  
    25422725    fname = Name+'_cbb.pdb'
    25432726    fl = open(fname,'w')
    2544     fl.write('REMARK    this file is generated using GSASII\n')
    25452727    fl.write('REMARK    Boundary Conditions:%6.2f  0.0  0.0  0.0%7.2f  0.0  0.0  0.0%7.2f\n'%(
    25462728             Cell[0],Cell[1],Cell[2]))
    2547     fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
    2548             Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
    25492729    fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
    25502730    fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
    25512731    fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
     2732    fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
     2733            Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
    25522734
    25532735    Natm = np.core.defchararray.count(np.array(Atcodes),'+')
    … …  
    25572739        NPM = RMCPdict['Natoms']
    25582740        for iat,atom in enumerate(Atoms):
    2559             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
    2560             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %-2s\n'%(       
     2741            XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
     2742            fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    25612743                    1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    25622744            nat += 1
    … …  
    25662748                if atom[1] == atm:
    25672749                    XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
    2568                     fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %-2s\n'%(       
     2750                    fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    25692751                            1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    25702752                    nat += 1
    … …  
    37123894            reflDict[hash('%5d%5d%5d'%(ref[0],ref[1],ref[2]))] = iref
    37133895    fbaName = os.path.splitext(prfName)[0]+'.fba'
    3714     try: # patch for FileNotFoundError not in Python 2.7
    3715         FileNotFoundError
    3716     except NameError:
    3717         FileNotFoundError = Exception
    3718     try:
     3896    if os.path.isfile(fbaName):
    37193897        fba = open(fbaName,'r')
    3720     except FileNotFoundError:
     3898    else:
    37213899        return False
    37223900    fba.readline()
Note: See TracChangeset for help on using the changeset viewer.