Changeset 4518 for trunk/GSASIIpwd.py
- Timestamp:
- Jul 12, 2020 8:59:57 PM (3 years ago)
- File:
-
- 1 edited
-
trunk/GSASIIpwd.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIpwd.py
r4483 r4518 20 20 import os.path 21 21 import subprocess as subp 22 import datetime as dt 22 23 import copy 23 24 … … 33 34 34 35 import GSASIIpath 36 filversion = "$Revision$" 35 37 GSASIIpath.SetVersionNumber("$Revision$") 36 38 import GSASIIlattice as G2lat … … 2511 2513 return fname 2512 2514 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 2625 def 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: 2528 2633 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) 2626 2635 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' 2631 2647 restart = '%s_restart.pdb'%pName 2632 2648 Files = RMCPdict['files'] … … 2634 2650 rundata = '' 2635 2651 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") 2636 2654 rundata += ''' 2637 2655 # fullrmc imports (all that are potentially useful) 2656 #import matplotlib.pyplot as plt 2638 2657 import numpy as np 2639 2658 import time 2640 from fullrmc.sincConvolution import sincConvolution 2641 from fullrmc.Globals import LOGGER 2659 from fullrmc.Core import Collection 2642 2660 from fullrmc.Engine import Engine 2643 from fullrmc.Constraints.PairDistributionConstraints import PairDistributionConstraint 2644 from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint 2645 from fullrmc.Constraints.DistanceConstraints import InterMolecularDistanceConstraint2661 import fullrmc.Constraints.PairDistributionConstraints as fPDF 2662 from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint 2663 from fullrmc.Constraints.DistanceConstraints import DistanceConstraint 2646 2664 from fullrmc.Constraints.BondConstraints import BondConstraint 2647 2665 from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint 2648 2666 from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint 2649 2667 from 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 2670 FRESH_START = {} 2652 2671 time0 = 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' 2667 2681 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 2685 ENGINE = Engine(path=None) 2686 if 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 2674 2703 for File in Files: 2675 2704 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 2687 2722 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() 2753 else: 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) 2719 2770 if RMCPdict['FitScale']: 2720 rundata += ' c onstraint.set_adjust_scale_factor((10, 0.01, 100.))\n'2771 rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n' 2721 2772 rundata += ' elif "StructureFactor" in strcons:\n' 2722 rundata += ' c onstraint.set_variance_squared(wtDict["Struct-"+constraint.weighting])\n'2773 rundata += ' c.set_variance_squared(wtDict["Struct-"+c.weighting])\n' 2723 2774 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' 2755 2786 2756 2787 rundata += 'for _ in range(%d):\n'%RMCPdict['Cycles'] 2788 if BondList and RMCPdict['Swaps']: rundata += setBondConstraints(' ') 2789 if AngleList and RMCPdict['Swaps']: rundata += setAngleConstraints(' ') 2757 2790 rundata += ' ENGINE.set_groups_as_atoms()\n' 2758 2791 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) 2760 2795 rundata += ' for swaps in SwapGen:\n' 2761 2796 rundata += ' AB = swaps.split("-")\n' 2762 2797 rundata += ' ENGINE.set_groups_as_atoms()\n' 2763 2798 rundata += ' for g in ENGINE.groups:\n' 2764 rundata += ' if a llNames[g.indexes[0]]==AB[0]:\n'2799 rundata += ' if aN[g.indexes[0]]==AB[0]:\n' 2765 2800 rundata += ' g.set_move_generator(SwapGen[swaps][0])\n' 2766 rundata += ' elif a llNames[g.indexes[0]]==AB[1]:\n'2801 rundata += ' elif aN[g.indexes[0]]==AB[1]:\n' 2767 2802 rundata += ' g.set_move_generator(SwapGen[swaps][1])\n' 2768 2803 rundata += ' sProb = SwapGen[swaps][2]\n' … … 2770 2805 rundata += ' ENGINE.set_groups_as_atoms()\n' 2771 2806 rundata += ' ENGINE.run(restartPdb="%s",numberOfSteps=10000*(1.-sProb), saveFrequency=1000)\n'%restart 2772 rundata += 'ENGINE.close()\n'2807 #rundata += 'ENGINE.close()\n' 2773 2808 rundata += 'print("ENGINE run time %.2f s"%(time.time()-time0))\n' 2774 2809 rfile = open(rname,'w') … … 2778 2813 return rname 2779 2814 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 cell2797 Natm = np.count_nonzero(Natm-1)2798 Satoms = []2799 for i in range(len(Atoms)//Natm):2800 ind = i*Natm2801 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 = 02808 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.02818 atom = atoms[0]2819 atom[1] = 'Va'2820 Natoms.append(atom)2821 i += ldup2822 else:2823 i += 12824 else:2825 Natoms = Atoms2826 2827 XYZ = np.array([atom[3:6] for atom in Natoms]).T2828 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 = 02844 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 += 12851 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 Cartesian2856 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 += 12859 fl.close()2860 return fname2815 # 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 2861 2896 2862 2897 def MakePdparse(RMCPdict): … … 3252 3287 Reflectometry as a function of kz for a set of slabs. 3253 3288 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`. 3255 3290 This is :math:`\\tfrac12 Q_z`. 3256 3291 :param depth: float[m] (Ang).
Note: See TracChangeset
for help on using the changeset viewer.
