# mDyn - molecular dynamics program Purpose of program 1) relaxation of a molecular 3D structure via optimization of molecular conformational energy via local optimization methods 2) pseudo-global relaxation of conformational energy via simulation of thermal fluctuations of atoms via method of molecular dynamics in an explicit or implicit solvent using user defined positional or atom-atom distance restraints. 3) simulation of equilibrium molecular dynamics of protein 4) BLIND DOCKING OF FLEXIBLE LIGAND : a) calculation of low-resolution binding sites b) docking for user defined initial site of Protein c) automatic blind docking of Ligand on Protein via hierarchical method including refinement of ligand position/orientation and conformation via molecular dynamic simulated annealing coupled with force field deformation of ligand-protein interactions # USE of program 1)DOS installation make dir myMdynQ09/ put program mMdynQ09.exe to this directory subDirectories myMdynQ09/dat have to contain all file.dat - Lib data for program myMdynQ09/src fortran code files MAKE enviromenr variable > set MDYN09HOME="fullPath to myMdyn09/dat" myMdynQ09/doc directory contains Mannual and *.txt files explaining how to run program and undestand results of calculations RUN program from separate /jobX directory 2) Linux installation a) copy mDynQ09.tgz to myMdyn09/ directory b) extract archive > tar -zxvf mDynQ09.tgz c) set env variable $MDYN09HOMED > . ./setMDynQ09.sh 3) Compilation cd to source dir $MDYN09HOME/src do commands: g77 -c -O3 *.f g77 *.o -O3 -o mDynQ09 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 4) RUN program by the command ../$MDYN09HOME/mDynQ09 -i inProtcol -c inPDB [-cL 1etr.L.h.0.pdb] [-mdR mdRestXYZVin] [-mv moveRes] [-r1 inRestrainA1 ] [-r2 inRestrainA2] [-rB rigBodyFile] [-sa saProtocol] [-mn molName] [-mdX mdFinalPDB] [-tL ligMolTopoDat] [-bsX bsXYZinPDBfile] -o runOutFile [-er errorFile] in parenthesis [ ] are uxilarry files. The auxilary files will be used by program if the main command file defines the respective task. Command line DESCRIPTION: -i inProtcol : file MdynPar.inp defines protocol for the mDyn particular Run -c inPDB : file of the initial molecular structure as molec.pdb file in the PDB format -mdR mdRestXYZVin : XYZ+Velocity file to REstart MD from the last snapshot file XYZV , see exaple t5 1arb.mdXYZVfin0001.pdb it is USED with $mdRestart keyword in command file inProtcol NOTE! the initial XYZ will be taken from mdRestXYZVin file ! the PDB file inPDB is not USED with the key -mdR -r1 inRestrainA1 : file defines of positional restraints for atoms of the molecule -r2 inRestrainA2 : file defines atom-atom distance restraints -rB rigBodyFile : file defines rigid body segments of the main chain of protein -mv moveRes : file defines List of moving Residues -sa saProtocol : file defines simulated annealing protocol -mn molName : character set defining molecula name. molName. will be attached to RESULT files -cL 1etr.L.h.0.pdb - Ligand allAtom structure -tL ligMolTopoDat : file of LigandMolecTopology created by program mTopoHQ -bsX bsXYZinPDBfile: file of low-resolution binding site positions to be refined -o runOutFile : run output file -mdX mdFinalPDB : final PDB file of the Energy/MD optimization Current status of program run is printed on the standart output device (consol) or can be redirected to user defined file or can be defined in the argument line: -er errorFile : error message file : they are dublicated in the runOutFile # if file name definition in the argument line is missing for a file than the default name is used for this file NOTE! if the command line does not include a key -X , while the command file defines task which needs an additional data file coupled with -X keyword, than program try to find default (standart) name data file in the current directory. Default names: # inProtcol = ./MdynPar.inp inPDB = ./molec.pdb mdRestXYZVfile = ./mdXYZVin.pdb moveRes = ./moveRes.inp inRestrainA1 = ./restrAt1.inp' inRestrainA2 = ./restrAt2.inp' rigBodyFile = ./rigBody.inp saProtocol = ./SAprotocol.inp molName = space runOutFile = ./mDynSB.out errorFile = ./mDynSB.err mdFinalPDB = ./molMdFin.pdb # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # INPUT file description # inProtcol = ./MdynPar.inp The nain command file consist of lines with command keyword. Keyword start with $ sign in the first position of line One Keyword in line #example of MdynPar.inp file and keyword description # MdynPar.inp $OUTfull ! full extended output of program run ! by default OUTshort option is ON #Initial PDB data quality $Hread ! read INPUT pdb file with Hydrogens # DEfinition of OPtimized segments of protein: $fullProtMD ! full molecule is flexible $MovingRes ! defines List of opimized segments #FORCE FIELD MODIFICATIONS: # $shake=2 ! all valence bobds are fixed by shake method $shake=1 ! all H-X bonds are fixed ! default shake=0, all bonds are flexible, RECOMENDED # $zeroRot ! exclude translation and rotation of the molecule as rigid body $hBond128 = 2.0 ! scaling coeff for H-bonds ! default=1.0 is standart force field $harmAt1PosRst=0.25 !invoke restraintsA1 type = positional harmonic restraints for atom position with harmConst (kcal/A^2). program need a special file -r1 restrA1File which defines restrained segmants of protein see additional description $distRestrA2 !invoke restraintsA2 type atom-atom distances for user defined pairs of atoms in the file -r2 restrA2File (see additional description) $rigBody !invoke optimization with frozen internal structure of protein main chain for user defined segments of sequence need file -rB rigidBodySegment (see additional description) $compactForce = 0.5 ! invoke additional compactization forces ! to accelerate protein folding # $aSoftCore = 0.5 !invoke SOFTNES for the van der waals atom-atom potential ! at the small (contact) atom-atom distances ! Use of the softCore VDW potential helps to optimize ! BAD molecular structures with many spartial atom-atom clashes ! values range 0 - 1 from very Soft to standart VDW #SOLVATION MODEL $SolvMod = GShell # # # OPIMIZATION PROTOCOL: $engCalc ! do energy calculation $engOptim ! do energy optimization by local Optimizer $nOptStep=1 !max N optim steps # #PROTOCOL for Molecular Dynamics: $doMDyn ! do MolDynamics $MDSA !do MolecularDynamis SimAnnealing needs SAprotocolFile -sa saProtocol File, see additional description # #PROTOCOL of MD equilibration: # $initMDTemp=50.00 !initial Temperature to start MolDyn $bathMDTemp=50.00 !thermostat temperature of thermostat i.e. target temperature $runMDnstep=2000 !number of time-steps for MD simulation $mdTimeStep=0.002 # $NTV=1 ! MD ensemble definition # # # MD Trajectory writing: $nwtra=500 $WRpdb ! write snarshort structures in the PDB format ! default WRpdbq OPTion is ON : extended PDB format ! PDB + Qatom # END # NOTE that parameter file formatted, i.e. $ sign should be the firs character of the line ---------------------------------------------------------------------------- KEYWORD DESCRIPTION: #OUTPUT DETAILES: $OUTfull ! full extended output of program run ! by default OUTshort option is ON # # INPUT PDB FILE DETAILES: $Hread ! defines that all Hydrogens will be read from input molecule structure -c inPDB file otherwise the ALL HYDrogens will be restored by the program, i.e. all H atoms will be deleted and added according to molecular topology for RESidues. Using Library in the ./dat/h_add.dat NOTE! it is recommended start to works with a new protein without option $Hread even if the PDB file has all hydrogen atoms, because the hydrogen atom names for protein side chains have multiple definition in the PDB data base. It is better if mDyn program will add all hydrogens to the heavy atoms. #DEFINITION OF OPTIMIZED RESIDUES: $fullProtMD !defines FULL (i.e. ALL atoms) of the USER molecule will be free to move in energy relaxation or molDyn $MovingRes ! logical keyWord defines that only a defined set of RESidue are free this keyWord is coupled with file -mv moveRes in the argument line to start the program default name for moveRes file is ./moveRes.inp #EXAMPLE of ./moveRes.inp #1arb aaaaaaIIIIiiii # MOVRES 1 10 !line defines first and last resudue of moving segments integers devided by space MOVRES 45 76 MOVRES 115 260 end !end or END should be last line if the file ************ #FORCE FIELD DEFINITION: $hBond128 = 2.0 ! scaling coeff for H-bonds $aSoftCore = 0.5 !invoke van der waals atom-atom potential with modified ! SoftCore at the small (contact) atom-atom distances ! SoftCore modification is used for energyOPtimization and MD equilibration stages. ! Use of the softCore VDW potential helps to optimize ! BAD structures with many starical atom-atom clashes ! values range 0 - 1 from very Soft to standart VDW $harmAt1PosRst=0.25 ! digital keyWord define RESidue segments with 1 atom position harmonic restrants. 0.25 = harmonic restrain Constant K restrEnergy = 0.5*K(r - r0)**2, the reference position r0 = initialXYZinput.pdb - positions from the initial INPut PDB file which defines INItial structure of molecule this keyWord is coupled with file -r1 inRestrainA1 of the argument line to start the program mdyn default name for inRestrain file is ./restrAt1.inp #EXAMPLE of inRestrainA1 file: #harmonically restrained RESidue segments #xxxxxIIIIiiiiaaAAA #(6x,2i4,a40) RESTA1 1 63 PBB ! line starts from keyWord RESTAT numbers=first/last residue of segment ! PBB (only protein backbone atoms are restrained, i.e. side chains are free) RESTA1 78 120 ALL ! ALL (all atoms are restrained) ! integers and words are devided by space end # --------------------------------------------------- $distRestrA2 ! defines optimization/MD with atom-atom dist RestrainA2 ! needs file [-r2 inRestrainA2] in command line -r2 inRestrainA2 : default name : restrAt2.inp # EXAMPLE of inRestrainA2 file: #harmonically restrained Atom-Atom distances #xxxxxx #keyword atom1 atom2 distA HarmConst(kcal/mol*A^2) RESTA2 ND2 ASN 222 : OG1 THR 219 = 7.0 1.5 RESTA2 O GLY 170 : OG1 THR 219 = 8.0 2.5 RESTA2 OH TYR 109 : OG1 THR 111 = 7.5 3.0 END #---------------------------------------------------- $rigBody !defines optimizatiom/MD considering some segments of the main chain ! as a rigid body. ! The List of rigid segments of the main chain is user defined. ! Each segment will keep rigid internal structure of the protein main chain, ! has rotatational and translational degrees of freedom. ! The side chains of the rigid segments are flexible. #Needs file rigidBody.inp #EXAMPLE of rigidBody.inp file: # RIGB01 11 16 !line defines first and last resudue of moving segments integers devided by space RIGB02 47 59 RIGB03 77 99 end !end or END should be last line if the file # - - - - - - - - - - - - - - - - - - - - - - - - - $compactForce = 0.25 ! define additional compactization forces for protein atoms ! Recomended forceParameter = 0.1 - 1.0 # -------------------------------------------------- $shake=2 ! invoke shake subroutine to keep bonds fixed. =1 -bonds with Hydrogen, =2 all bonds ---------------------------------------------------- #Definition of the SOLVation model: there are 4 variants of Implicit models 1 variant of Explicit model #: $SolvMod = GShell ! implicit Gaussian Shell solvation model $SolvMod = GShell + WBrg ! implicit Gaussian Shell solvation model + WaterBridges between polar atoms ! WaterBridges descride solvent mediated interactions trough stong bound water ! molecules via implicit model of water bridges $SolvMod = GBorn ! implicit Generalized Born model + SAS HydroPhobic solvation $SolvMod = GBorn + WBrg ! implicit Generalized Born model + SAS HydroPhobic solvation + WaterBridges $Solv = ExWshell 4.5 [A] ! explicit water shell of 4.5 Angst around protein; ! recomended thikness 3.0 - 6.0 A --------------------------------------------------- $mdRestart ! restart molDynamics from a snapshot [molName.]mdXYZVfin000N.pdb the file [molName.]mdXYZVfin000N.pdb should be copied to the file mdyn Restart file mdXYZVin.pdb $doMDyn ! do molecular dynamics $MDSA ! do Molecular Dynamical Simulated Annealing ! coupled with file -sa SAprotocol which define protocol of the simulated annealing #EXAMPLE of Aprotocol.inp file #SA protocol #nSAstep 2 #(f10.1,1x,f8.1,1x,3(f6.1,1x) # nMDstep tempTg SCvdW wfHb128BB wfhB128BS SAPROT 100000 500.0 0.8 1.0 1.0 !line starts from keyword SAPROT SAPROT 100000 100.0 1.0 1.0 1.0 END # nMDstep - number of md timeStep tempTg - target temperature in K, this temperature will be reach during ntimeMX steps SCvdW - parameter 0 - 1 to define softness of the van der waals potential. Soft potential modifies Potential Energy Surface and decrease barriers of conformational transitions wfHb128BB, wfhB128BS - (1 - 0) scaling factors for BackBone-BackBone and BackBone-SideChain Hydrogen Bond energy #-------------------------------------------------- # # OPIMIZATION PROTOCOL: $engCalc ! do energy calculation $engOptim ! do energy optimization by local Optimizer $nOptStep=1 !max N optim steps # #PROTOCOL for Molecular Dynamics: $doMDyn ! do MolDynamics $MDSA !do MolecularDynamis SimAnnealing needs SAprotocolFile -sa saProtocol File, #MD EQUILIBRATION: $initMDTemp=50.00 !defines initial temperature to start MD ! recommended low temperature < 50K ! temperature can be steadelly increased to the 300K and higher ! USING $MDSA option $bathMDTemp=50.00 ! bath temperature in the MD equilibration run $runMDnstep=2000 ! number of MD time steps in the equilibration run $mdTimeStep=0.002 ! value of the MD time step in ps, ! recomended 0.001 - 0.002 $NTV=1 ! ansemble NTV=0/1 ! =1 md run with constant T #MD TRAJECTORY WRITING $nwtra=500 ! structure XYZ (snapshot) will be written !as a series of molMdResXXXX.pdb files $WRpdb ! write snapshort structures in the PDB format ! default is WRpdbq OPTion is ON : extended PDB format ! PDB + Qatom column #* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # -c inPDB file - standart pdb file #EXAMPLE of inPDB file: **************************************************************************************** NOTE! it is recommended to start to work with a new protein without option $Hread even if the PDB file has all hydrogen atoms, because the hydrogen atom names for protein side chains have multiple definition in the PDB data. It is better if mDyn program will add all hydrogens to the heavy atoms. ******************************************************************************************* REMARK: PDB: ATOM 1 N GLY A 1 11.726 -10.369 10.598 ATOM 2 H1 GLY A 1 11.921 -11.015 9.807 ATOM 3 H2 GLY A 1 12.518 -10.395 11.271 ATOM 4 H3 GLY A 1 10.852 -10.663 11.079 ATOM 5 CA GLY A 1 11.567 -9.015 10.090 ATOM 6 HA2 GLY A 1 10.772 -8.977 9.420 ATOM 7 HA3 GLY A 1 12.439 -8.710 9.612 ATOM 8 C GLY A 1 11.280 -8.099 11.303 ATOM 9 O GLY A 1 11.256 -8.584 12.493 ATOM 10 N VAL A 2 11.060 -6.876 11.020 ATOM 11 H VAL A 2 11.066 -6.574 10.025 etc. TER ! CHAIN TERmination ATOM 1302 N GLY A 94 10.957 -15.678 12.832 ATOM 1303 H GLY A 94 10.735 -14.663 12.877 ATOM 1303 H GLY A 94 10.735 -14.663 12.877 ATOM 1304 CA GLY A 94 10.193 -16.559 11.950 ATOM 1305 HA2 GLY A 94 9.428 -16.004 11.516 ATOM 1306 HA3 GLY A 94 9.784 -17.323 12.525 ATOM 1307 C GLY A 94 11.016 -17.184 10.843 ... etc. TER ! CHAIN TERmination END ! file END * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * # # PDB mDyn trajectory file description: # Program mDyn generate a series of snapshot files, e.g., 1arb.molMdRes0nnn.pdb (test/t4) the molMdResXXXX.pdb file (see example) contains all atomic coordinates and additional information in the REMARK: lines #### REMARK: Md result : MdTime(ps): 2.4940 REMARK: $nstep: 1247 REMARK: $nRecPDB: 5 REMARK: RMSD(x0): 0.43 <- RMSD all atom REMARK: badBond: n,erAv(A) : 0 0.000 <- number and error Average for bond length in Angstrem REMARK: badAng : n,erAv(grd): 8 9.42 <- number and error Average for bond angles in grad # ENERGY TERMS for the given structure REMARK: $ENERGY: :Kcal REMARK: eVbondDef: 100.89315 <-bond deformation energy REMARK: eVangDef : 441.63705 <-angle deformation energy REMARK: eImpDef : 35.68147 <-Improper torsion agle [planarity] energy REMARK: eTorsDef : 691.25769 <-torsion potentioal energy REMARK: engVDWR1 : -1031.16211 <- van der waals energy for cutoff R1=8 A REMARK: ehBHxY128: -608.70599 <- H-bondinds energy REMARK: engCOULR1: -816.25323 <- COULOMBIC for distances < cutoff R1 REMARK: engCOULR2: -4.47208 <- COULOMBIC for distances Rij, R1< rij < R2=14 A REMARK: restr1Eng: 2.37470 <- atomA1positional restrain eng REMARK: eRstA2Eng: 137.93436 <- atom-atom distances restrain eng REMARK: eCompcEng: 3.04940 <- artificial compactization energy REMARK: eGeoDef : 1269.46936 <- total ValenceGeometryDeformation eng REMARK: engCOUL : -820.72534 <- total COIUL eng REMARK: engSolvGS: -311.56546 <- GaussianShell solvation model eng REMARK: engSolHP : 0.00000 <- Hydrophocic cavity solvation energy REMARK: bornPolz : 0.00000 <- Generalized Born model polarization eng REMARK: engSolWatBrg: 0.00000 <- Water Bridges solvation energy REMARK: engSolTot: -311.56546 <-total solvation energy REMARK: engPOTENT: -1502.68958 <- total potential energy REMARK: kinEng : 517.13666 <- thermal kinetic energy REMARK: Temp(K) : 87.45163 <- effective temperature (K) of atomic motions REMARK: engTotal : -985.55292 <- total energy = Potential + Kinetic REMARK: REMARK: PDB: X Y Z Q atom ATOM 1 N GLY A 1 11.289 -9.838 10.638 0.29430 ATOM 2 H1 GLY A 1 11.186 -10.306 9.749 0.16420 ATOM 3 H2 GLY A 1 12.210 -9.985 11.026 0.16420 etc. TER END REMARK: GEOMETRY ANALYSYS: REMARK: BAD BOND LIST: REMARK: BAD BOND 1 C ARG 156 -N ILE 157 : 1.3990 err: 0.0640 REMARK: BAD ANGL LIST: # ANGLE between atoms at1 at2 at3 angl(grad) REMARK: BAD ANGLE 1 C GLY 1 -N VAL 2 -CA VAL 2 : 129.771 err: 8.071 REMARK: BAD ANGLE 2 N SER 5 -CA SER 5 -CB SER 5 : 102.698 err: -7.802 REMARK: BAD ANGLE 3 O SER 5 -C SER 5 -N CYX 6 : 115.473 err: -7.527 REMARK: BAD ANGLE 4 CB ILE 8 -CA ILE 8 -C ILE 8 : 122.456 err: 12.356 REMARK: BAD ANGLE 5 C ASP 9 -N VAL 10 -CA VAL 10 : 128.773 err: 7.073 etc. REMARK: BAD ANGLE 51 C ASP 259 -N GLY 260 -CA GLY 260 : 128.458 err: 6.758 REMARK: END