"""high-level initialization routines. These routines provide high level interfaces to initialize potential terms and to set up minimization and dynamics calculations. """ from atomSel import AtomSel #default topology/parameter files parameters={} parametersInitialized={} topology={} topologyInitialized={} parameters['protein'] = "protein.par" parametersInitialized['protein']=[] topology['protein'] = "protein.top" topologyInitialized['protein']=[] parameters['nucleic'] = "nucleic.par" parametersInitialized['nucleic']=[] topology['nucleic'] = "nucleic.top" topologyInitialized['nucleic']=[] parameters['water'] = "tip3p.parameter" parametersInitialized['water']=[] topology['water'] = "tip3p.topology" topologyInitialized['water']=[] def initRandomSeed(seed=None): """set the initial random seed. If the seed argument is omitted, it is set to seconds from the epoch mod 10^7 """ if globals().has_key('initialRandomSeed_'): raise Exception("random seed already initialized") global initialRandomSeed_ import time if seed==None: seed = int(time.time() % 1e7) import xplor seed = xplor.p_comm.distribute(seed) pass initialRandomSeed_ = seed import simulationWorld simWorld = simulationWorld.SimulationWorld_world() simWorld.setRandomSeed( initialRandomSeed_ ) if simWorld.logLevel() != 'none': print "random seed initialized to ", initialRandomSeed_ pass return def initialRandomSeed(): """return the initial random seed value, or that from .random if initRandomSeed has not been called. """ if globals().has_key('initialRandomSeed_'): return initialRandomSeed_ else: import simulationWorld simWorld = simulationWorld.SimulationWorld_world() return simWorld.random.seed() pass from xplorSimulation import getXplorSimulation def initParams(files, reset=0, weak_omega=0, simulation=0): """file is a structure type or filename or a list of a combination of these two. valid structure types are: protein, nucleic, and water If the argument is a filename, the current directory is first searched for the file, and then the TOPPAR directory is searched. XPLOR parameters are documented . If reset is true, all previous parameter settings are discarded If weak_omega is true, improper force constants associated with the peptide bond are reduced by 1/2 to allow a small amount of flexibility. """ import os from os import environ as env xSim = getXplorSimulation(simulation) command = xSim.command simName = xSim.name() if reset: command("param reset end") for key in parametersInitialized.keys(): if simName in parametersInitialized[key]: parametersInitialized[key].remove(simName) pass pass pass if weak_omega: command('eval ($weak_omega=1)') def addParams(file): if parameters.has_key(file): if simName in parametersInitialized[file]: return parametersInitialized[file].append(simName) file = parameters[file] pass if os.access(file,os.F_OK): pass elif os.access(env['TOPPAR'] + '/' + file,os.F_OK): file = 'TOPPAR:' + file else: mess = """initParams: could not find name %s as a structure type or as file in the current dir, or in TOPPAR""" % file raise Exception(mess) command("param @%s end" %file) return if type(files)==type("string"): files=[files] pass for file in files: addParams(file) pass return def initTopology(files, reset=0, simulation=0): """file is a structure type or filename or a list of a combination of these two. valid structure types are: protein, nucleic, and water First the current directory is searched for the file, and then TOPPAR is searched. XPLOR topology is documented . If reset is true, all previous topology settings are discarded """ import os from os import environ as env xSim = getXplorSimulation(simulation) command = xSim.command simName = xSim.name() if reset: command("rtf reset end") for key in topologyInitialized.keys(): if simName in topologyInitialized[key]: topologyInitialized[key].remove(simName) pass pass def addTopology(file): if topology.has_key(file): if simName in topologyInitialized[file]: return topologyInitialized[file].append(simName) file = topology[file] pass if os.access(file,os.F_OK): pass elif os.access(env['TOPPAR'] + '/' + file,os.F_OK): file = 'TOPPAR:' + file else: mess = """initTopology: could not find name %s as a structure type or as file in the current dir, or in TOPPAR""" % file raise Exception(mess) command("rtf @%s end" %file) return if type(files)==type("string"): files=[files] pass for file in files: addTopology(file) pass return def initStruct(files=0, erase=1, simulation=0): """read XPLOR PSF files. the file argument is a filename or list of filenames to read. First the current directory is searched for the file, and then TOPPAR is searched. Alternatively, a psf entry (starting with PSF...) can be directly passed as the files argument. Any pre-existing structure information is erased unless the erase argument is cleared. """ xSim = getXplorSimulation(simulation) command = xSim.command if erase: command("struct reset end") if not files: return import re, os if type(files)==type("string"): if re.search("PSF\s*\n",files): command("struct %s end" % files) return else: files=[files] pass pass for file in files: try: os.stat(file) command("struct @%s end" %file) except OSError: from os import environ as env tfile = env['TOPPAR'] + '/' + file try: os.stat(tfile) command("struct @TOPPAR:%s end" %file) except OSError: raise Exception("initStruct: could not find file " + file +\ " in current dir, or in TOPPAR") pass pass return pdbLocation="ftp://ftp.rcsb.org/pub/pdb/data/structures/divided/pdb/" def downloadPDB(entry): """download a compressed PDB file from the RCSB database return the PDB record as a string. """ entry = entry.lower() if len(entry)!=4: raise Exception("entry must be a 4 character string:" + `entry`) middle = entry[1:3] url = pdbLocation + middle + '/pdb' + entry + '.ent.Z' from simulationTools import mktemp tmpFile = mktemp(prefix='downloadPDB_%s'%entry) zFilename=tmpFile + ".Z" try: import urllib open(zFilename,'w').write(urllib.urlopen(url).read()) except: raise Exception("error fetching ",url) #this is a bit ugly- but I haven't bothered to code LZW decompression import os os.system('uncompress -f %s' % zFilename) ret = open(tmpFile).read() os.unlink(tmpFile) return ret def loadPDB(file='', string='', entry='', simulation=None, verbose=-1): """load a PDB entry from file, string, or from the RCSB database. autogenerate the psf and initialize coordinates. """ sim = simulation import simulation if sim: oldCurrentSimulation = simulation.currentSimulation() simulation.makeCurrent(sim) pass logfile=open("/dev/null","w",0) from simulationWorld import world import sys if verbose<0 and world().logLevel()!="none": verbose=True if verbose==True: logfile=sys.stdout if file: print >>logfile,'loading pdb file: %s ' % file, ; logfile.flush() pdbEntry = open(file).read() elif string: pdbEntry = string elif entry: print >>logfile,'loading pdb entry: %s ' % entry, ; logfile.flush() pdbEntry = downloadPDB(entry) print >>logfile,' [downloaded]', ; logfile.flush() else: raise Exception("entry or file argument must be specified.") import psfGen psfGen.pdbToPSF(pdbEntry) print >>logfile,' [psf]', ; logfile.flush() import protocol protocol.initCoords(string=pdbEntry,verbose=verbose) print >>logfile,' [coords]' if sim: simulation.makeCurrent(oldCurrentSimulation) pass return def initCoords(files=[],string="", model=-1, verbose=-1, erase=False, useChainID=True, maxUnreadEntries=(20,0.2), selection=0): """ Initialize coordinates from one of more pdb file, or from a string containing a PDB entry. if model is specified, the specified PDB MODEL record will be read. if erase is set to True, then all atom positions are cleared before reading the pdb record. If useChainID is True (the default), then the the chain ID PDB field is used as the segment name (if set), If the chain ID field is blank it is ignored regardless. PDB ATOM records which do not exactly match the expected input (the atom name is not exactly the same) are attempted to be matched using the function matchInexactAtomEntry. If more than maxUnreadEntries[0] or maxUnreadEntries[1]*(num residues) unmatchable ATOM records are encountered, an exception is thrown. """ if not selection: selection = AtomSel("all") pass if type(selection)==type("string"): selection = AtomSel(selection) pass xSim = getXplorSimulation(selection.simulation()) command = xSim.command if verbose<0: verbose=0 import simulationWorld simWorld = simulationWorld.SimulationWorld_world() if simWorld.logLevel() != 'none': verbose=1 pass pass if erase: # FIX: this should be an option to PDBTool::read command("coor INIT SELE=(all) end") pass if type(files)==type("string"): files=[files] pass if not files and not string: raise Exception("initCoords: a file or a string must be specified.") from pdbTool import PDBTool pdb = PDBTool("",selection,useChainID) pdb.setVerbose(verbose) unknownEntries=[] for file in files: pdb.setFilename(file) unknownEntries += pdb.read(model) pass if string: pdb.setFilename('') pdb.setContents(string) unknownEntries += pdb.read(model) pass unreadEntries=0 for entry in unknownEntries: if not matchInexactAtomEntry(entry,pdb,verbose): unreadEntries += 1 pass pass from selectTools import numResidues numResids = numResidues( selection ) if (unreadEntries>maxUnreadEntries[0] and unreadEntries>maxUnreadEntries[1]*numResids): raise Exception("too many unreadable ATOM entries: %d" % unreadEntries) if verbose: if unreadEntries: print "initCoords: unable to read %d pdb ATOM entries" % \ unreadEntries, print " (nonpseudoatom)" pass unknown = AtomSel("not known") if len(unknown): print "initCoords: still %d unknown atomic coordinates" % \ len(unknown) pass pass return # residue names of pseudo atoms which may or may not be present # in a coordinate record. pseudoResNames=['ANI'] def matchInexactAtomEntry(entry,pdbTool,verbose=0): """ try to match a PDB record to an atom if the name is not exactly the same. The coordinates will only be overwritten if they are unknown. The current algorithm for matching is as follows: match if there's only a single undefined atom whose name starts with the same character. else match if there's only a single undefined atom whose name starts with and ends with the same character. else match if the first name of the name is a digit present in the correct atom name and the second character matches the first character of the correct name. else match if the two names have two matching characters. else match up O --> OT1 OXT --> OT2 If the match is not unique, it is not made. The residue number and segment id must match. Entries with nonnumeric residue numbers are ignored. For atoms with residue names corresponding to those of pseudo atoms (e.g. ANI), this function prints a warning on missing atoms, but returns a True value. """ serial = entry[ 6 :12] name = entry[ 12:16] altLoc = entry[ 16:17] resName = entry[ 17:20] chainID = entry[ 21:22] try: resid = entry[ 22:26] except ValueError: if verbose: print 'nonnumeric residue number:' ,entry return 1 resSeq = entry[ 22:26] iCode = entry[ 26:27] x = float(entry[ 30:38]) y = float(entry[ 38:46]) z = float(entry[ 46:54]) occupancy = entry[ 54:60] tempFactor = entry[ 60:66] segID = entry[ 72:76] element = entry[ 76:78] charge = entry[ 78:80] if pdbTool.useChainID(): if chainID!=' ': segID = chainID+' ' pass readAtoms = pdbTool.readAtoms() name = name.strip() from atomSel import AtomSel resAtoms = AtomSel('not known and (resid %s and segid "%s")' % (resSeq,segID) ) resAtoms = filter(lambda a: not readAtoms[a.index()], resAtoms) matchAtom=None from psfGen import renameResidues resName = renameResidues([resName])[0] if resAtoms: recResName = resAtoms[0].residueName() if recResName != resName: raise Exception("residue name mismatch %s != %s\n" % (resName, recResName)+ " entry:" + entry) pass if not matchAtom: hn = filter(lambda a:a.atomName()=="HN",resAtoms) if hn and name=='H': matchAtom = hn[0] pass pass if not matchAtom: # match atom name that begin with same character for atom in resAtoms: if atom.atomName().startswith(name[0]): if matchAtom: matchAtom = None break else: matchAtom = atom pass pass pass pass if not matchAtom: # match atom name that begin and end with same character atom = filter(lambda a: (a.atomName()[0] ==name[0] and a.atomName()[-1]==name[0-1] ) , resAtoms) if len(atom)==1: matchAtom = atom[0] pass pass if not matchAtom: for atom in resAtoms: if ((name[0].isdigit() and atom.atomName().find(name[0])>=0) and name[1] == atom.atomName()[0]): if matchAtom: matchAtom = None break else: matchAtom = atom pass pass pass pass if not matchAtom: for atom in resAtoms: count=0 for c in atom.atomName(): if name.find(c)>=0: count += 1 pass if count>=2: if matchAtom: matchAtom = None break else: matchAtom = atom pass pass pass pass #rules for strange nucleic base HN atom naming if not matchAtom and name=="HN'": h21 = filter(lambda a:a.atomName()=="H21",resAtoms) h41 = filter(lambda a:a.atomName()=="H41",resAtoms) h61 = filter(lambda a:a.atomName()=="H61",resAtoms) if h21 and not readAtoms[h21[0].index()]: matchAtom = h21[0] elif h41 and not readAtoms[h41[0].index()]: matchAtom = h41[0] elif h61 and not readAtoms[h61[0].index()]: matchAtom = h61[0] pass pass #special rules for the O-terminus if not matchAtom: ot1 = filter(lambda a:a.atomName()=="OT1",resAtoms) if name=='O' and ot1 and not readAtoms[ot1[0].index()]: matchAtom = ot1[0] pass pass if not matchAtom: ot2 = filter(lambda a:a.atomName()=="OT2",resAtoms) if name=='OXT' and ot2 and not readAtoms[ot2[0].index()]: matchAtom = ot2[0] pass pass if not altLoc in pdbTool.allowedAltLoc(): matchAtom=False if matchAtom: if verbose: print "matchInexactAtomEntry:", \ "matching entry %s %s %s to atom %s %s %s" % \ (segID,resSeq,name, matchAtom.segmentName(),matchAtom.residueNum(), matchAtom.atomName()) matchAtom.setPos( [x,y,z] ) try: occupancy = float(occupancy) except ValueError: occupancy = 0 pass try: tempFactor = float(tempFactor) except ValueError: tempFactor = 0 pass pdbTool.setAux1(matchAtom, occupancy) pdbTool.setAux2(matchAtom, tempFactor) pass else: if verbose: print "matchInexactAtomEntry:", \ "found no match for entry %s %s %s" % (segID,resSeq,name) pass # print "matchInexactAtomEntry:", \ # "found no match for entry %s %s %s %s" % (segID,resSeq,resName, # name) pass if resName in pseudoResNames: if verbose: print " pseudo atom" matchAtom=True pass return matchAtom != None def addDisulfideBond(sel1,sel2): """ deprecated. Calls psfGen.addDisulfideBond """ import psfGen return psfGen.addDisulfideBond(sel1,sel2) def initNBond(cutnb=4.5, rcon=1.0, nbxmod=3, selStr="known", tolerance=0.5, repel=0.9, onlyCA=0, simulation=0): """standard initialization of the non-bonded repel potential. The XPLOR nonbonded potential term is described note that cutnb should be greater than rmax + 2*tolerance, where rmax is the largest vdw radius. selStr specifies which atoms to use in nonbonded energy calculations. If onlyCA is True, this string gets set to "name CA" - but this option is deprecated. """ noSelStr="resname ANI" #exclude these nonbonded interactions if onlyCA: selStr="name CA" xSim = getXplorSimulation(simulation) xSim.command(""" constraints interaction (%s and (not (%s))) (%s and (not (%s))) weights * 1 vdw 1 end interaction (not (%s) and (not (%s))) (not (%s)) weights * 1 vdw 0 end end""" % (selStr,noSelStr,selStr,noSelStr,selStr,noSelStr,noSelStr) ) xSim.command(""" parameters nbonds atom nbxmod %d wmin = 0.01 ! warning off cutnb = %f ! nonbonded cutoff tolerance %f repel= %f ! scale factor for vdW radii = 1 ( L-J radii) rexp = 2 ! exponents in (r^irex - R0^irex)^rexp irex = 2 rcon=%f ! actually set the vdW weight end end """ % (nbxmod,cutnb,tolerance,repel,rcon) ) def initRamaDatabase(type='protein', simulation=0): """standard initialization of the . """ xSim = getXplorSimulation(simulation) outputState=xSim.disableOutput() xSim.command(""" eval ($krama=1.) rama nres=10000 end """) if type=='protein': xSim.command(""" rama @QUARTS:2D_quarts_new.tbl @QUARTS:3D_quarts_new.tbl @QUARTS:forces_torsion_prot_quarts_intra.tbl end @QUARTS:setup_quarts_torsions_intra_2D3D.tbl """) xSim.enableOutput(outputState) elif type=='nucleic': xSim.command(""" evaluate ($knuc=1.0) rama @QUARTS:nucleic_deltor_quarts2d.tbl @QUARTS:nucleic_deltor_quarts3d.tbl @QUARTS:nucleic_deltor_quarts4d.tbl @QUARTS:force_nucleic_quarts2d.tbl @QUARTS:force_nucleic_quarts3d.tbl @QUARTS:force_nucleic_quarts4d.tbl end @QUARTS:setup_nucleic_2d3d.tbl @QUARTS:setup_nucleic_4d.tbl """) xSim.enableOutput(outputState) else: xSim.enableOutput(outputState) raise Exception("initRamaDatabase: unknown database type: "+ type) return def initDihedrals(filenames=[], string="", scale=1, useDefaults=1, simulation=0): """initialize the XPLOR . parameters are: filenames - either a single filename, or a sequence of filenames of dihedral restraint assignment tables. string - assignment table as a plain string. scale - scale factor (defaults to 1). useDefaults - use the default sidechain restraints (default: TRUE) these force chi2 angles of PHE, TYR, ASP, and chi3 of GLU to obey IUPAC naming conventions. """ xSim = getXplorSimulation(simulation) outputState=xSim.disableOutput() xSim.command(""" restraints dihed reset scale %f nass = 10000 end""" % scale) if type(filenames)==type("string"): filenames = [filenames] for file in filenames: xSim.command("restraints dihed @%s end" % file) pass if string: xSim.command("restraints dihed %s end" % string) if useDefaults: xSim.command(""" ! ! phe_angles.tbl ! ! constrains the chi2 angles of phe, tyr, asp, and glu ! in a protein to -90..90. This so that the dihedral angles follow ! IUPAC naming convention Biochemistry 9, 3471 (1970). ! ! JJK 3/10/97 ! for $res in id (name ca and (resn phe or resn tyr)) loop ang restraints dihedral assign (byresidue id $res and name ca) (byresidue id $res and name cb) (byresidue id $res and name cg) (byresidue id $res and name cd1) 1.0 0.0 90.0 2 end end loop ang for $res in id (name ca and resn asp) loop ang restraints dihedral assign (byresidue id $res and name ca) (byresidue id $res and name cb) (byresidue id $res and name cg) (byresidue id $res and name od1) 1.0 0.0 90.0 2 end end loop ang for $res in id (name ca and resn glu) loop ang restraints dihedral assign (byresidue id $res and name cb) (byresidue id $res and name cg) (byresidue id $res and name cd) (byresidue id $res and name oe1) 1.0 0.0 90.0 2 end end loop ang """) pass xSim.enableOutput(outputState) return def initCollapse(sel ="all", Rtarget=-1, scale =1): """initialize the XPLOR . parameters are: sel - string or atom selection specifying atoms in include in the Rgyr calculation. If this is omitted, all atoms are included. Rtarget - target radius. If this is omitted, the target Rgyr is calculated as Rgyr = 2.2*numResidues^0.38 -1 [angstrom] scale - scale factor. If omitted, it defaults to 1. Note that the per-assignment scale is always 100, so the energy is actually scaled by 100*scale. """ if type(sel)==type("string"): sel = AtomSel(sel) xSim = getXplorSimulation( sel.simulation() ) from selectTools import numResidues numResidues = numResidues(sel) if Rtarget<0: Rtarget = (2.2 * numResidues**0.38 -1) pass xSim.command(""" collapse assign (%s) 100.0 %f scale %f end""" % (sel.string(), Rtarget, scale)) return def initCarb13(filenames=[], scale=0.5): """initialize the XPLOR . arguments are: filenames - either a single filename, or a sequence of filenames of chemical shift assignment tables. In addition to these files C13SHIFTS:rcoil_c13.tbl and C13SHIFTS:expected_edited.tbl are always added. """ xSim = getXplorSimulation() if type(filenames)==type("string"): filenames = [filenames] restraints="" nres=0 for file in filenames: for line in open(file).readlines(): if line.strip().lower().startswith('assi'): nres += 1 restraints += line pass pass outputState=xSim.disableOutput() xSim.command(""" carbon phistep=180 psistep=180 nres=%d class all force %f potential harmonic set echo=off mess=off end @C13SHIFTS:rcoil_c13.tbl !rcoil shifts @C13SHIFTS:expected_edited_c13.tbl !13C shift database %s set echo=$prev_echo mess=$prev_messages end end""" % (nres,scale,restraints)) xSim.enableOutput(outputState) return def initHBDA(filename, scale=500, simulation=0): """initialize the XPLOR hbda potential term """ xSim = getXplorSimulation(simulation) lines=open(filename).readlines() nres = len(filter(lambda x: x.strip().lower().startswith('assi'),lines)) outputState=xSim.disableOutput() xSim.command(""" hbda nres %d class back @%s force %f end """ % (nres,filename,scale)) xSim.enableOutput(outputState) return def initMinimize(ivm, potList=0, printInterval=10, numSteps=500, maxCalls=20000, dEPred=0.001): """initialize an IVM object (from the module) for Powell minimization. In addition to the function arguments, the following IVM parameters are initialized: constrainLengths maxDeltaE eTolerance gTolerance """ ivm.resetReuse() # ivm.setConstrainLengths(0) ivm.setMaxDeltaE( 10000 ) ivm.setStepType("powell") ivm.setNumSteps( numSteps ) ivm.setMaxCalls( maxCalls ) ivm.setPrintInterval( printInterval ) ivm.setETolerance(1e-7) ivm.setGTolerance(1e-8) ivm.setDEpred(dEPred) if potList: ivm.setPotList( potList ) return def initDynamics(ivm, bathTemp=-1, finalTime=0.2, numSteps=0, stepsize=0.001, potList=0, printInterval=50, initVelocities=0, eTol_factor=0.001, eTol_minimum=0 ): """ initialize an IVM object (from the module) for PC6 (6th order predictor-corrector) dynamics. The default value of bathTemp is the ivm's current value. In addition to the function arguments, the following IVM parameters are initialized: constraintLengths maxDeltaE responseTime eTolerance adjustStepsize scaleVel resetCMInterval additionally, if initVelocities!=0, the initial velocities will be randomized. eTolerance is set by the following formula: eTolerance = eTol_factor * bathTemp + eTol_minimum """ from atomAction import randomizeVelocities ivm.resetReuse() eTol_temp=0 if bathTemp<0: bathTemp=ivm.bathTemp() ivm.setBathTemp(bathTemp) if initVelocities: randomizeVelocities( bathTemp ) eTol_temp=eTol_factor*ivm.bathTemp() # ivm.setConstrainLengths(0) ivm.setMaxDeltaE( 10000 ) ivm.setStepType("pc6") ivm.setResponseTime(5) ivm.setStepsize( stepsize ) #initial stepsize value ivm.setETolerance( eTol_temp + eTol_minimum) ivm.setPrintInterval( printInterval ) ivm.setAdjustStepsize(1) ivm.setScaleVel(1) ivm.setResetCMInterval( 10 ) ivm.setFinalTime(finalTime) ivm.setNumSteps(numSteps) if potList: ivm.setPotList( potList ) return def torsionTopology(ivm,fixedOmega=0,oTensors=[]): """configure the .IVM tolopogy for standard torsion angle setup: group rigid sidechains and break proline, ribose rings in the appropriate places. Disulfide bonds are also broken. This function sets all groupings and hinge types which are not already specified, so it must be called last in the setup of an IVM's topology. If fixedOmega is set, also fix protein omega backbone angles. oTensors is a list of .VarTensor objects used for alignment tensors """ import varTensorTools varTensorTools.topologySetup(ivm,oTensors) #check that all VarTensors are configured from atomSel import AtomSel aniOAtoms = AtomSel("resname ANI and name OO") oTensorOAtoms = map(lambda x:x.oAtom(), oTensors) for oAtom in aniOAtoms: if not oAtom.index() in map(lambda x:x.index(), oTensorOAtoms): print "torsionTopology: WARNING:", print "oTensor containing %s is not specified" % oAtom.string() pass pass import selectTools if fixedOmega: selectTools.IVM_groupRigidBackbone(ivm) selectTools.IVM_groupRigidSidechain(ivm) selectTools.IVM_breakProlines(ivm) selectTools.IVM_breakRiboses(ivm) selectTools.IVM_breakDisulfides(ivm) ivm.autoTorsion() return def cartesianTopology(ivm, sel="known",oTensors=[]): """configure the .IVM tolopogy for Cartesian dynamics/minimization - for the specified selection. This consists of breaking topological bonds, and specifying that all atoms are tree ``bases.'' This function should be called after any custom changes are made to the ivm's topology setup, but before torsionTopology(). oTensors is a list of .VarTensor objects used for alignment tensors """ import varTensorTools varTensorTools.topologySetup(ivm,oTensors) #check that all VarTensors are configured from atomSel import AtomSel import types if type(sel) == types.StringType: sel = AtomSel(sel) pass aniOAtoms = AtomSel("resname ANI and name OO") oTensorOAtoms = map(lambda x:x.oAtom(), oTensors) for oAtom in aniOAtoms: if not oAtom.index() in map(lambda x:x.index(), oTensorOAtoms): print "torsionTopology: WARNING:", print "oTensor containing %s is not specified" % oAtom.string() pass pass ivm.breakAllBondsIn("%s and not (resname ANI)" % sel.string()) ivm.setBaseAtoms(sel) return # decide how useful this code really is # it is used by eginput/dna_refi/ensemble.py def covalentMinimize(sel=0, numSteps=100, dEpred=1. ): """perform gradient optimization including only covalent terms (bonds, angles, impropers) """ from atomSel import AtomSel if isinstance(sel,str): sel = AtomSel(sel) if not sel: sel = AtomSel("known") import simulation oldCurrentSimulation = simulation.currentSimulation() sim = sel.simulation() xSim = getXplorSimulation(sim) simulation.makeCurrent( xSim ) xplorCommand = xSim.fastCommand outputState=xSim.disableOutput() prevPrintFile = xplorCommand("set print=off end","prev_print_file")[0] savedMass=[] indices=sel.indices() for atom in AtomSel('all'): savedMass.append(atom.mass()) if not atom.index() in indices: atom.setMass(0.) pass pass from xplorPot import XplorPot pots = map(lambda name:XplorPot(name),["bond","angl","impr"]) xplorCommand("""constraints interaction (attr mass>0) (attr mass=0) interaction (attr mass>0) (attr mass>0) end""") from ivm import IVM ivm=IVM() ivm.setVerbose(0) #don't print info messages ivm.fix("attr mass = 0") cartesianTopology(ivm,sel=sel) from atomAction import SetProperty AtomSel("all",sel.simulation()).apply(SetProperty("mass",100.)) initMinimize(ivm,potList=pots,numSteps=numSteps, dEPred=dEpred,printInterval=10) ivm.run() xplorCommand("set print=%s end"%prevPrintFile) xSim.enableOutput(outputState) def restoreContext(): for atom in AtomSel('all'): #restore masses atom.setMass(savedMass[ atom.index() ]) pass oldCurrentSimulation.sync() simulation.makeCurrent( oldCurrentSimulation ) pass restoreContext() return def fixupCovalentGeom(sel=0, rigidRegions=(), useVDW=0, useDynamics=1, dynRatio=5, maxIters=40, verbose=0, maxViols =0, bondTol =0.01, angleTol =2., torsionTol=2., extraTerms=[] ): """given the atoms in selection sel, perform, minimization and (optionally) dynamics to remove bond, angle, and improper violations - so that there total number is less than or equal to maxViols. rigidRegions is a sequence of selections which specify those regions to move as a rigid unit. If useVDW is set, the nonbonded term will be used 1/4 of the time. if useDynamics is set, dynamics will be used 1/dynRatio of the time. maxIters is the total maximum number of iterations. tolerances for bond, bond angle and improper torsion angles can also be specified using optional arguments. additional terms to satify may be specified with the extraTerms argument. Terms should be specified as a tuple of (pot,tol), where pot is a .Pot term and tol is the desired tolerance to be obtained. If this function is not successful in satisfying all covalent restraints, it will throw the exception protocol.CovalentViolation. This function changes the XPLOR constraints settings. If verbose>0, intermediate status is printed. Increase verbose for more verbosity. """ from atomSel import AtomSel if isinstance(sel,str): sel = AtomSel(sel) if not sel: sel = AtomSel("known") import simulation oldCurrentSimulation = simulation.currentSimulation() sim = sel.simulation() xSim = getXplorSimulation(sim) simulation.makeCurrent( xSim ) xplorCommand = xSim.fastCommand outputState=xSim.disableOutput() prevPrintFile = xplorCommand("set print=off end","prev_print_file")[0] savedMass=[] indices=sel.indices() for atom in AtomSel('all'): savedMass.append(atom.mass()) if not atom.index() in indices: atom.setMass(0.) pass pass def covalentViols(): ret = map(lambda (name,thresh): int(xplorCommand("print threshold %f %s end"%(thresh,name), "violations")[0]), (("bonds",bondTol), ("angles",angleTol), ("impropers",torsionTol))) for term in extraTerms: ret.append( term.violations() ) return ret from xplorPot import XplorPot potCombinations = map(lambda list:map(lambda name:XplorPot(name),list), [["bond"], ["bond","angl"], ["bond","impr"], ["bond","angl","impr"]]) if useVDW: initNBond(nbxmod=2,cutnb=6.5,tolerance=2.0,repel=0.8,rcon=0.01) potCombinations.append( potCombinations[-1]+[XplorPot('VDW')] ) pass for term in extraTerms: potCombinations.append( potCombinations[-1]+[term] ) if useVDW: potCombinations.append( potCombinations[-3]+[term] ) pass pass xplorCommand("""constraints interaction (attr mass>0) (attr mass=0) interaction (attr mass>0) (attr mass>0) end""") import random minNumViols=1e30 from ivm import IVM ivm=IVM() ivm.setVerbose(0) #don't print info messages ivm.fix("attr mass = 0") for aSel in rigidRegions: ivm.group(aSel) pass cartesianTopology(ivm,sel=sel) from atomAction import SetProperty actions = ["min"]*(dynRatio-1) if useDynamics: AtomSel("all",sel.simulation()).apply(SetProperty("mass",100.)) sel.apply(SetProperty("fric",10.)) actions.append("dyn") pass oldPots=None for iter in range(0,maxIters): viols=covalentViols() if verbose: print "iter", iter, " violations:",viols numViols=reduce(lambda x,y:x+y,viols) if numViols<=maxViols: break if numViols1: printInterval=1 pass if verbose>2: ivm.setVerbose( ivm.printResetCM | ivm.printVelFromCartCost | ivm.printTemperature | ivm.printEnergy | ivm.printCMVel | ivm.printStepDebug | ivm.printStepInfo | ivm.printNodeDef | ivm.printLoopDebug | ivm.printLoopInfo ) pass if verbose>3: ivm.setVerbose(ivm.verbose() | ivm.printNodePos) pass if verbose>4: ivm.setVerbose(ivm.verbose() | ivm.printNodeTheta) pass if action=="min": energy=abs(potList.calcEnergy()) initMinimize(ivm,numSteps=500,printInterval=printInterval, dEPred=energy/10., potList=potList) else: initDynamics(ivm,numSteps=500,stepsize=0.001, printInterval=printInterval, bathTemp=1000,initVelocities=True, potList=potList) pass ivm.run() pass xplorCommand("set print=%s end"%prevPrintFile) xSim.enableOutput(outputState) def restoreContext(): for atom in AtomSel('all'): #restore masses atom.setMass(savedMass[ atom.index() ]) pass oldCurrentSimulation.sync() simulation.makeCurrent( oldCurrentSimulation ) pass if numViols>maxViols: print "fixupCovalentGeom: Covalent geometry still violated at exit." print "using best structure:" xSim.setAtomPosArr( minCoords ) viols = covalentViols() mString= " (violations: bond: %d angle: %d improper: %d" % \ tuple(viols[:3]) for term in extraTerms: print "Violated restraints in term: ",term.instanceName() if hasattr(term,"restraints"): for r in term.restraints(): if r.violated(): print "%40s %10f"% (r.name(), r.diff()) pass pass pass mString += " %s: %d" % (term.instanceName(), term.violations()) pass mString += ")\n" print mString mess = "Covalent geometry still violated after fixupCovalentGeom\n" mess += mString if not useVDW: mess += "Try increasing maxIters or enabling the useVdw flag." else: mess += "Try increasing maxIters." pass restoreContext() raise CovalentViolation(mess) restoreContext() return class CovalentViolation(Exception): def __init__(s,mess): Exception.__init__(s,mess) pass def addUnknownAtoms(dyn_stepsize=0.02, dyn_numStepMul=1, verbose=0, simulation=0): """add in unknown atoms so that covalent and vdw terms are satisfied. This routine is slow, but it is rather robust. dyn_stepsize specifies the timestep size during the MD phase. Reduce this if you have convergence problems. dyn_numStepMul is a multiplier for the number of molecular dynamics steps taken. Increase this to get better convergence. if verbose=True, details of the minimization procedure are printed. """ xSim = getXplorSimulation(simulation) outputState=xSim.disableOutput() if not verbose: xSim.command("set print=off end") dyn_numStep = 500 * dyn_numStepMul dyn_ramp_numStep = 100 * dyn_numStepMul xSim.command(" evaluate ($timestep = %f)" % dyn_stepsize) xSim.command(" evaluate ($ramp_nstep = %f)" % dyn_ramp_numStep) xSim.command(" evaluate ($nstep = %f)" % dyn_numStep) xSim.command(r""" vector do (fbeta=10) (not known) {*Friction coefficient for MD heatbath.*} vector do (q=mass) (all) {* Save the real masses for later *} vector do (mass=100) (not known) {*Uniform heavy masses to speed*} vector do (fbeta=0) (known) {*Friction coefficient for MD heatbath.*} vector do (mass=0) (known) {*Uniform heavy masses to speed*} eval ($knoe=0.1) constraints fix (known) end vector do (vx = 0) (known) vector do (vy = 0) (known) vector do (vz = 0) (known) evaluate ($init_t = 1000 ) {* Initial simulated annealing temperature.*} vector do (vx = maxwell($init_t)) (attr mass>0) vector do (vy = maxwell($init_t)) (attr mass>0) vector do (vz = maxwell($init_t)) (attr mass>0) vector do (x=(random()-0.5)*20) (attr mass>0) vector do (y=(random()-0.5)*20) (attr mass>0) vector do (z=(random()-0.5)*20) (attr mass>0) !try bonds first flags exclude * include bond end constraints interaction (attr mass>0) (attr mass=0) interaction (attr mass>0) (attr mass>0) end !mini powell ! drop=1e5 ! nprint=1 ! tolg=1e-5 ! nstep=1000 !end flags exclude * include bond angle dihe cdihe impr vdw end !energy end !mini powell ! drop=1e5 ! nprint=1 ! tolg=1e-5 ! nstep=1000 !end evaluate ($kbon = 0.00005 ) {* Bonds. *} evaluate ($kang = 0.00005 ) {* Angles. *} !constraints ! interaction (all) (all) ! weights bond $kbon angl $kang impr $kimp vdw 0 elec 0 end !end !dynamics verlet ! nstep=5000 timestep=$timestep iasvel=current ! tcoupling=true tbath=$init_t nprint=50 iprfrq=0 !end while ($kbon < 0.01) loop stage1 evaluate ($kbon = min(0.25, $kbon * 1.25)) evaluate ($kang = $kbon) evaluate ($kimp = $kbon/10) noe scale * $knoe end !restraints dihed scale 0. end constraints interaction (attr mass>0) (attr mass=0) weights bond $kbon angl $kang impr $kimp vdw 5e-4 elec 0 end interaction (attr mass>0) (attr mass>0) weights bond $kbon angl $kang impr $kimp vdw 5e-4 elec 0 end end dynamics verlet nstep=$ramp_nstep timestep=$timestep iasvel=current tcoupling=true tbath=$init_t nprint=50 iprfrq=0 end end loop stage1 parameter {* Parameters for the repulsive energy term. *} nbonds repel=0.9 {* Initial value for repel - modified later. *} nbxmod=-3 {* Initial value for nbxmod - modified later. *} wmin=0.01 cutnb=4.5 ctonnb=2.99 ctofnb=3. tolerance=0.5 end end ! add vdw and slowly increase its weight parameter nbonds atom cutnb 100 tolerance 45 repel=1.2 rexp=2 irexp=2 rcon=1.0 nbxmod 4 end end flags exclude * include bond angle impr dihe cdihe vdw end constraints interaction (attr mass>0) (attr mass=0) weights bond $kbon angl $kang impr $kimp vdw 0.002 elec 0 end interaction (attr mass>0) (attr mass>0) weights bond $kbon angl $kang impr $kimp vdw 0.002 elec 0 end end dynamics verlet nstep=$nstep timestep=$timestep iasvel=current tcoupling=true tbath=$init_t nprint=50 iprfrq=0 end constraints interaction (attr mass>0) (attr mass=0) weights bond $kbon angl $kang impr $kimp vdw 0.005 elec 0 end interaction (attr mass>0) (attr mass>0) weights bond $kbon angl $kang impr $kimp vdw 0.005 elec 0 end end dynamics verlet nstep=$nstep timestep=$timestep iasvel=current tcoupling=true tbath=$init_t nprint=50 iprfrq=0 end constraints interaction (attr mass>0) (attr mass=0) weights bond $kbon angl $kang impr $kimp vdw 0.01 elec 0 end interaction (attr mass>0) (attr mass>0) weights bond $kbon angl $kang impr $kimp vdw 0.01 elec 0 end end dynamics verlet nstep=$nstep timestep=$timestep iasvel=current tcoupling=true tbath=$init_t nprint=50 iprfrq=0 end constraints interaction (attr mass>0) (attr mass=0) weights bond $kbon angl $kang impr $kimp vdw 0.5 elec 0 end interaction (attr mass>0) (attr mass>0) weights bond $kbon angl $kang impr $kimp vdw 0.5 elec 0 end end dynamics verlet nstep=$nstep timestep=$timestep iasvel=current tcoupling=true tbath=$init_t nprint=50 iprfrq=0 end !flags exclude * include bond angle end energy end mini powell drop=100 nprint=1 tolg=1e-5 nstep=10000 end constraints interaction (attr mass>0) (attr mass=0) weights bond $kbon angl $kang impr $kimp vdw 0.1 elec 0 end interaction (attr mass>0) (attr mass>0) weights bond $kbon angl $kang impr $kimp vdw 0.1 elec 0 end end mini powell drop=100 nprint=1 tolg=1e-5 nstep=10000 end vector do (mass=q) (all) {* Return the masses to sane values *} """) xSim.command("set print=OUTPUT end") xSim.enableOutput(outputState) return def genExtendedStructure(pdbFilename=0, sel=0, verbose=0, maxFixupIters=500 ): """ This assigns X, Y, and Z coordinates to each atom, and then calls .fixupCovalentGeom() to correct the covalent geometry, using maxFixupIters as the maxIters argument to that function. The Y and Z coordinates are random (but small enough (within a range of -0.5 to 0.5) to allow bonded atoms to form their bonds) and the X coordinate is the atom number divided by 10. This will result in an extended configuration along the X axis. If pdbFilename is nonnull, the file will be read if it exists, and that structure wil be used as a starting point (an attempt will be made to fix the the covalent geom). Whether it exists or not, it will be written to before returning. """ from atomSel import AtomSel if not sel: sel = AtomSel("not resname ANI") if type(sel)==type('string'): sel = AtomSel(sel) xSim = getXplorSimulation(sel.simulation()) from pdbTool import PDBTool import os if pdbFilename and os.path.exists(pdbFilename): PDBTool(pdbFilename,sel).read() else: from atomAction import SetProperty import random from vec3 import Vec3 for atom in sel: atom.setPos( Vec3(float(atom.index())/10, random.uniform(-0.5,0.5), random.uniform(-0.5,0.5)) ) pass pass if verbose: print "fixing covalent geometry..." pass fixupCovalentGeom(useVDW=1,maxIters=maxFixupIters,sel=sel,dynRatio=10, verbose=verbose) if pdbFilename: PDBTool(pdbFilename,sel).write() return