"""high-level initialization routines. These routines provide high level interfaces to initialize potential terms and to set up minimization and dynamics calculations. """ from xplor import command from atomSel import AtomSel #default topology/parameter files parameters={} parametersInitialized={} topology={} topologyInitialized={} parameters['protein'] = "protein.par" parametersInitialized['protein']=0 topology['protein'] = "protein.top" topologyInitialized['protein']=0 parameters['nucleic'] = "nucleic.par" parametersInitialized['nucleic']=0 topology['nucleic'] = "nucleic.top" topologyInitialized['nucleic']=0 parameters['water'] = "tip3p.parameter" parametersInitialized['water']=0 topology['water'] = "tip3p.topology" topologyInitialized['water']=0 def initParams(files, reset=0, weak_omega=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 if reset: command("param reset end") for key in parametersInitialized.keys(): parametersInitialized[key]=0 pass pass if weak_omega: command('eval ($weak_omega=1)') def addParams(file): if parameters.has_key(file): if parametersInitialized[file]: return parametersInitialized[file]=1 file = parameters[file] pass if file in os.listdir('.'): pass elif file in os.listdir(env['TOPPAR']): file = 'TOPPAR:' + file else: raise """initParams: could not find name %s as a structure type or as file in the current dir, or in TOPPAR""" % file 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): """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 if reset: command("rtf reset end") for key in topologyInitialized.keys(): topologyInitialized[key]=0 pass pass def addTopology(file): if topology.has_key(file): if topologyInitialized[file]: return topologyInitialized[file]=1 file = topology[file] pass if file in os.listdir('.'): pass elif file in os.listdir(env['TOPPAR']): file = 'TOPPAR:' + file else: raise """initTopology: could not find name %s as a structure type or as file in the current dir, or in TOPPAR""" % file 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): """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. """ 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 "initStruct: could not find file " + file +\ " in current dir, or in TOPPAR"; pass pass return def initCoords(files=[],string=""): """ Initialize coordinates from one of more pdb file, or from a string containing a PDB entry. """ if type(files)==type("string"): files=[files] pass if not files and not string: raise "initCoords: a file or a string must be specified." from pdbTool import PDBTool pdb = PDBTool() for file in files: pdb.setFilename(file) pdb.read() pass pdb.setContents(string) pdb.read() return def addDisulfideBond(sel1,sel2): """ add a disulfide bond between the two atom selections Should be called after PSF information is generated. """ from atomSel import AtomSel if isinstance(sel1,str): sel1 = AtomSel(sel1) if isinstance(sel2,str): sel2 = AtomSel(sel2) if len(AtomSel('segid "%s" and resid %d and name SG'% (sel1[0].segmentName(),sel1[0].residueNum())))==0: raise "No sulfur atom in "+sel1.string() if len(AtomSel('segid "%s" and resid %d and name SG'% (sel2[0].segmentName(),sel2[0].residueNum())))==0: raise "No sulfur atom in "+sel2.string() import xplor xplor.command("""patch DISU reference=1=( segid "%s" and resid %d ) reference=2=( segid "%s" and resid %d ) end""" % (sel1[0].segmentName(),sel1[0].residueNum(), sel2[0].segmentName(),sel2[0].residueNum())) return def initNBond(cutnb=4.5, rcon=1.0, nbxmod=3, tolerance=0.5, repel=0.9, onlyCA=0): """standard initialization of the non-bonded repel potential. The XPLOR nonbonded potential term is described If onlyCA is True, the nonbonded potential only acts between CA atoms. """ if onlyCA: command(""" constraints interaction (not name ca) (all) weights * 1 angl 0.4 impr 0.1 vdw 0 end interaction (name ca) (name ca) weights * 1 angl 0.4 impr 0.1 vdw 1.0 end end""") else: command(""" constraints interaction (all) (all) weights * 1 end end""") pass 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'): """standard initialization of the . """ (mess,echo) = command("set message off echo off end", ('prev_messages','prev_echo')) command(""" eval ($krama=1.) rama nres=10000 end """) if type=='protein': 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 set message %s echo %s end """ % (mess,echo) ) elif type=='nucleic': 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 set message %s echo %s end """ % (mess, echo)) else: raise "initRamaDatabase: unknown database type:", type return def initDihedrals(filenames=[], scale=1, useDefaults=1): """initialize the XPLOR . parameters are: filenames - either a single filename, or a sequence of filenames of dihedral restraint assignment tables. scale - scale factor (defaults to 1). useDefaults - use the default sidechain restraints (default: TRUE) """ command("set echo off mess off end") command(""" restraints dihed reset scale %f nass = 10000 end""" % scale) if type(filenames)==type("string"): filenames = [filenames] for file in filenames: command("restraints dihed @%s end" % file) pass if useDefaults: command(""" ! ! phe_angles.tbl ! ! constrains the chi2 angles of phe, tyr, asp, and glu ! in a protein to 0..180 degrees ! ! 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 90.0 30.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) 10.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) 10.0 0.0 90.0 2 end end loop ang """) pass command("set echo on mess on end") 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) from selectTools import numResidues numResidues = numResidues(sel) if Rtarget<0: Rtarget = (2.2 * numResidues**0.38 -1) pass command(""" collapse assign (%s) 100.0 %f scale %f end""" % (sel.string(), Rtarget, scale)) return def initHBDA(filename, scale=500): """initialize the XPLOR hbda potential term """ lines=open(filename).readlines() nres = len(filter(lambda x: x.strip().lower().startswith('assi'),lines)) command(""" set echo off mess off end hbda nres %d class back @%s force %f set echo $prev_echo mess $prev_messages end end """ % (nres,filename,scale)) 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 ): """ initialize an IVM object (from the module) for PC6 (6th order predictor-corrector) dynamics. 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 """ from atomAction import randomizeVelocities ivm.resetReuse() if bathTemp>=0: ivm.setBathTemp(bathTemp) if initVelocities: randomizeVelocities( bathTemp ) pass # ivm.setConstrainLengths(0) ivm.setMaxDeltaE( 10000 ) ivm.setStepType("pc6") ivm.setResponseTime(5) ivm.setStepsize( stepsize ) #initial stepsize value ivm.setETolerance( ivm.bathTemp()/1000 ) 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): """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. """ 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"): """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(). """ ivm.breakAllBondsIn(sel) ivm.setBaseAtoms(sel) return def covalentMinimize(numSteps=100, dEpred=1. ): """perform gradient optimization including only covalent terms (bonds, angles, impropers) """ import xplor xplor.command(""" flags exclude * include bonds angles impr end mini powell nstep=%d step=%f end """ % (numSteps,dEpred)) def fixupCovalentGeom(sel=0, useVDW=0, useDynamics=1, dynRatio=5, maxIters=40, verbose=0, ): """given the atoms in selection sel, perform, minimization and (optionally) dynamics to remove all bond, angle, and improper violations. 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. This function changes the XPLOR constraints settings. """ from atomSel import AtomSel import xplor (old_mess,old_echo) = xplor.command("set mess=off print=off echo=off end", ('prev_messages','prev_echo')) if isinstance(sel,str): sel = AtomSel(sel) if not sel: sel = AtomSel("known") xplor.command("""constraints fix (not (%s)) interaction (%s) (not (%s)) interaction (%s) (%s) end""" % (sel.string(),sel.string(),sel.string(), sel.string(),sel.string())) xplor.command("flags exclude * include scripting end") def covalentViols(): ret = map(lambda (name,thresh): int(xplor.command("print threshold %f %s end"%(thresh,name), "violations")[0]), (("bonds",0.01), ("angles",5.0), ("impropers",5.0))) return ret potCombinations = [("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(("bond","angl","impr","vdw")) pass import random from xplorPot import XplorPot potCombinations = map(lambda list:map(lambda name:XplorPot(name),list), potCombinations) from atomAction import SetProperty actions = ["min"]*(dynRatio-1) if useDynamics: sel.apply(SetProperty("mass",100.)) sel.apply(SetProperty("fric",10.)) actions.append("dyn") pass minNumViols=1e30 for iter in range(0,maxIters): viols=covalentViols() if verbose: print iter, viols numViols=reduce(lambda x,y:x+y,viols) if numViols==0: break if numViols0: print "fixupCovalentGeom: Covalent geometry still violated at exit." print "using best structure:" xplor.simulation.setAtomPosArr( minCoords ) print " (%d bond viols, %d angle viols, %d improper viols)\n" % \ tuple(covalentViols()) mess = "Covalent geometry still violated after fixupCovalentGeom\n" mess += ("(%d bond viols, %d angle viols, %d improper viols)\n" % tuple(viols)) if not useVDW: mess += "Try increasing maxIters or enabling the useVdw flag." else: mess += "Try increasing maxIters." pass raise mess pass return def addUnknownAtoms(dyn_stepsize=0.02, dyn_numStepMul=1, verbose=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. """ import xplor (old_mess,old_echo) = xplor.command("set mess=off echo=off end", ('prev_messages','prev_echo')) if not verbose: xplor.command("set print=off end") dyn_numStep = 500 * dyn_numStepMul dyn_ramp_numStep = 100 * dyn_numStepMul xplor.command(" evaluate ($timestep = %f)" % dyn_stepsize) xplor.command(" evaluate ($ramp_nstep = %f)" % dyn_ramp_numStep) xplor.command(" evaluate ($nstep = %f)" % dyn_numStep) xplor.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 *} """) xplor.command("set mess=%s print=OUTPUT echo=%s end"%(old_mess,old_echo)) return def genExtendedStructure(pdbFilename=0, sel=0, verbose=0, ): """ This assigns X, Y, and Z coordinates to each atom, and then calls .fixupCovalentGeom() to correct the covalent geometry. 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("all") if type(sel)==type('string'): sel = AtomSel(sel) 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=500,sel=sel,dynRatio=10, verbose=verbose) if pdbFilename: PDBTool(pdbFilename,sel).write() return