Changeset 4518 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Jul 12, 2020 8:59:57 PM (16 months ago)
Author:
toby
Message:

start on fullrmc update for 5.0; fix for rigid bodies

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r4483 r4518  
    2020import os.path
    2121import subprocess as subp
     22import datetime as dt
    2223import copy
    2324
     
    3334
    3435import GSASIIpath
     36filversion = "$Revision$"
    3537GSASIIpath.SetVersionNumber("$Revision$")
    3638import GSASIIlattice as G2lat
     
    25112513    return fname
    25122514
    2513 def 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:
     2515# def FindBonds(Phase,RMCPdict):
     2516#     generalData = Phase['General']
     2517#     cx,ct,cs,cia = generalData['AtomPtrs']
     2518#     atomData = Phase['Atoms']
     2519#     Res = 'RMC'
     2520#     if 'macro' in generalData['Type']:
     2521#         Res = atomData[0][ct-3]
     2522#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
     2523#     Pairs = RMCPdict['Pairs']   #dict!
     2524#     BondList = []
     2525#     notNames = []
     2526#     for FrstName in AtDict:
     2527#         nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0]
     2528#         Atyp1 = AtDict[FrstName]
     2529#         if 'Va' in Atyp1:
     2530#             continue
     2531#         for nbr in nbrs:
     2532#             Atyp2 = AtDict[nbr[0]]
     2533#             if 'Va' in Atyp2:
     2534#                 continue
     2535#             try:
     2536#                 bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:]
     2537#             except KeyError:
     2538#                 bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:]
     2539#             if any(bndData):
     2540#                 if bndData[0] <= nbr[1] <= bndData[1]:
     2541#                     bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n'
     2542#                     revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n'
     2543#                     if bondStr not in BondList and revbondStr not in BondList:
     2544#                         BondList.append(bondStr)
     2545#         notNames.append(FrstName)
     2546#     return Res,BondList
     2547
     2548# def FindAngles(Phase,RMCPdict):
     2549#     generalData = Phase['General']
     2550#     Cell = generalData['Cell'][1:7]
     2551#     Amat = G2lat.cell2AB(Cell)[0]
     2552#     cx,ct,cs,cia = generalData['AtomPtrs']
     2553#     atomData = Phase['Atoms']
     2554#     AtLookup = G2mth.FillAtomLookUp(atomData,cia+8)
     2555#     AtDict = {atom[ct-1]:atom[ct] for atom in atomData}
     2556#     Angles = RMCPdict['Angles']
     2557#     AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles}
     2558#     AngleList = []
     2559#     for MidName in AtDict:
     2560#         nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True)
     2561#         if len(nbrs) < 2: #need 2 neighbors to make an angle
     2562#             continue
     2563#         Atyp2 = AtDict[MidName]
     2564#         for i,nbr1 in enumerate(nbrs):
     2565#             Atyp1 = AtDict[nbr1[0]]
     2566#             for j,nbr3 in enumerate(nbrs[i+1:]):
     2567#                 Atyp3 = AtDict[nbr3[0]]
     2568#                 IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]]
     2569#                 try:
     2570#                     angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)]
     2571#                 except KeyError:
     2572#                     try:
     2573#                         angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)]
     2574#                     except KeyError:
     2575#                         continue
     2576#                 XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3))
     2577#                 calAngle = G2mth.getRestAngle(XYZ,Amat)
     2578#                 if angData[0] <= calAngle <= angData[1]:
     2579#                     angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n'
     2580#                     revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n'
     2581#                     if angStr not in AngleList and revangStr not in AngleList:
     2582#                         AngleList.append(angStr)
     2583#     return AngleList
     2584
     2585# def GetSqConvolution(XY,d):
     2586
     2587#     n = XY.shape[1]
     2588#     snew = np.zeros(n)
     2589#     dq = np.zeros(n)
     2590#     sold = XY[1]
     2591#     q = XY[0]
     2592#     dq[1:] = np.diff(q)
     2593#     dq[0] = dq[1]
     2594   
     2595#     for j in range(n):
     2596#         for i in range(n):
     2597#             b = abs(q[i]-q[j])
     2598#             t = q[i]+q[j]
     2599#             if j == i:
     2600#                 snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i]
     2601#             else:
     2602#                 snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i]
     2603#         snew[j] /= np.pi*q[j]
     2604   
     2605#     snew[0] = snew[1]
     2606#     return snew
     2607
     2608# def GetMaxSphere(pdbName):
     2609#     try:
     2610#         pFil = open(pdbName,'r')
     2611#     except FileNotFoundError:
     2612#         return None
     2613#     while True:
     2614#         line = pFil.readline()
     2615#         if 'Boundary' in line:
     2616#             line = line.split()[3:]
     2617#             G = np.array([float(item) for item in line])
     2618#             G = np.reshape(G,(3,3))**2
     2619#             G = nl.inv(G)
     2620#             pFil.close()
     2621#             break
     2622#     dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)]
     2623#     return min(dspaces)
     2624   
     2625def MakefullrmcRun(pName,Phase,RMCPdict):
     2626    BondList = {}
     2627    for k in RMCPdict['Pairs']:
     2628        if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0:
     2629            BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2])
     2630    AngleList = []
     2631    for angle in RMCPdict['Angles']:
     2632        if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0:
    25282633            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 
    2546 def 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 
    2583 def 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 
    2606 def 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    
    2623 def MakefullrmcRun(pName,Phase,RMCPdict):
    2624     Res,BondList = FindBonds(Phase,RMCPdict)
    2625     AngleList = FindAngles(Phase,RMCPdict)
     2634        AngleList.append(angle)
    26262635    rmin = RMCPdict['min Contact']
    2627     rmax = GetMaxSphere(RMCPdict['atomPDB'])
    2628     if rmax is None:
    2629         return None
    2630     rname = pName+'-run.py'
     2636    cell = Phase['General']['Cell'][1:7]
     2637    bigcell = np.array(cell)*np.array(RMCPdict['SuperCell']+[1,1,1])
     2638    bigG = G2lat.cell2Gmat(bigcell)[0]
     2639    rmax = min([0.5/np.sqrt(G2lat.calc_rDsq2(H,bigG)) for H in np.eye(3)])
     2640    SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0]
     2641    cx,ct,cs,cia = Phase['General']['AtomPtrs']
     2642    atomsList = []
     2643    for atom in Phase['Atoms']:
     2644        el = ''.join([i for i in atom[ct] if i.isalpha()])
     2645        atomsList.append([el] + atom[cx:cx+4])
     2646    rname = pName+'-fullrmc.py'
    26312647    restart = '%s_restart.pdb'%pName
    26322648    Files = RMCPdict['files']
     
    26342650    rundata = ''
    26352651    rundata += '#### fullrmc %s file; edit by hand if you so choose #####\n'%rname
     2652    rundata += '# created in '+__file__+" v"+filversion.split()[1]
     2653    rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n")
    26362654    rundata += '''
    26372655# fullrmc imports (all that are potentially useful)
     2656#import matplotlib.pyplot as plt
    26382657import numpy as np
    26392658import time
    2640 from fullrmc.sincConvolution import sincConvolution
    2641 from fullrmc.Globals import LOGGER
     2659from fullrmc.Core import Collection
    26422660from fullrmc.Engine import Engine
    2643 from fullrmc.Constraints.PairDistributionConstraints import PairDistributionConstraint
    2644 from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint
    2645 from fullrmc.Constraints.DistanceConstraints import InterMolecularDistanceConstraint
     2661import fullrmc.Constraints.PairDistributionConstraints as fPDF
     2662from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint
     2663from fullrmc.Constraints.DistanceConstraints import DistanceConstraint
    26462664from fullrmc.Constraints.BondConstraints import BondConstraint
    26472665from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint
    26482666from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint
    26492667from fullrmc.Generators.Swaps import SwapPositionsGenerator
    2650 from fullrmc.debugStuff import *
    2651 InvokeDebugOpts()
     2668
     2669### When True, erases an existing enging to provide a fresh start
     2670FRESH_START = {}
    26522671time0 = time.time()
    2653 SwapGen = {}
    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
     2672'''.format(RMCPdict['ReStart'][0])
     2673   
     2674    rundata += '# setup structure\n'
     2675    rundata += 'cell = ' + str(cell) + '\n'
     2676    rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n'
     2677    rundata += 'atomList = ' + str(atomsList).replace('],','],\n  ') + '\n'
     2678    rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n'
     2679
     2680    rundata += '\n# initialize engine\n'
    26672681    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'
     2682
     2683    rundata += '''\n# check Engine exists if so (and not FRESH_START) load it
     2684# otherwise build it
     2685ENGINE = Engine(path=None)
     2686if not ENGINE.is_engine(engineFileName) or FRESH_START:
     2687    ## create structure
     2688    ENGINE = Engine(path=engineFileName, freshStart=True)
     2689    ENGINE.build_crystal_set_pdb(symOps     = SymOpList,
     2690                                 atoms      = atomList,
     2691                                 unitcellBC = cell,
     2692                                 supercell  = supercell)
     2693    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
     2694'''   
     2695    import atmdata
     2696    rundata += '# conversion factors (may be needed)\n'
     2697    rundata += '    sumCiBi2 = 0.\n'
     2698    for elem in Phase['General']['AtomTypes']:
     2699        rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
     2700        rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
     2701    rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
     2702    # settings that require a new Engine
    26742703    for File in Files:
    26752704        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]
     2705        if not os.path.exists(filDat[0]): continue
     2706        sfwt = 'neutronCohb'
     2707        if 'Xray' in File:
     2708            sfwt = 'atomicNumber'
     2709        if 'G(r)' in File:
     2710            rundata += '    GR = np.loadtxt("%s").T\n'%filDat[0]
     2711            if filDat[3] == 0:
     2712                rundata += '''    # read and xform G(r) as defined in RMCProfile
     2713# see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
     2714                rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
     2715                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
     2716            elif filDat[3] == 1:
     2717                rundata += '    # This is G(r) as defined in PDFFIT\n'
     2718                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
     2719            elif filDat[3] == 2:
     2720                rundata += '    # This is g(r)\n'
     2721                rundata += '    GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
    26872722            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)
     2723                raise ValueError('Invalid G(r) type: '+str(filDat[3]))
     2724            rundata += '    ENGINE.add_constraints([GofR])\n'
     2725            rundata += '    GofR.set_limits((None, rmax))\n'
     2726            wtDict['Pair-'+sfwt] = filDat[1]
     2727        elif 'F(Q)' in File:
     2728            rundata += '    SOQ = np.loadtxt("%s").T\n'%filDat[0]
     2729            if filDat[3] == 0:
     2730                rundata += '    # Read & xform F(Q) as defined in RMCProfile\n'
     2731                rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
     2732                rundata += '    SOQ[1] += 1\n'
     2733            elif filDat[3] == 1:
     2734                rundata += '    # This is S(Q) as defined in PDFFIT\n'
     2735            if filDat[4]:
     2736                rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax={:.3f})\n'.format(rmax)
     2737            rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
     2738            rundata += '    ENGINE.add_constraints([SofQ])\n'
     2739        else:
     2740            print('What is this?')
     2741        wtDict['Struct-'+sfwt] = filDat[1]
     2742    rundata += '    ENGINE.add_constraints(DistanceConstraint(defaultLowerDistance={}))\n'.format(RMCPdict['min Contact'])
     2743    rundata += '''    B_CONSTRAINT   = BondConstraint()
     2744    ENGINE.add_constraints(B_CONSTRAINT)
     2745    B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[
     2746'''
     2747    for pair in BondList:
     2748        e1,e2 = pair.split('-')
     2749        rundata += '            ("element","{}","{}",{},{}),\n'.format(
     2750                                    e1.strip(),e2.strip(),*BondList[pair])
     2751    rundata += '''             ])
     2752    ENGINE.save()
     2753else:
     2754    ENGINE = ENGINE.load(path=engineFileName)
     2755'''               
     2756#    if RMCPdict['Swaps']:
     2757#        rundata += '#set up for site swaps\n'
     2758#        rundata += 'aN = ENGINE.allNames\n'
     2759#        rundata += 'SwapGen = {}\n'
     2760#        for swap in RMCPdict['Swaps']:
     2761#            rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0]
     2762#            rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1]
     2763#            rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2])
     2764    rundata += '# setup/change constraints - can be done without restart\n'   
     2765    rundata += 'for c in ENGINE.constraints:  # loop over predefined constraints\n'
     2766    rundata += '    strcons = str(type(c))\n'
     2767    rundata += '    if type(c) is fPDF.PairDistributionConstraint:\n'
     2768    rundata += '        c.set_variance_squared(wtDict["Pair-"+c.weighting])\n'
     2769    rundata += '        c.set_limits((None,%.3f))\n'%(rmax)
    27192770    if RMCPdict['FitScale']:
    2720         rundata += '        constraint.set_adjust_scale_factor((10, 0.01, 100.))\n'
     2771        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
    27212772    rundata += '    elif "StructureFactor" in strcons:\n'
    2722     rundata += '        constraint.set_variance_squared(wtDict["Struct-"+constraint.weighting])\n'
     2773    rundata += '        c.set_variance_squared(wtDict["Struct-"+c.weighting])\n'
    27232774    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'
     2775        rundata += '        c.set_adjust_scale_factor((10, 0.01, 100.))\n'
     2776#    if AngleList and not RMCPdict['Swaps']: rundata += setAngleConstraints()
     2777    # if len(RMCPdict['Torsions']):         # Torsions currently commented out in GUI
     2778    #     rundata += 'for c in ENGINE.constraints:  # look for Dihedral Angle Constraints\n'
     2779    #     rundata += '    if type(c) is DihedralAngleConstraint:\n'
     2780    #     rundata += '        c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight']
     2781    #     rundata += '        c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res
     2782    #     for torsion in RMCPdict['Torsions']:
     2783    #         rundata += '    %s\n'%str(tuple(torsion))
     2784    #     rundata += '        ]})\n'           
     2785    rundata += '\n# setup runs for fullrmc\n'
    27552786
    27562787    rundata += 'for _ in range(%d):\n'%RMCPdict['Cycles']
     2788    if BondList and RMCPdict['Swaps']: rundata += setBondConstraints('    ')
     2789    if AngleList and RMCPdict['Swaps']: rundata += setAngleConstraints('    ')
    27572790    rundata += '    ENGINE.set_groups_as_atoms()\n'
    27582791    rundata += '    ENGINE.run(restartPdb="%s",numberOfSteps=10000, saveFrequency=1000)\n'%restart
    2759     if len(RMCPdict['Swaps']):
     2792    if RMCPdict['Swaps']:
     2793        if BondList: rundata += setBondConstraints('    ',clear=True)
     2794        if AngleList: rundata += setAngleConstraints('    ',clear=True)
    27602795        rundata += '    for swaps in SwapGen:\n'
    27612796        rundata += '        AB = swaps.split("-")\n'
    27622797        rundata += '        ENGINE.set_groups_as_atoms()\n'
    27632798        rundata += '        for g in ENGINE.groups:\n'
    2764         rundata += '            if allNames[g.indexes[0]]==AB[0]:\n'
     2799        rundata += '            if aN[g.indexes[0]]==AB[0]:\n'
    27652800        rundata += '                g.set_move_generator(SwapGen[swaps][0])\n'
    2766         rundata += '            elif allNames[g.indexes[0]]==AB[1]:\n'
     2801        rundata += '            elif aN[g.indexes[0]]==AB[1]:\n'
    27672802        rundata += '                g.set_move_generator(SwapGen[swaps][1])\n'
    27682803        rundata += '            sProb = SwapGen[swaps][2]\n'
     
    27702805        rundata += '        ENGINE.set_groups_as_atoms()\n'
    27712806        rundata += '        ENGINE.run(restartPdb="%s",numberOfSteps=10000*(1.-sProb), saveFrequency=1000)\n'%restart
    2772     rundata += 'ENGINE.close()\n'
     2807    #rundata += 'ENGINE.close()\n'
    27732808    rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n'
    27742809    rfile = open(rname,'w')
     
    27782813    return rname
    27792814
    2780 def MakefullrmcPDB(Name,Phase,RMCPdict):
    2781     generalData = Phase['General']
    2782     Atseq = RMCPdict['atSeq']
    2783     Dups,Fracs = findDup(Phase['Atoms'])
    2784     Sfracs = [np.cumsum(fracs) for fracs in Fracs]
    2785     ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
    2786     Supercell = RMCPdict['SuperCell']
    2787     Cell = generalData['Cell'][1:7]
    2788     Trans = np.eye(3)*np.array(Supercell)
    2789     newPhase = copy.deepcopy(Phase)
    2790     newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
    2791     newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
    2792     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
    2793     Atoms = newPhase['Atoms']
    2794 
    2795     if ifSfracs:
    2796         Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
    2797         Natm = np.count_nonzero(Natm-1)
    2798         Satoms = []
    2799         for i in range(len(Atoms)//Natm):
    2800             ind = i*Natm
    2801             Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
    2802         Natoms = []
    2803         for satoms in Satoms:
    2804             for idup,dup in enumerate(Dups):
    2805                 ldup = len(dup)
    2806                 natm = len(satoms)
    2807                 i = 0
    2808                 while i < natm:
    2809                     if satoms[i][0] in dup:
    2810                         atoms = satoms[i:i+ldup]
    2811                         try:
    2812                             atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
    2813                             Natoms.append(atom)
    2814                         except IndexError:      #what about vacancies?
    2815                             if 'Va' not in Atseq:
    2816                                 Atseq.append('Va')
    2817                                 RMCPdict['aTypes']['Va'] = 0.0
    2818                             atom = atoms[0]
    2819                             atom[1] = 'Va'
    2820                             Natoms.append(atom)
    2821                         i += ldup
    2822                     else:
    2823                        i += 1
    2824     else:
    2825         Natoms = Atoms
    2826 
    2827     XYZ = np.array([atom[3:6] for atom in Natoms]).T
    2828     XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
    2829     Cell = newPhase['General']['Cell'][1:7]
    2830     A,B = G2lat. cell2AB(Cell)
    2831     fname = Name+'_cbb.pdb'
    2832     fl = open(fname,'w')
    2833     fl.write('REMARK    Boundary Conditions:%6.2f  0.0  0.0  0.0%7.2f  0.0  0.0  0.0%7.2f\n'%(
    2834              Cell[0],Cell[1],Cell[2]))
    2835     fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
    2836     fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
    2837     fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
    2838     fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
    2839             Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
    2840 
    2841     Natm = np.core.defchararray.count(np.array(Atcodes),'+')
    2842     Natm = np.count_nonzero(Natm-1)
    2843     nat = 0
    2844     if RMCPdict['byMolec']:
    2845         NPM = RMCPdict['Natoms']
    2846         for iat,atom in enumerate(Natoms):
    2847             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
    2848             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    2849                     1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    2850             nat += 1
    2851     else:
    2852         for atm in Atseq:
    2853             for iat,atom in enumerate(Natoms):
    2854                 if atom[1] == atm:
    2855                     XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
    2856                     fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
    2857                             1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
    2858                     nat += 1
    2859     fl.close()
    2860     return fname
     2815# def MakefullrmcPDB(Name,Phase,RMCPdict):
     2816#     generalData = Phase['General']
     2817#     Atseq = RMCPdict['atSeq']
     2818#     Dups,Fracs = findDup(Phase['Atoms'])
     2819#     Sfracs = [np.cumsum(fracs) for fracs in Fracs]
     2820#     ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs])
     2821#     Supercell = RMCPdict['SuperCell']
     2822#     Cell = generalData['Cell'][1:7]
     2823#     Trans = np.eye(3)*np.array(Supercell)
     2824#     newPhase = copy.deepcopy(Phase)
     2825#     newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1]
     2826#     newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans.T)
     2827#     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True)
     2828#     Atoms = newPhase['Atoms']
     2829
     2830#     if ifSfracs:
     2831#         Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
     2832#         Natm = np.count_nonzero(Natm-1)
     2833#         Satoms = []
     2834#         for i in range(len(Atoms)//Natm):
     2835#             ind = i*Natm
     2836#             Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3))
     2837#         Natoms = []
     2838#         for satoms in Satoms:
     2839#             for idup,dup in enumerate(Dups):
     2840#                 ldup = len(dup)
     2841#                 natm = len(satoms)
     2842#                 i = 0
     2843#                 while i < natm:
     2844#                     if satoms[i][0] in dup:
     2845#                         atoms = satoms[i:i+ldup]
     2846#                         try:
     2847#                             atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
     2848#                             Natoms.append(atom)
     2849#                         except IndexError:      #what about vacancies?
     2850#                             if 'Va' not in Atseq:
     2851#                                 Atseq.append('Va')
     2852#                                 RMCPdict['aTypes']['Va'] = 0.0
     2853#                             atom = atoms[0]
     2854#                             atom[1] = 'Va'
     2855#                             Natoms.append(atom)
     2856#                         i += ldup
     2857#                     else:
     2858#                        i += 1
     2859#     else:
     2860#         Natoms = Atoms
     2861
     2862#     XYZ = np.array([atom[3:6] for atom in Natoms]).T
     2863#     XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
     2864#     Cell = newPhase['General']['Cell'][1:7]
     2865#     A,B = G2lat. cell2AB(Cell)
     2866#     fname = Name+'_cbb.pdb'
     2867#     fl = open(fname,'w')
     2868#     fl.write('REMARK    Boundary Conditions:%6.2f  0.0  0.0  0.0%7.2f  0.0  0.0  0.0%7.2f\n'%(
     2869#              Cell[0],Cell[1],Cell[2]))
     2870#     fl.write('ORIGX1      1.000000  0.000000  0.000000        0.00000\n')
     2871#     fl.write('ORIGX2      0.000000  1.000000  0.000000        0.00000\n')
     2872#     fl.write('ORIGX3      0.000000  0.000000  1.000000        0.00000\n')
     2873#     fl.write('CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n'%(
     2874#             Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5]))
     2875
     2876#     Natm = np.core.defchararray.count(np.array(Atcodes),'+')
     2877#     Natm = np.count_nonzero(Natm-1)
     2878#     nat = 0
     2879#     if RMCPdict['byMolec']:
     2880#         NPM = RMCPdict['Natoms']
     2881#         for iat,atom in enumerate(Natoms):
     2882#             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
     2883#             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
     2884#                     1+nat%NPM,atom[0],1+nat//NPM,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
     2885#             nat += 1
     2886#     else:
     2887#         for atm in Atseq:
     2888#             for iat,atom in enumerate(Natoms):
     2889#                 if atom[1] == atm:
     2890#                     XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
     2891#                     fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
     2892#                             1+nat,atom[0],1+nat,XYZ[0],XYZ[1],XYZ[2],atom[1].lower()))
     2893#                     nat += 1
     2894#     fl.close()
     2895#     return fname
    28612896   
    28622897def MakePdparse(RMCPdict):
     
    32523287    Reflectometry as a function of kz for a set of slabs.
    32533288
    3254     :param kz: float[n] (1/Ang). Scattering vector, :math:`2\pi\sin(\\theta)/\lambda`.
     3289    :param kz: float[n] (1/Ang). Scattering vector, :math:`2\\pi\\sin(\\theta)/\\lambda`.
    32553290        This is :math:`\\tfrac12 Q_z`.       
    32563291    :param depth:  float[m] (Ang).
Note: See TracChangeset for help on using the changeset viewer.