Changeset 4415 for trunk


Ignore:
Timestamp:
May 7, 2020 12:17:26 PM (3 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.

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIimage.py

    r4350 r4415  
    571571    dist = data['distance']/npcosd(tilt)
    572572    x0 = data['distance']*nptand(tilt)
     573    x0x = x0*npcosd(data['rotation'])
     574    x0y = x0*npsind(data['rotation'])
    573575    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
    574576    distsq = data['distance']**2
    575577    dx = x-data['center'][0]
    576578    dy = y-data['center'][1]
    577     G = ((dx-x0)**2+dy**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
     579    G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
    578580    Z = np.dot(np.dstack([dx.T,dy.T,np.zeros_like(dx.T)]),MN).T[2]
    579581    xyZ = dx**2+dy**2-Z**2   

  • trunk/GSASIIlattice.py

    r4268 r4415  
    298298    return newAtoms,Codes
    299299   
    300 def TransformPhase(oldPhase,newPhase,Trans,Uvec,Vvec,ifMag):
     300def TransformPhase(oldPhase,newPhase,Trans,Uvec,Vvec,ifMag,Force=True):
    301301    '''Transform atoms from oldPhase to newPhase
    302302    M' is inv(M)
    … …  
    325325    SGData = newPhase['General']['SGData']
    326326    invTrans = nl.inv(Trans)
    327     newAtoms,atCodes = FillUnitCell(oldPhase)
     327    newAtoms,atCodes = FillUnitCell(oldPhase,Force)
    328328    newAtoms,atCodes = ExpandCell(newAtoms,atCodes,cx,Trans)
    329329    if ifMag:
    … …  
    343343        newPhase['Draw Atoms'] = []
    344344    for atom in newAtoms:
    345         xyz = TransformXYZ(atom[cx:cx+3]+Uvec,invTrans.T,Vvec)%1.
    346         atom[cx:cx+3] = np.around(xyz,6)%1.
     345        xyz = TransformXYZ(atom[cx:cx+3]+Uvec,invTrans.T,Vvec)
     346        if Force:
     347            xyz = np.around(xyz,6)%1.
     348        atom[cx:cx+3] = xyz
    347349        if atom[cia] == 'A':
    348350            atom[cia+2:cia+8] = TransformU6(atom[cia+2:cia+8],Trans)
    … …  
    476478    return phase
    477479
    478 def FillUnitCell(Phase):
     480def FillUnitCell(Phase,Force=True):
    479481    Atoms = copy.deepcopy(Phase['Atoms'])
    480482    atomData = []
    … …  
    489491    for iat,atom in enumerate(Atoms):
    490492        XYZ = np.array(atom[cx:cx+3])
    491         xyz = XYZ%1.
     493        xyz = XYZ
     494        if Force:
     495            xyz %= 1.
    492496        if atom[cia] == 'A':
    493497            Uij = atom[cia+2:cia+8]
    494             result = G2spc.GenAtom(xyz,SGData,False,Uij,True)
     498            result = G2spc.GenAtom(xyz,SGData,False,Uij,False)
    495499            for item in result:
    496                 if item[0][2] >= .95: item[0][2] -= 1.
     500#                if item[0][2] >= .95: item[0][2] -= 1.
    497501                atom[cx:cx+3] = item[0]
    498502                atom[cia+2:cia+8] = item[1]

  • trunk/GSASIImath.py

    r4401 r4415  
    492492   
    493493    :param list atomData: Atom table to be used
     494    :param int  indx: pointer to position of atom id in atom record (typically cia+8)
    494495   
    495496    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
    496497   
    497498    '''
    498     atomLookUp = {}
    499     for iatm,atom in enumerate(atomData):
    500         atomLookUp[atom[indx]] = iatm
    501     return atomLookUp
     499    return {atom[indx]:iatm for iatm,atom in enumerate(atomData)}
    502500
    503501def GetAtomsById(atomData,atomLookUp,IdList):
    … …  
    636634    for j in IndB[0]:
    637635        if j != Orig:
    638             if AtNames[j] != notName:
     636            if AtNames[j] not in notName:
    639637                Neigh.append([AtNames[j],dist[j],True])
    640638                Ids.append(Atoms[j][cia+8])
    … …  
    697695    return bond,std,meanDisp,stdDisp,A,V,vecDisp
    698696   
    699 def FindAllNeighbors(phase,FrstName,AtNames,notName=''):
     697def FindAllNeighbors(phase,FrstName,AtNames,notName='',Orig=None,Short=False):
    700698    General = phase['General']
    701699    cx,ct,cs,cia = General['AtomPtrs']
    … …  
    713711    radiusFactor = DisAglCtls['Factors'][0]
    714712    AtInfo = dict(zip(AtTypes,Radii)) #or General['BondRadii']
    715     Orig = atNames.index(FrstName)
     713    if Orig is None:
     714        Orig = atNames.index(FrstName)
    716715    OId = Atoms[Orig][cia+8]
    717716    OType = Atoms[Orig][ct]
    … …  
    739738                        else:
    740739                            Topstr = ' +(%4d)'%(Top)
    741                         Neigh.append([AtNames[iA]+Topstr,atTypes[iA],dist[iU],dx[iU]])
     740                        if Short:
     741                            Neigh.append([AtNames[iA],dist[iU],True])
     742                        else:
     743                            Neigh.append([AtNames[iA]+Topstr,atTypes[iA],dist[iU],dx[iU]])
    742744                        Ids.append(Atoms[iA][cia+8])
    743745    return Neigh,[OId,Ids]
    … …  
    19521954        vec[i] /= dist
    19531955    angle = acosd(np.sum(vec[0]*vec[1]))
    1954 #    GSASIIpath.IPyBreak()
    19551956    return angle
    19561957

  • trunk/GSASIIphsGUI.py

    r4405 r4415  
    6464import numpy as np
    6565import numpy.linalg as nl
     66import numpy.fft as fft
    6667import atmdata
    6768
    … …  
    27702771        nTet = 0
    27712772        nElse = 0
    2772         for atom in data['Atoms']:
     2773        for iatm,atom in enumerate(data['Atoms']):
    27732774            if atom[ct] == Oatoms:
    2774                 results = G2mth.FindAllNeighbors(data,atom[ct-1],atNames)[0]      #slow step
     2775                results = G2mth.FindAllNeighbors(data,atom[ct-1],atNames,Orig=iatm)[0]      #slow step
    27752776                if len(results) == 4:
    27762777                    bond,std,meanDisp,stdDisp,A,V,dVec = G2mth.FindTetrahedron(results)
    … …  
    30063007                    dlg.Destroy()
    30073008                elif Atoms.GetColLabelValue(c) == 'Uiso':       #this needs to ask for value
    3008                     pass                                        #& then change all 'I' atoms
     3009                    return                                        #& then change all 'I' atoms; now do nothing
    30093010                else:
    30103011                    return
    … …  
    39773978        indx = GetSelectedAtoms()
    39783979        DisAglCtls = {}
    3979         if len(indx) == 1:
     3980        if indx is not None and len(indx) == 1:
    39803981            generalData = data['General']
    39813982            if 'DisAglCtls' in generalData:
    … …  
    45034504                        start += 1
    45044505                Xlab = 'Q'
    4505                 if 'G(R)' in fileItem[2]:
     4506                if 'G(R)' in fileItem[2].upper():
    45064507                    Xlab = 'R'
    45074508                G2plt.PlotXY(G2frame,[XY.T,],labelX=Xlab,
    45084509                    labelY=fileItem[2],newPlot=True,Title=fileItem[0],
    45094510                    lines=True)
     4511               
     4512            def OnCorrChk(event):
     4513                Obj = event.GetEventObject()
     4514                fil = Indx[Obj.GetId()]
     4515                RMCPdict['files'][fil][3] = not RMCPdict['files'][fil][3]
    45104516                           
    45114517            Indx = {}
    … …  
    45224528                mainSizer.Add(wx.StaticText(G2frame.FRMC,
    45234529                    label=' NB: fullrmc data files must be 2 columns; all other lines preceeded by "#". Edit before use.'),0,WACV)
    4524                 ncol = 4
    4525                 Heads = ['Name','File','Weight','Plot']
     4530                Heads = ['Name','File','Weight','Plot','Corr']
    45264531            fileSizer = wx.FlexGridSizer(ncol,5,5)
    45274532            Formats = ['RMC','GUDRUN','STOG']
    … …  
    45494554                    Indx[plotBtn.GetId()] = fil
    45504555                    fileSizer.Add(plotBtn,0,WACV)
     4556                    if G2frame.RMCchoice == 'fullrmc':
     4557                        lab = 'Apply sinc convolution? '
     4558                        if 'G(r)' in fil:
     4559                            lab = 'Apply r*G(r) conversion? '
     4560                        corrChk = wx.CheckBox(G2frame.FRMC,label=lab)
     4561                        #patch
     4562                        if len(RMCPdict['files'][fil]) < 4:
     4563                            RMCPdict['files'][fil].append(True)
     4564                        #end patch
     4565                        corrChk.SetValue(RMCPdict['files'][fil][3])
     4566                        Indx[corrChk.GetId()] = fil
     4567                        corrChk.Bind(wx.EVT_CHECKBOX,OnCorrChk)
     4568                        fileSizer.Add(corrChk,0,WACV)
    45514569                else:
    45524570                    RMCPdict['files'][fil][0] = 'Select'
    45534571                    fileSizer.Add((5,5),0)
    45544572                    fileSizer.Add((5,5),0)
    4555                     if G2frame.RMCchoice == 'RMCProfile':                       
    4556                         fileSizer.Add((5,5),0)
     4573                    fileSizer.Add((5,5),0)
    45574574            return fileSizer
    45584575           
    … …  
    46014618                    Pairs += pair
    46024619                Pairs = {pairs:[0.0,0.0,0.0] for pairs in Pairs}
    4603                 files = {'Neutron real space data; G(r): ':['Select',0.05,'G(r)'],
    4604                           'Neutron reciprocal space data; F(Q): ':['Select',0.05,'F(Q)'],
    4605                           'Neutron reciprocal space data; S(Q): ':['Select',0.05,'S(Q)'],
    4606                           'Xray real space data; G(r): ':['Select',0.01,'G(r)'],
    4607                           'Xray reciprocal space data; F(Q): ':['Select',0.01,'F(Q)'],}
     4620                files = {'Neutron real space data; G(r): ':['Select',0.05,'G(r)',True],
     4621                          'Neutron reciprocal space data; F(Q): ':['Select',0.05,'F(Q)',True],
     4622                          'Xray real space data; G(r): ':['Select',0.01,'G(r)',True],
     4623                          'Xray reciprocal space data; F(Q): ':['Select',0.01,'F(Q)',True],}
    46084624                data['RMC']['fullrmc'] = {'SuperCell':[1,1,1],'Box':[10.,10.,10.],'aTypes':aTypes,'byMolec':False,
    4609                     'Natoms':1,'atSeq':atSeq,'Pairs':Pairs,'files':files,'ReStart':[False,False],
     4625                    'Natoms':1,'atSeq':atSeq,'Pairs':Pairs,'files':files,'ReStart':[False,False],'Cycles':1,
    46104626                    'Swaps':[],'useBVS':False,'FitScale':False,'AveCN':[],'FxCN':[],'Angles':[],'Angle Weight':1.e-5,
    4611                     'moleculePdb':'Select','targetDensity':1.0,'maxRecursion':10000,
    4612                     'atomPDB':''}
     4627                    'moleculePdb':'Select','targetDensity':1.0,'maxRecursion':10000,'Torsions':[],'Torsion Weight':1.e-5,
     4628                    'atomPDB':'','Bond Weight':1.e-5,'min Contact':1.5}
    46134629            RMCPdict = data['RMC']['fullrmc']
    46144630
    … …  
    46864702                wx.CallAfter(UpdateRMC)
    46874703               
     4704            def OnAddTorsion(event):
     4705                RMCPdict['Torsions'].append(['','','','',0.,0.,0.,0.,0.,0.])
     4706                wx.CallAfter(UpdateRMC)
     4707               
    46884708            def GetAngleSizer():
    46894709               
    … …  
    47224742                    angleSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,angle,4,min=0.,max=180.,OnLeave=SetRestart1,size=(50,25)),0,WACV)
    47234743                return angleSizer
     4744           
     4745            def GetTorsionSizer():
     4746               
     4747                def OnDelTorsion(event):
     4748                    Obj = event.GetEventObject()
     4749                    angle = Indx[Obj.GetId()]
     4750                    del RMCPdict['Torsions'][angle]
     4751                    wx.CallAfter(UpdateRMC)
     4752                   
     4753                def OnTorsionAtSel(event):
     4754                    Obj = event.GetEventObject()
     4755                    torsion,i = Indx[Obj.GetId()]
     4756                    RMCPdict['Torsions'][torsion][i] = Obj.GetStringSelection()
     4757                                           
     4758                def SetRestart1(invalid,value,tc):
     4759                    RMCPdict['ReStart'][1] = True
     4760               
     4761                Indx = {}
     4762                torsionSizer = wx.FlexGridSizer(11,5,5)
     4763                fxcnLabels = [' ','Atom-A','Atom-B','Atom-C','Atom-D',' min angle1',' max angle1',' min angle2',' max angle2',' min angle3',' max angle3']
     4764                for lab in fxcnLabels:
     4765                    torsionSizer.Add(wx.StaticText(G2frame.FRMC,label=lab),0,WACV)
     4766                for ifx,torsion in enumerate(RMCPdict['Torsions']):
     4767                    delBtn = wx.Button(G2frame.FRMC,label='Delete')
     4768                    delBtn.Bind(wx.EVT_BUTTON,OnDelTorsion)
     4769                    Indx[delBtn.GetId()] = ifx
     4770                    torsionSizer.Add(delBtn,0,WACV)
     4771                    for i in [0,1,2,3]:
     4772                        atmSel = wx.ComboBox(G2frame.FRMC,choices=atNames,style=wx.CB_DROPDOWN|wx.TE_READONLY)
     4773                        atmSel.SetStringSelection(torsion[i])
     4774                        atmSel.Bind(wx.EVT_COMBOBOX,OnTorsionAtSel)
     4775                        Indx[atmSel.GetId()] = [ifx,i]
     4776                        torsionSizer.Add(atmSel,0,WACV)
     4777                    for i in  [4,5,6,7,8,9]:
     4778                        torsionSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,torsion,i,min=0.,max=360.,OnLeave=SetRestart1,size=(50,25)),0,WACV)
     4779                return torsionSizer
    47244780#patches
    47254781            if 'useBVS' not in RMCPdict:
    … …  
    47374793            if 'Angles' not in RMCPdict:
    47384794                RMCPdict.update({'Angles':[],'Angle Weight':1.e-5,'Bond Weight':1.e-5,'Torsions':[],'Torsion Weight':1.e-5})
     4795            if 'Cycles' not in RMCPdict:
     4796                RMCPdict['Cycles'] = 1
     4797            if 'min Contact' not in RMCPdict:
     4798                RMCPdict['min Contact'] = 1.5
    47394799#end patches
    47404800
    47414801            generalData = data['General']
     4802            cx,ct,cs,cia = generalData['AtomPtrs']
     4803            atomData = data['Atoms']
     4804            atNames = [atom[ct-1] for atom in atomData]
     4805            ifP1 = False
     4806            if generalData['SGData']['SpGrp'] == 'P 1':
     4807                ifP1 = True               
    47424808            ifBox = False
    47434809            if 'macromolecular' in generalData['Type']:
    … …  
    47474813                lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Big box dimensions, %s:'%Angstr),0,WACV)               
    47484814                lineSizer.Add(GetBoxSizer(),0,WACV)
    4749             else:
     4815            elif ifP1:
    47504816                lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Lattice multipliers:'),0,WACV)
    47514817                lineSizer.Add(GetSuperSizer(),0,WACV)
    … …  
    47564822                lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Num. atoms per molecule '),0,WACV)
    47574823                lineSizer.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Natoms',min=1,size=[40,25]),0,WACV)
     4824            else:
     4825                lineSizer.Add(wx.StaticText(G2frame.FRMC,label=' Starting phase symmetry must be P 1; transform structure first'))
    47584826            mainSizer.Add(lineSizer,0,WACV)
    47594827            if ifBox:
    … …  
    47714839                molecSizer.Add(makePDB,0,WACV)               
    47724840                mainSizer.Add(molecSizer,0,WACV)
    4773             else:
     4841            elif ifP1:
    47744842                makePDB = wx.Button(G2frame.FRMC,label='Make big box PDB')
    47754843                makePDB.Bind(wx.EVT_BUTTON,OnMakePDB)
    4776                 mainSizer.Add(makePDB,0,WACV)               
     4844                mainSizer.Add(makePDB,0,WACV)
     4845            else:       #Abort because starting phase symmetry isn't P 1
     4846                SetPhaseWindow(G2frame.FRMC,mainSizer)
     4847                return
    47774848               
    47784849            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    47794850            mainSizer.Add(wx.StaticText(G2frame.FRMC,label=' fullrmc run file preparation:'),0,WACV)
     4851            resLine = wx.BoxSizer(wx.HORIZONTAL)
    47804852            restart = wx.CheckBox(G2frame.FRMC,label=' Restart fullrmc Engine? (will clear old result!) ')
    47814853            restart.SetValue(RMCPdict['ReStart'][0])
    47824854            restart.Bind(wx.EVT_CHECKBOX,OnReStart)
    4783             mainSizer.Add(restart,0,WACV)
     4855            resLine.Add(restart,0,WACV)
     4856            resLine.Add(wx.StaticText(G2frame.FRMC,label=' Computation cycles: '),0,WACV)
     4857            resLine.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Cycles',min=1,size=[60,25]),0,WACV)
     4858            mainSizer.Add(resLine,0,WACV)
    47844859               
    47854860            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    … …  
    48014876            distBox.Add(wx.StaticText(G2frame.FRMC,label=' Distance constraints, weight: :'),0,WACV)       
    48024877            distBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Bond Weight',min=0.,max=100.,size=(50,25)),0,WACV)
     4878            distBox.Add(wx.StaticText(G2frame.FRMC,label=' min contact dist: '),0,WACV)
     4879            distBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'min Contact',min=0.,max=4.,size=(50,25)),0,WACV)           
    48034880            mainSizer.Add(distBox,0,WACV)
    48044881            mainSizer.Add(GetPairSizer(RMCPdict),0,WACV)
    … …  
    48144891                mainSizer.Add(GetAngleSizer(),0,WACV)
    48154892               
     4893            torBox = wx.BoxSizer(wx.HORIZONTAL)
     4894            torAdd = wx.Button(G2frame.FRMC,label='Add')
     4895            torAdd.Bind(wx.EVT_BUTTON,OnAddTorsion)
     4896            torBox.Add(torAdd,0,WACV)
     4897            torBox.Add(wx.StaticText(G2frame.FRMC,label=' A-B-C-D torsion angle restraints, weight: '),0,WACV)
     4898            torBox.Add(G2G.ValidatedTxtCtrl(G2frame.FRMC,RMCPdict,'Torsion Weight',min=0.,max=100.,size=(50,25)),0,WACV)
     4899            mainSizer.Add(torBox,0,WACV)
     4900            if len(RMCPdict['Torsions']):
     4901                mainSizer.Add(GetTorsionSizer(),0,WACV)
     4902
    48164903            G2G.HorizontalLine(mainSizer,G2frame.FRMC)
    48174904            mainSizer.Add(FileSizer(RMCPdict),0,WACV)
    … …  
    52205307            G2frame.OnFileSaveas(event)
    52215308        if G2frame.RMCchoice == 'fullrmc':
     5309            DisAglCtls = {}
     5310            if 'DisAglCtls' in generalData:
     5311                DisAglCtls = generalData['DisAglCtls']
     5312            dlg = G2G.DisAglDialog(G2frame,DisAglCtls,generalData,Reset=False)
     5313            if dlg.ShowModal() == wx.ID_OK:
     5314                DisAglCtls = dlg.GetData()
     5315                if 'H' not in DisAglCtls['AtomTypes']:
     5316                    DisAglCtls['AtomTypes'].append('H')
     5317                    DisAglCtls['AngleRadii'].append(0.5)
     5318                    DisAglCtls['BondRadii'].append(0.5)
     5319            dlg.Destroy()
     5320            generalData['DisAglCtls'] = DisAglCtls
    52225321            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    52235322            RMCPdict = data['RMC']['fullrmc']
    52245323            fname = G2pwd.MakefullrmcRun(pName,data,RMCPdict)
    5225             print('fullrmc file %s build completed'%fname)
     5324            if fname is None:
     5325                wx.MessageDialog(G2frame,' Big box pdb file is missing; fullrmc will not run','Missing pdb file',wx.OK).ShowModal()
     5326            else:   
     5327                print('fullrmc file %s build completed'%fname)
    52265328        else:
    52275329            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    … …  
    52905392# TBD - remove filedialog & use defined run.py file name here
    52915393            rname = pName+'-run.py'
    5292             dlg = wx.FileDialog(G2frame, 'Choose fullrmc python file to execute', G2G.GetImportPath(G2frame),
    5293                 wildcard='fullrmc python file (*.py)|*.py',style=wx.FD_CHANGE_DIR)
    5294             try:
    5295                 if dlg.ShowModal() == wx.ID_OK:
    5296                     rname = dlg.GetPath()
    5297                 else:
    5298                     return
    5299             finally:
    5300                 dlg.Destroy()
     5394            # dlg = wx.FileDialog(G2frame, 'Choose fullrmc python file to execute', G2G.GetImportPath(G2frame),
     5395            #     wildcard='fullrmc python file (*.py)|*.py',style=wx.FD_CHANGE_DIR)
     5396            # try:
     5397            #     if dlg.ShowModal() == wx.ID_OK:
     5398            #         rname = dlg.GetPath()
     5399            #     else:
     5400            #         return
     5401            # finally:
     5402            #     dlg.Destroy()
    53015403            import subprocess as sb
    53025404            batch = open('fullrmc.bat','w')
    … …  
    54135515                                    if 'total' not in rdfDict:
    54145516                                        print('No data yet - wait for a save')
     5517                                        ENGINE.close()
    54155518                                        if Proc is not None:
    54165519                                            Proc.resume()
    … …  
    54665569                                            G2plt.PlotBarGraph(G2frame,bondAngs,Xname=r'%s Angle, deg'%Aname,Title='%s Bond angles for %s'%(Aname,pName),
    54675570                                                PlotName='%s Angles for %s'%(Aname,pName),maxBins=20)
    5468                                     elif 'ImproperAngleConstraint' in sitem:
     5571                                    elif 'DihedralAngleConstraint' in sitem:
    54695572                                        impangles = 180.*item.get_constraint_value()['angles']/np.pi
    54705573                                        impangleList = item.anglesList[:4]
    … …  
    54765579                                        for Aname in impangleSet:
    54775580                                            impAngs = [angle[1] for angle in impAngles if angle[0]==Aname]
    5478                                             G2plt.PlotBarGraph(G2frame,impAngs,Xname=r'%s $Improper Angle, deg$'%Aname,Title='%s Improper angles for %s'%(Aname,pName),
    5479                                                 PlotName='%s Improper Angles for %s'%(Aname,pName),maxBins=20)
     5581                                            G2plt.PlotBarGraph(G2frame,impAngs,Xname=r'%s $Dihedral Angle, deg$'%Aname,Title='%s Dihedral angles for %s'%(Aname,pName),
     5582                                                PlotName='%s Dihedral Angles for %s'%(Aname,pName),maxBins=20)
    54805583                                    elif 'AtomicCoordinationNumber' in sitem or 'InterMolecular' in sitem:
    54815584                                        print(sitem+' not plotted')

  • trunk/GSASIIplot.py

    r4408 r4415  
    290290    fil.write(line+'\n')
    291291
    292 def MPLsubplots(figure, nrows=1, ncols=1, sharex=False, sharey=False,
    293                  squeeze=True, subplot_kw=None, gridspec_kw=None):
    294         """
    295         Add a set of subplots to this figure.
    296        
    297         This is a copy of Figure.figure.subplots from matplotlib.figure.py.
    298         Included here because this appears only in later versions of
    299         matplotlib. When v2.2 is standard, this should be removed from GSAS-II.
    300        
    301         :param int nrows, ncols: default: 1
    302             Number of rows/cols of the subplot grid.
    303         :param bool/str sharex, sharey: True/False or {'none', 'all', 'row', 'col'}, default: False
    304 
    305             Controls sharing of properties among x (`sharex`) or y (`sharey`)
    306             axes:
     292# def MPLsubplots(figure, nrows=1, ncols=1, sharex=False, sharey=False,
     293#                  squeeze=True, subplot_kw=None, gridspec_kw=None):
     294#         """
     295#         Add a set of subplots to this figure.
     296       
     297#         This is a copy of Figure.figure.subplots from matplotlib.figure.py.
     298#         Included here because this appears only in later versions of
     299#         matplotlib. When v2.2 is standard, this should be removed from GSAS-II.
     300       
     301#         :param int nrows, ncols: default: 1
     302#             Number of rows/cols of the subplot grid.
     303#         :param bool/str sharex, sharey: True/False or {'none', 'all', 'row', 'col'}, default: False
     304
     305#             Controls sharing of properties among x (`sharex`) or y (`sharey`)
     306#             axes:
    307307           
    308                 - True or 'all': x- or y-axis will be shared among all
    309                   subplots.
    310                 - False or 'none': each subplot x- or y-axis will be
    311                   independent.
    312                 - 'row': each subplot row will share an x- or y-axis.
    313                 - 'col': each subplot column will share an x- or y-axis.
     308#                 - True or 'all': x- or y-axis will be shared among all
     309#                   subplots.
     310#                 - False or 'none': each subplot x- or y-axis will be
     311#                   independent.
     312#                 - 'row': each subplot row will share an x- or y-axis.
     313#                 - 'col': each subplot column will share an x- or y-axis.
    314314               
    315             When subplots have a shared x-axis along a column, only the x tick
    316             labels of the bottom subplot are visible.  Similarly, when
    317             subplots have a shared y-axis along a row, only the y tick labels
    318             of the first column subplot are visible.
    319         :param bool squeeze: default: True
    320 
    321             - If True, extra dimensions are squeezed out from the returned
    322               axis object:
     315#             When subplots have a shared x-axis along a column, only the x tick
     316#             labels of the bottom subplot are visible.  Similarly, when
     317#             subplots have a shared y-axis along a row, only the y tick labels
     318#             of the first column subplot are visible.
     319#         :param bool squeeze: default: True
     320
     321#             - If True, extra dimensions are squeezed out from the returned
     322#               axis object:
    323323             
    324                 - if only one subplot is constructed (nrows=ncols=1), the
    325                   resulting single Axes object is returned as a scalar.
    326                 - for Nx1 or 1xN subplots, the returned object is a 1D numpy
    327                   object array of Axes objects are returned as numpy 1D
    328                   arrays.
    329                 - for NxM, subplots with N>1 and M>1 are returned as a 2D
    330                   arrays.
     324#                 - if only one subplot is constructed (nrows=ncols=1), the
     325#                   resulting single Axes object is returned as a scalar.
     326#                 - for Nx1 or 1xN subplots, the returned object is a 1D numpy
     327#                   object array of Axes objects are returned as numpy 1D
     328#                   arrays.
     329#                 - for NxM, subplots with N>1 and M>1 are returned as a 2D
     330#                   arrays.
    331331                 
    332             - If False, no squeezing at all is done: the returned Axes object
    333               is always a 2D array containing Axes instances, even if it ends
    334               up being 1x1.
     332#             - If False, no squeezing at all is done: the returned Axes object
     333#               is always a 2D array containing Axes instances, even if it ends
     334#               up being 1x1.
    335335             
    336         :param dict subplot_kw: default: {}
    337             Dict with keywords passed to the
    338             :meth:`matplotlib.figure.Figure.add_subplot` call used to create
    339             each subplots.
    340         :param dict gridspec_kw: default: {}
    341             Dict with keywords passed to the
    342             :class:`matplotlib.gridspec.GridSpec` constructor used to create
    343             the grid the subplots are placed on.
     336#         :param dict subplot_kw: default: {}
     337#             Dict with keywords passed to the
     338#             :meth:`matplotlib.figure.Figure.add_subplot` call used to create
     339#             each subplots.
     340#         :param dict gridspec_kw: default: {}
     341#             Dict with keywords passed to the
     342#             :class:`matplotlib.gridspec.GridSpec` constructor used to create
     343#             the grid the subplots are placed on.
    344344           
    345         :returns: A single Axes object or array of Axes objects with
    346             the added axes.  The dimensions of the resulting array can be
    347             controlled with the squeeze keyword, see above.
     345#         :returns: A single Axes object or array of Axes objects with
     346#             the added axes.  The dimensions of the resulting array can be
     347#             controlled with the squeeze keyword, see above.
    348348           
    349         See also: pyplot.subplots: pyplot API; docstring includes examples.
    350        
    351         """
    352 
    353         # for backwards compatibility
    354         if isinstance(sharex, bool):
    355             sharex = "all" if sharex else "none"
    356         if isinstance(sharey, bool):
    357             sharey = "all" if sharey else "none"
    358         share_values = ["all", "row", "col", "none"]
    359         if sharex not in share_values:
    360             # This check was added because it is very easy to type
    361             # `subplots(1, 2, 1)` when `subplot(1, 2, 1)` was intended.
    362             # In most cases, no error will ever occur, but mysterious behavior
    363             # will result because what was intended to be the subplot index is
    364             # instead treated as a bool for sharex.
    365             if isinstance(sharex, int):
    366                 warnings.warn(
    367                     "sharex argument to subplots() was an integer. "
    368                     "Did you intend to use subplot() (without 's')?")
    369 
    370             raise ValueError("sharex [%s] must be one of %s" %
    371                              (sharex, share_values))
    372         if sharey not in share_values:
    373             raise ValueError("sharey [%s] must be one of %s" %
    374                              (sharey, share_values))
    375         if subplot_kw is None:
    376             subplot_kw = {}
    377         if gridspec_kw is None:
    378             gridspec_kw = {}
    379 
    380         # if figure.get_constrained_layout():
    381         #     gs = GridSpec(nrows, ncols, figure=figure, **gridspec_kw)
    382         # else:
    383         #     # this should turn constrained_layout off if we don't want it
    384         #     gs = GridSpec(nrows, ncols, figure=None, **gridspec_kw)
    385         gs = mpl.gridspec.GridSpec(nrows, ncols, **gridspec_kw)
    386 
    387         # Create array to hold all axes.
    388         axarr = np.empty((nrows, ncols), dtype=object)
    389         for row in range(nrows):
    390             for col in range(ncols):
    391                 shared_with = {"none": None, "all": axarr[0, 0],
    392                                "row": axarr[row, 0], "col": axarr[0, col]}
    393                 subplot_kw["sharex"] = shared_with[sharex]
    394                 subplot_kw["sharey"] = shared_with[sharey]
    395                 axarr[row, col] = figure.add_subplot(gs[row, col], **subplot_kw)
    396 
    397         # turn off redundant tick labeling
    398         if sharex in ["col", "all"]:
    399             # turn off all but the bottom row
    400             for ax in axarr[:-1, :].flat:
    401                 ax.xaxis.set_tick_params(which='both',
    402                                          labelbottom=False, labeltop=False)
    403                 ax.xaxis.offsetText.set_visible(False)
    404         if sharey in ["row", "all"]:
    405             # turn off all but the first column
    406             for ax in axarr[:, 1:].flat:
    407                 ax.yaxis.set_tick_params(which='both',
    408                                          labelleft=False, labelright=False)
    409                 ax.yaxis.offsetText.set_visible(False)
    410 
    411         if squeeze:
    412             # Discarding unneeded dimensions that equal 1.  If we only have one
    413             # subplot, just return it instead of a 1-element array.
    414             return axarr.item() if axarr.size == 1 else axarr.squeeze()
    415         else:
    416             # Returned axis array will be always 2-d, even if nrows=ncols=1.
    417             return axarr
     349#         See also: pyplot.subplots: pyplot API; docstring includes examples.
     350       
     351#         """
     352
     353#         # for backwards compatibility
     354#         if isinstance(sharex, bool):
     355#             sharex = "all" if sharex else "none"
     356#         if isinstance(sharey, bool):
     357#             sharey = "all" if sharey else "none"
     358#         share_values = ["all", "row", "col", "none"]
     359#         if sharex not in share_values:
     360#             # This check was added because it is very easy to type
     361#             # `subplots(1, 2, 1)` when `subplot(1, 2, 1)` was intended.
     362#             # In most cases, no error will ever occur, but mysterious behavior
     363#             # will result because what was intended to be the subplot index is
     364#             # instead treated as a bool for sharex.
     365#             if isinstance(sharex, int):
     366#                 warnings.warn(
     367#                     "sharex argument to subplots() was an integer. "
     368#                     "Did you intend to use subplot() (without 's')?")
     369
     370#             raise ValueError("sharex [%s] must be one of %s" %
     371#                              (sharex, share_values))
     372#         if sharey not in share_values:
     373#             raise ValueError("sharey [%s] must be one of %s" %
     374#                              (sharey, share_values))
     375#         if subplot_kw is None:
     376#             subplot_kw = {}
     377#         if gridspec_kw is None:
     378#             gridspec_kw = {}
     379
     380#         # if figure.get_constrained_layout():
     381#         #     gs = GridSpec(nrows, ncols, figure=figure, **gridspec_kw)
     382#         # else:
     383#         #     # this should turn constrained_layout off if we don't want it
     384#         #     gs = GridSpec(nrows, ncols, figure=None, **gridspec_kw)
     385#         gs = mpl.gridspec.GridSpec(nrows, ncols, **gridspec_kw)
     386
     387#         # Create array to hold all axes.
     388#         axarr = np.empty((nrows, ncols), dtype=object)
     389#         for row in range(nrows):
     390#             for col in range(ncols):
     391#                 shared_with = {"none": None, "all": axarr[0, 0],
     392#                                "row": axarr[row, 0], "col": axarr[0, col]}
     393#                 subplot_kw["sharex"] = shared_with[sharex]
     394#                 subplot_kw["sharey"] = shared_with[sharey]
     395#                 axarr[row, col] = figure.add_subplot(gs[row, col], **subplot_kw)
     396
     397#         # turn off redundant tick labeling
     398#         if sharex in ["col", "all"]:
     399#             # turn off all but the bottom row
     400#             for ax in axarr[:-1, :].flat:
     401#                 ax.xaxis.set_tick_params(which='both',
     402#                                          labelbottom=False, labeltop=False)
     403#                 ax.xaxis.offsetText.set_visible(False)
     404#         if sharey in ["row", "all"]:
     405#             # turn off all but the first column
     406#             for ax in axarr[:, 1:].flat:
     407#                 ax.yaxis.set_tick_params(which='both',
     408#                                          labelleft=False, labelright=False)
     409#                 ax.yaxis.offsetText.set_visible(False)
     410
     411#         if squeeze:
     412#             # Discarding unneeded dimensions that equal 1.  If we only have one
     413#             # subplot, just return it instead of a 1-element array.
     414#             return axarr.item() if axarr.size == 1 else axarr.squeeze()
     415#         else:
     416#             # Returned axis array will be always 2-d, even if nrows=ncols=1.
     417#             return axarr
    418418
    419419class _tabPlotWin(wx.Panel):   
    … …  
    29992999        Plot.set_visible(False)         #hide old plot frame, will get replaced below
    30003000        GS_kw = {'height_ratios':[4, 1],}
    3001         try:
    3002             Plot,Plot1 = Page.figure.subplots(2,1,sharex=True,gridspec_kw=GS_kw)
    3003         except AttributeError: # figure.Figure.subplots added in MPL 2.2
    3004             Plot,Plot1 = MPLsubplots(Page.figure, 2, 1, sharex=True, gridspec_kw=GS_kw)
     3001        # try:
     3002        Plot,Plot1 = Page.figure.subplots(2,1,sharex=True,gridspec_kw=GS_kw)
     3003        # except AttributeError: # figure.Figure.subplots added in MPL 2.2
     3004        #     Plot,Plot1 = MPLsubplots(Page.figure, 2, 1, sharex=True, gridspec_kw=GS_kw)
    30053005        Plot1.set_ylabel(r'$\mathsf{\Delta(I)/\sigma(I)}$',fontsize=16)
    30063006        Plot1.set_xlabel(xLabel,fontsize=16)
    … …  
    76427642        Plot.set_visible(False)
    76437643        GS_kw = {'width_ratios':[1,2],}
    7644         try:
    7645             Plot1,Plot = Page.figure.subplots(1,2,gridspec_kw=GS_kw)
    7646         except AttributeError: # figure.Figure.subplots added in MPL 2.2
    7647             Plot1,Plot = MPLsubplots(Page.figure, 1,2, gridspec_kw=GS_kw)
     7644        # try:
     7645        Plot1,Plot = Page.figure.subplots(1,2,gridspec_kw=GS_kw)
     7646        # except AttributeError: # figure.Figure.subplots added in MPL 2.2
     7647        #     Plot1,Plot = MPLsubplots(Page.figure, 1,2, gridspec_kw=GS_kw)
    76487648        Plot1.set_title('Line scan at azm= %6.1f'%Data['linescan'][1])
    76497649        Plot1.set_xlabel(r'$\mathsf{2\Theta}$',fontsize=12)

  • 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()

  • trunk/exports/G2export_PDB.py

    r3245 r4415  
    141141                'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
    142142            seq = 0
     143            notProt = True
    143144            for atom in Atoms:
    144145                if atom[ct-3] in AA3letter and int(atom[ct-4]) != seq:
     146                    notProt = False
    145147                    if atom[ct-2] not in seqList:
    146148                        seqList[atom[ct-2]] = []
    … …  
    172174                    Biso = atom[cia+1]*8.*np.pi**2
    173175                xyz = np.inner(A,np.array(atom[cx:cx+3]))
    174                 if atom[ct-3] in AA3letter:
     176                if atom[ct-3] in AA3letter or notProt:
    175177                    self.Write(fmt.format('ATOM  ',iatom,atom[ct-1],atom[ct-3].strip(),    \
    176178                        atom[ct-2].strip(),atom[ct-4].rjust(4),xyz[0],xyz[1],xyz[2],atom[cx+3], \

  • trunk/imports/G2phase.py

    r4213 r4415  
    5959#            return False
    6060        for i,l in enumerate(fp):
    61             if l.startswith('ATOM'):
     61            if l.startswith('ATOM') or l.startswith('HETATM'):
    6262                fp.close()
    6363                return True
    … …  
    7979        Phase = {}
    8080        Title = os.path.basename(filename)
     81        RES = Title[:3]
    8182       
    8283        Compnd = ''
    … …  
    8889        cell = None
    8990        Dummy = True
     91        Anum = 0
    9092        while S:
    9193            self.errors = 'Error reading at line '+str(line)
    … …  
    134136                    cell = [20.0,20.0,20.0,90.,90.,90.]
    135137                    Volume = G2lat.calc_V(G2lat.cell2A(cell))
    136                     AA,AB = G2lat.cell2AB(cell)                   
     138                    AA,AB = G2lat.cell2AB(cell)
     139                    Anum = 1                   
    137140                XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
    138141                XYZ = np.inner(AB,XYZ)
    … …  
    143146                if Dummy and S[12:17].strip() == 'CA':
    144147                    Type = 'C'
    145                 Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
    146                     S[12:17].strip(),Type.strip().capitalize(),'',XYZ[0],XYZ[1],XYZ[2],
     148                Aname = S[12:17].strip()
     149                if Anum:
     150                    Aname += '%d'%Anum
     151                if S[17:20].upper() != 'UNL':
     152                    RES = S[17:20].upper()
     153                Atom = [S[22:27].strip(),RES,S[20:22],
     154                    Aname,Type.strip().capitalize(),'',XYZ[0],XYZ[1],XYZ[2],
    147155                    float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
    148156                if S[16] in [' ','A','B']:
    … …  
    150158                    Atom.append(ran.randint(0,sys.maxsize))
    151159                    Atoms.append(Atom)
     160                if Anum:
     161                    Anum += 1
    152162            elif 'ANISOU' in S[:6]:
    153163                Uij = S[30:72].split()
Note: See TracChangeset for help on using the changeset viewer.