MolDyn description


I. Instalation and RUN

See "Instalation and RUN" topic.

II. Program flow and Basic algorithms of the program

1. Main program

 Main Program file  :       MDynSBmain.f

Start from the call of the input parameters

1.	call inputMDSApar    
reads the main Input file 
filenam = './MdynPar.inp'        ! in current job_dir

the file has the fixed name and located in the current job directory
the main input file  MdynPar.inp defines main parameters of the job (see chapter input file description)

2.	call initMolecTopSeq01  

reads  a defined  molecular PDB file, which can be defined in the  MdynPar.inp file
or  has the standard name   ./molec.pdb     and located in the current job directory ./ ;
defines  residue sequence

     3.  call initMolecTopSeq02       
calculates  12neighbour list (covalent bonds connecting atoms) using a predefined topology 
information about resdues stored in the $MDSBHOME/dat 

the pair12 list array: pair12List(*)  is the basic molecular topology information.
Based on the pair12List(*)  the all other lists are calculated, namely
Bonded triplets and  quartets to form list of  covalent angles,  torsion angles, improper torsion angles. 
The list of triplets and quartets are calculated via tree  algorithm

      Call       vbondListPDB2(atomXYZ,
     &           natom,atomNumb,atomName,resName,chName,resNumb,
     &           nres,resNameRes,chNameRes,
     &           atomNameEx,startAtInRes,
     &           nmoveatom,moveAtomList, 
     &           pair12List,startPairL12,nPairL12,np12MAX,
     &           pair13List,startPairL13,nPairL13,np13MAX,
     &           pair14List,startPairL14,nPairL14,np14MAX,
     &           bond12List,nbond12,
     &           trip123List,nTrip123,np123MAX,
     &           quar1234List,nQuar1234,np1234MAX,
     &           quarImp1234L,nImp1234,nImp1234MAX)

the call of the subroutine initMolecTopPDB  results in the complete definition of the molecular topology 
from the input molec.pdb  3D structure.

4.	call initFFieldParam

Initialization of the force field parameters for the bond, angle, torsion angle, improper angle deformations, 
 van der waals non bond interactions and  atomic point charges for the electrostatic interactions.  
 For bond, angle, torsion and improper angles a respective list of parameters are generated and stored in the arrays.

A list All force  field parameters are based on the amber94 force field parameter set [Cornell  1995]. 
 Molecular mechanical energy is based on the standard  equations for the force field of second generation 
 amber94 [Cornell  1995].  
Decoding of the atom names (residue names) to the forceField atom name is based on the look up table     
ffAtomTypeFile =  $MDSBHOME/dat/atmAAmberff.dat   

5.	Extraction of the data from Library file
All search of the proper names in the look up table of the MDynSB program are based on 
the hashing of  a records in the look up table, i.e. conversion of the table in numerically 
sequential order. If   several records of the look up table have the same hash number (degenerated case), 
they are placed  in a linkedLis for this hash number. 

Force field parameters are taken from the file:
ffParFile = $MDSBHOME/dat/bsparBATV.dat
code fragment to initialize force field parameters
c get ff-atom code from atomNames
	 call defFFatomName (ffAtomTypeFile,
     &              natom,atomNameEx,ResName,chName,
     &              ffAtomName,atomQ)
c define bondDef parameters for pair12List()
        call getBondDefPar(ffParFile,
     &              natom,atomNameEx,ResName,chName,ffAtomName,
     &              bond12List,nbond12,bond12ParL)
c c define valence angles def parameters
        call getVangDefPar(ffParFile,
     &              natom,atomNameEx,ResName,chName,ffAtomName,
     &              trip123List,nTrip123,ang123ParL)

c define Improper angle def parameters
        call getImpDefPar(ffParFile,
     &              natom,atomNameEx,ResName,chName,ffAtomName,
     &              quarImp1234L,nImp1234,impAng1234ParL)  

c define torsion parameters
        call getTorsPar(ffParFile,
     &        natom,atomNameEx,ResName,chName,ffAtomName,
     &        quar1234List,nQuar1234,quar1234ParL,quar1234nPar)
c assign atomMass and vdwParameters
        call  getVDWatMass(ffParFile,
     &              natom,atomNameEx,ResName,chName,ffAtomName,
     &              nVDWtype,atomVDWtype,atomVDW12ab,atomMass) 
c all FField Parameters are defined

6.	call initSolvatGSmod   

Defines  atomic parameters of the current structure for  the Gaussian Shell implicit solvation model [Lazaridis, 1999]. 
A parameters of the GS model are stored in the files:

7.	call initMDStart(tempT0)

Initialize MD calculation: 

Calculate the Initial nonBondPair lists
c generate three nonbonded atom pair Lists: van der Waals, Coulombic and solvation model.

        makeVdW = 1
        makeCL = 1
        makeSL = 1
        call initNonBondList(atomXYZ,makeVdW,makeCL,makeSL) 

Calculates the forces on atoms for initial atomic coordinates 
initial forces on atoms 
        fcall = 0
	call initAllForce(fcall,atomXYZ,makeVdW,makeCL,makeSL,
     &              eVbondDef,vbdefForce,
     &              eVangDef,vAngdefForce,
     &              eImpDef,impDefForce,
     &              eTorsDef,torsAngForce,
     &              engVDWR1,vdwForceR1,
     &              engCOULR1,coulForceR1,
     &              engCOULR2,coulForceR2,
     &              restr1Eng,restr1AtForce,
     &              molSolEn, atomSolEn,atomSolFr)

Calculates initial atomic velocities, which are distributed according to Maxwell law 

                  probability(vi) = ( ) exp(-mivi2/kT)

	call initVelocity(temp,natom,
     &       nmoveatom,moveAtomList,atomMass,atomVel0)

8.	Run MD 

The subroutine mdRun perform MD run for a given  number of time steps ntimeMX

	call mdRun(ntimeMX,ntime0,ntime,ntimeR1,ntimeR2,
     &             ntimeF1,ntimeF2,ntimeF3,deltat, 
     &             tempTg,tauTRF,atype,optra,wtra,nwtra,cltra)

9.	 Simulated Annealing optimization

	call simAnnealing(nSAstep,SAProtcol)
with user defined SAProtocol(nstep,T) consisted of nSAstep.

Each step of the SA is  MD run of  nstep with particular temperature T.

III. Details of the atomic force calculation

All atoms of the molecular system consists of two sets of fixed and moving atoms.
The force are calculated only for the moving atom set.

1. Covalent bond deformation

For covalent bond deformation we use the GROMOS functional form

rij = ri - rj
bn = rij .
This functional form is equivalent to the usual harmonic function for a small deformations but a computationally is more effective.
Force on atom i due to bond n

Total bond deformation force on atom i is the sum over all bonds n involving the atom i.
The calculation of the force fin is doing by

subroutine vbonddefenf(xyz1,xyz2,bondPar,edef,f1,f2) (see file vdefenforce.f)

2. Covalent angle deformation

The covalent angle deformation energy function has the form

This functional form is equivalent to the usual harmonic function for the angles for a small angle deformation but a computationally is more effective. The angle 2n ( at the j ) is between atoms i-j-k . The cosine of the angle 2n

The forces on atoms i,j,k due to the deformation of the angle 2n

respectively force on atom k

force on atom j is given from the conservation of the total force acting on three atoms

The covalent angle deformation energy and force are calculated in subroutine

subroutine vangldefenf(xyz1,xyz2,xyz3,angPar,
     &                         edef,f1,f2,f3)  

(see file vdefenforce.f)

3. Torsion angle energy and force

The total torsion energy is a sum over a set of torsion angles for the four atoms i-j-k-l with a rotation around bond j-k ,

where torsion energy for bond j-k can have several torsion barriers with different multiplicity.
Torsion angle N is defined as

The forces on atoms i,j,k,l due to the single term of eq.(8b) are

The torsion energy and force are calculated via

	subroutine torsanglenf(xyz1,xyz2,xyz3,xyz4,nTorsH,
     &                          torsPar,eTors,f1,f2,f3,f4)

c torsPar(4*nTorsH) = {pass,Vt/2/pass,cos(delta),nFi },
c  eTors = sum{ Ki*[1+cos(delti)cos(i*Ftors)] }; i=1,..,nTorsH

Torsion parameters are taken from the LibData = bsparBATV.dat
The extraction of the torsion parameters from LibData = bsparBATV.dat for all quartets is done by

       subroutine getTorsPar(ffParFile,
     &              natom,atomNameEx,ResName,chName,ffAtomName,
     &              quar1234L,nQuar1234,quar1234Par,quar1234nPar)        
c InPut:
c       ffParFile - ffParameters file 
c       natom,atomNameEx,ResName,chName : PDB info
c       ffAtomName(ia) - FFatomName to search table
c       the quar1234L(i),i=1,..,nQuar1234 : the QuartetList
c RESULT: quar1234Par(16*nQuar1234) - torsionFF parameters for list 
c         of  quartets 
c         pass,Vt/2,delta,nFi - (printed) for each torsHarmonics,
c         pass,Vt/2/pass,cos(delta),nFi - finally in array
c        4- torsionHarmanics is possible.
c        quar1234nPar(iQuart) - number of torsHarmonics for the torsAngl

4. Improper Torsion Angle (out of plane) deformation

The improper torsion angle deformation keeps the four atoms 1-2-3-4 (i-j-k-l ) in specified geometry. The first atom in the improper quartet is a planar or (tetrahedral) atom. For example atoms Ci-CAi-N(i+1)-Oi are kept planar. The out of plane potential

CA-N-C-CB are kept in the tetrahedral configuration (L-amino acid) or CA-C-N-CB (D-amino acid) if CA in the united atom (CH) presentation.
The out of plane angle is defined for j-i-k four atoms with i is the planar (tetrahedral)
angle between to planes (i-j-k) and (j-k-l) with rotation angle around j-k, other words the
torsion angle in the sequence i-j-k-l


The forces on atoms i,j,kl due to a single term Vn

finally from the third Newton law

The improper energy and forces for a given improper quartet of atoms are calculated by the subroutine

c improper torsion energy force
        subroutine imprtorsanglenf(xyz1,xyz2,xyz3,xyz4,impPar,
     &                         eImpt,f1,f2,f3,f4)

c ImptPar(2) = K1, ksi0

5. Covalent back-bond deformation calculation

All valence back-bond deformation are calculated in the file initAllForce.f 

subroutine initAllForce(fcall,atomXYZ,
     &              makeVdWs,makeCLs,makeSLs,
     &              eVbondDef,vbdefForce,
     &              eVangDef,vAngdefForce,
     &              eImpDef,impDefForce,
     &              eTorsDef,torsAngForce,
     &              engVDWR1,vdwForceR1,
     &              engCOULR1,coulForceR1,
     &              engCOULR2,coulForceR2,
     &              restr1Eng,restr1AtForce,
     &              molSolEn, atomSolEn, atomSolFr)
        include 'xyzPDBsize.h'
        include 'xyzPDBinfo.h'
        include 'pair1234array.h'
        include 'nbondPairVCS.h'
        include 'vdw12Par.h'
        include 'restrainInfo.h'
        include 'loopInfo.h'
        include 'movingAtom.h'
        include 'solvGSarray.h'
        include 'optionPar.h'
 . . . . . . . . . . . . . . . . . . . . . 

c all GeoDef forces are calculated at each step

 	call allAtVBondEForce(atomXYZ,
     &           natom,bond12List,nbond12,bond12ParL,
     &           eVbondDef,vbdefForce )
	call allAtVangEForce(atomXYZ,
     &           natom,trip123List,nTrip123,ang123ParL,
     &           eVangDef,vAngdefForce )
	call allAtImpTEForce(atomXYZ,
     &           natom,quarImp1234L,nImp1234,impAng1234ParL,
     &           eImpDef,impDefForce )
c torsionEnForces
        call allAtTorsEForce(atomXYZ,
     &           natom,quar1234List,nQuar1234,
     &           quar1234ParL,quar1234nPar,
     &           eTorsDef,torsAngForce )

 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

The deformation forces are calculated at each time step in the MD run.

6. Non bonded pair list calculation

The non bonded pair interactions are calculated for the pair list. Pair list for the central atom i is a sequence of atom numbers for atom within the radius R from the central atom. Three separate pair lists are calculated. The Van der Waals pair list(i) includes atom j if

    The pair list calculated with via the lattice algorithm:
  1. a) the atomic coordinates r1,,rN are projected on the cubic lattice, the integer coordinates of the atoms h1,,hN are obtained. The lattice size is quite small ~ 2 A, to include just one atom.

The linked list and all pairList (nnbPairLV, nnbPairLC, nnbPairLS) are calculated in the subroutine

        subroutine nonbondListVCS(rcutV,rcutC,rcutS,atomXYZ,atomQ,
     &           rbuffV,rbuffC,rbuffS,
     &           makeVdW,makeCL,makeS,
     &           natom,atomNumb,atomName,resName,chName,resNumb,
     &           nres,resNameRes,chNameRes,
     &           atomNameEx,startAtInRes,
     &           nmoveatom,moveAtomList,moveFlag,
     &           pair12List,startPairL12,nPairL12,
     &           pair13List,startPairL13,nPairL13,
     &           pair14List,startPairL14,nPairL14,
     &           nbpairListV,startnbPairLV,nnbPairLV,nnbpLVMAX,
     &           nbpairListC,startnbPairLC,nnbPairLC,nnbpLCMAX,
     &           nbpairListS,startnbPairLS,nnbPairLS,nnbpLSMAX)

fragment of code for the linked list calculation:

c distribute atoms over cells
c make linked list of atoms in cells
c headat(n) - head(incellN)
c linkList(ia) - linkedList
        do ia = 1,natom
c calculate cell numb
        ix = xyzi(1)/cellh+1
        iy = xyzi(2)/cellh+1
        iz = xyzi(3)/cellh+1
        if(ixm .lt. ix)ixm = ix
        if(iym .lt. iy)iym = iy
        if(izm .lt. iz)izm = iz
c cell number
        ncell = ix + (iy-1)*nsiz(1) + (iz-1)*nsiz(1)*nsiz(2)
        if(ncell .gt. ncell3MAX)then
        write(kanalp,*)'ERROR!:nonbondList: ncell3MAX is low !!'
        end if!
c make linked list        
        linkList(ia) =  headat(ncell) 
        headat(ncell) = ia
        end do !ia
c end of linked list calculation

The pair lists VDW and COULOMbic energy exclude 12, 13, 14 covalent bonded pairs. The Solvent model pairList 
include all 12,13, 14 pairs.
The pair list are calculated for the range respectively:
        rcutV2 = (rcutV + rbuffV)**2     ! range for List1 - 
                                                     VDWaals -  nbPairListV
        rcutV2m = (rcutV - rbuffC)**2    ! range for List2 - Coulombic twin
                                                        range - nbPairListC   
        rcutC2p = (rcutC + rbuffC)**2    ! range for List2
        rcutS2 = (rcutS + rbuffS)**2     ! range for SolvationGSList - 

see file nonbobdListVCS.f 

7. Non bonded force calculation

Van der waals forces are calculated for the non-bonded pair list nbpairListV()for atoms j within rij < RCUTV the cutoff radius for van der waals interactions. The modified potential 6-12 are used

the pair list for atom i includes atoms j > i, to count each pair interaction once. The force Fvdwi on atom i due to interaction with atoms in the pair list

The modified (smoothed) 6-12 potential prevents over-flow when atoms are too close and generates smooth driving forces to resolve clash problems between atoms in molecular dynamics simulations, see

	subroutine vdwenforceij(dij2,dij1,rij,A12,B12,evdw,fi)

The coulombic energy and forces for atom i are calculated for all pairs within the radius RCUTC.
The coulombic energy/forces for a central atom i are calculated for the classical coulombic law or as a coulombic interaction between two charges on the compensating background charge uniformly distributed within the sphere of radius RCUTC

has zero interaction energy and forces for the rij > RCUTC. This form of electrostatic interactions is better suitable to prevent energy conservation in the molecular dynamic calculation, see

        subroutine coulenforceij(var,rcutC,dij2,dij1,rij,qi,qj,ecoul,fi)

The nonbonded energy and force within short range RCUTV=R1 are calculated in the subroutine

c allAtNonBondEForce : VDW and COULOMBIC
        subroutine allAtVDWEForceR1(atomXYZ,atomQ,
     &           natom,nmoveatom,moveAtomList,
     &           nbpairListV,startnbPairLV,nnbPairLV,
     &           pair14List,startPairL14,nPairL14,
     &           nVDWtype,atomVDWtype,atomVDW12ab,
     &           rcutV,rcutC,engVDW,vdwForce,engCOULR1,coulForceR1)

for the pair list nbpairListV() and pair14List(). The last one includes all 1-4 neihgbours for which the amber force field uses the scaling factors for van der waals and coulombic interactions.
To increase performance of the van der waals energy/force calculations the table of coefficient A12, B12 for all atom types are precalculated and then right values A12/B12 for a given atom types in the pair ij are extracted from the vdw AB-parameter table

c get pointer to the AB table
        call vdw12TablePos(nVDWtype,t1,t2,t12)
        p4 = 4*t12
        A12 = atomVDW12ab(p4-3)
        B12 = atomVDW12ab(p4-2)

The long-range electrostatic forces within RCUTV < rij < RCUTC are calculated via the subroutine

        subroutine allAtVDWEForceR2(atomXYZ,atomQ,
     &           natom,nmoveatom,moveAtomList,
     &           nbpairListC,startnbPairLC,nnbPairLC,
     &           rcutR1,rcutR2,engCOULR2,coulForceR2)
c LongRamge -   RCUT1 < rij < RCUT2

The program keep separately the short-range and the long-range electrostatic energy and force.

8. Solvation energy/force calculation

The implicit solvation model - the Gaussian Shell model of Lazaridis & Karplus is used to calculate the solvation energy [POTEINS 35: 133-152, 1999]. The solvation free energy of the atom i

where sum is going over all neighbors of atom i which exclude volume Vj from the solvation volume around of the atom i. The function gi(r) describe the solvation energy density in the volume around the atom i and is approximated by the Gaussian function

The sum over all solvation forces fi is zero.
The solvation forces are calculated by subroutine

	         call SolventEnForces(natom, atomXYZ,
     &         atomName,startPairL12,nPairL12,pair12List,
     &         nbpairListS,startnbPairLS,nnbPairLS,
     &         atomSolPar, molSolEn, atomSolEn, atomSolFr)

IV. Details of MD run

An MD run is performed by subroutine

	subroutine mdRun(ntimeMX,ntime0,ntime,ntimeR1,ntimeR2,
     &                ntimeF1,ntimeF2,ntimeF3,deltat,
     &                tempTg,tauTRF,atype,optra,wtra,nwtra,cltra)	
c MD RUN propagates MDtraj from files in mdAtomXYZvel.h 
c                                     [ atomXYZ0(*),atomVel0(*) ]
c    call initMDStart(T) inits the MD start 
c                          from the INput atomXYZ(*)-->atom0XYZ(*)
c ntimeMX max number of time steps
c ntime0 - executed number of timesteps in the previous call
c ntime   executed number of timesteps in this call 
c ntimeR1, ntimeR2 - update frequency  for R1, R2 pairLists
c ntimeF1,ntimeF2 - update freq for R1=(vdw+coulR1), R2-coulR2 en/forces
c ntimeF3 -  SOLVation forces
c GeoEn/force ntimeFg=1 - standart
c deltat- timestep, temp - initial(temp) of MD run
c tempTg - target T for NTV ansemble[K]
c tauTRF - tau Relaxation Factor [ps]
c atype - ansamble type = 0/1 - NEV, NTV

The MD algorithm consist of a long loop over the time steps.
For each time step MD trajectory is propagated for the )t = 1-2 femto sec, as defined by user.

1. Pair lists

The pair lists are updated for each n-th timestep equal to ntimeR1, ntimeR2 for the short-range and for the twin-range long-range electrostatic energy calculations.

        call initNonBondList(atomXYZ0,makeVdW,makeCL,makeSL) 

2. The atomic forces

The atomic forces due to deformation of covalent structure and short-range non-bonded calculation are updated for the each ntimeF1-th time step, the long-range electrostatic are updated for the each ntimeF2-th step and solvation forces are updated for each ntimeF3-th time step.
{Note! In the current version the multiple time step for pair list update and md equation integration are equal. The general case is not tested !}

c update forces/energy
	call initAllForce(fcall,atomXYZ0,doVdWef,doCLef,doSLef,
     &              eVbondDef,vbdefForce,
     &              eVangDef,vAngdefForce,
     &              eImpDef,impDefForce,
     &              eTorsDef,torsAngForce, 
     &              engVDWR1,vdwForceR1,
     &              engCOULR1,coulForceR1,
     &              engCOULR2,coulForceR2,
     &              restr1Eng,restr1AtForce,
     &              molSolEn, atomSolEn, atomSolFr)

MD simulation can be done with a specified set of forces. The set of forces can be specified by the array fEngWF(*)

          eGeoDef   = fEngWF(1)*eVbondDef + fEngWF(2)*eVangDef 
     &                + fEngWF(3)*eImpDef + fEngWF(4)*eTorsDef
     &                + fEngWF(8)* restr1Eng
          engCOUL   = fEngWF(6)*engCOULR1 + fEngWF(7)*engCOULR2
          engPOTENT = eGeoDef + fEngWF(5)*engVDWR1 + engCOUL + 
     &                molSolEn*fEngWF(9)

3. Propogation of the trajectory

For one time step propagation of the MD trajectory is done by the subroutine

c make mdStep
	call mdTimeStepProp(nmoveatom,moveAtomList,deltat)

which uses multi step leap-frog algorithm to calculate velocities and positions at time (t+deltat).

4. Temperature control - Berendsen thermostat method

At each time step the temperature control routine performs calculation of the total kinetic energy of the moving atoms. The relaxation the average temperature of the atomic system to the specified value are give via the weak-coupling method or Berendsen method, which scale the velocity by the factor lambTR(t)

5. Trajectory writing

Trajectory is written for each nwtra time steps. The trajectory can be written for atomic positions (and for atomic velocietis) in the user specified file.

6. Blind Docking method: Exhaustive search of low-resolution binding sites with global optimization of ligand pose

Method description
Docking method is performed by subroutine runLigDock02 in the mDyn09 program procedure runLigDock02 perform docking of molecular ligand of size up to ~120 atoms.

The algorithm flow can be described as:

1) Calculation of the accessible surface of the protein. Calculation of a surface grid for probe sphere of radius ~ average atomic radius, and contact positions [bindSiteAt01(*)] with protein atoms. Calculation are done by subroutine surf_SAS04.

2) Calculation of a surface grid points for a probe ligand of radius of typical aromatic ring [benzene] gridsizeSAS ~ 3.0 A. The surface grid are calculated by clustering of surface contact positions bindSiteAt01(*) and the surface grid bindGridXYZSAS01(*) is generated. The contact score [nsasGridPoint(*)] equal to the number of contact atomic positions included in to the surface grid point bindGridXYZSAS01(*) is calculated.
The bindGridXYZSAS01(*) are sorted by descent of the contact score value nsasGridPoint(*) and presents an initial trial positions for refined docking of ligand.

3) Refined docking is performed via subroutine runLigDock01(ig,bindGridXYZSAS01loc). For each initial positions bindGridXYZSAS01(*) for ligand center.

The best docking variant for the ligand can be chosen as a file LigDockFinMMM.nnn.pdb with minimal potential energy engPOTENTLG:


1bty : benzamidine-trypsine complex

File   #LigBindGridOnSAS:         X       Y        Z     contactScore
ATOM      1  LBSt        1      16.536  26.130   8.764   11
ATOM      2  LBSt        2      29.319  14.972  16.378   11
ATOM      3  LBSt        3       6.595  15.454  32.366    9
ATOM      4  LBSt        4      28.049  26.396   3.572    9
ATOM      5  LBSt        5      37.370  14.662  29.278    8
ATOM      6  LBSt        6       9.605  28.662  39.481    7
ATOM      7  LBSt        7      18.280  35.574  15.402    7
ATOM      8  LBSt        8      30.648  34.679  44.060    7
ATOM      9  LBSt        9      34.040  33.767  21.484    7
ATOM     10  LBSt       10       5.056  19.922  18.987    6
ATOM     11  LBSt       11      25.308   5.865  13.437    6
ATOM     12  LBSt       12      13.241  31.812  30.019    6
ATOM     13  LBSt       13       6.174  15.317  15.623    6
ATOM     14  LBSt       14      15.230  11.995  39.322    6
ATOM     15  LBSt       15      42.858  27.966  33.933    6
ATOM     16  LBSt       16      39.046  14.805   5.421    5
ATOM     17  LBSt       17      24.676  37.002  14.221    5
ATOM     18  LBSt       18      39.100  25.116   6.122    5
ATOM     19  LBSt       19      25.156   6.498   5.813    5
ATOM     20  LBSt       20      14.736  13.757   2.279    5
ATOM     21  LBSt       21      35.933  31.703  11.547    5
ATOM     22  LBSt       22      45.035  21.844  22.099    5
ATOM     23  LBSt       23      12.210   8.874  28.161    5
ATOM     24  LBSt       24      11.197  11.080  32.573    5
ATOM     25  LBSt       25      25.549  16.554  -0.897    4
ATOM     26  LBSt       26      34.793   8.348  15.236    4
ATOM     27  LBSt       27      26.857   9.202  21.336    4
ATOM     28  LBSt       28      34.072  12.246  27.335    4

1) 1bty complex benzamidine on trypsine

Fig.1. Docking results for benzamidine on trypsine - 1bty complex.
A - contact Score (black square) for binding grid points vs refined potential energy of ligand binding (red diamonds).
B - minimum energy docking mode (red bonds), RMSD = 0.54 A for all non Hydrogen atoms ligand of the native binding mode.
CPK- green and violet are less favorable binding modes with low binding energy are shown in (A). CPK (pink) - native binding mode of benzamidine in 1bty.

2) 1dwb : thrombin + benzamidine complex

Fig.2 Docking results for benzamidine on thrombin.
A - Contact Score (black square) for binding grid points vs refined potential energy of ligand binding (red diamonds).
B (CPK blue) - minimum energy docking mode. Less favorable binding modes are shown - yellow, brown, green. CPK- (red) native benzamidine binding mode in 1dwb complex,

Minimum energy mode has RMSD = 0.27 A from the native binding mode of benzamidine.

3) Biotine - streptavidine complex - 1stp

Fig.3. Docking result for biotine on streptavidine , 1stp complex.
A - contact Score (black square) for binding grid points vs refined potential energy of ligand binding (red diamonds).
B - minimum energy docking mode structure of biotine - CPK, lines - native biotine in the 1stp complex.

Minimum energy mode has RMSD = 0.96 A from the native binding mode of biotine.

4) Trypsinogen/pancreatic trypsin inhibitor + Ile-Val peptide complex : 3tpi

Fig. 4. Docking result for ILE-VAL dipeptide on Trypsinogen/pancreatic trypsin inhibitor.
A - contact Score (black square) for binding grid points vs refined potential energy of ligand binding (red diamonds).
B - Lines are minimum energy docking modes of rank 1- 4 structures of ILE-VAL peptide - lines, CPK - native binding mode of biotine in the 1stp complex.

The best binding energy mode has RMSD = 0.46 A from the native binding mode of dipeptide ILE-VAL

Table 1. Energies of top ranked binding modes, and RMSD from the native binding mode.

Binding mode ePL, kcal/mole RMSD, A
Rank 1 -
-76.07 0.46
Rank2 -
-75.6 0.58
Rank3 -
-75.5 0.78
Rank4 -
-74.8 0.88

5) 1dwc complex of Human thrombin with thrombin-inhibitor MIT .

Fig. 5. Human thrombin - 296 residues; MIT - molecule includes 80 atoms.
A. Top Ranked calculated docking mode - red lines, CPK - native MIT in the native binding mode, RMSD = 0.2 A for calculated docking mode from the native.
B. 1dwc complex. Red lines is docked MIT ligand, CPK is the native mode.

6) 1hiv complex of HIV1 protease with inhibitor NOA

Fig.6. A. Two top ranked calculated binding modes of NOA in comparison with the NOA ligand in the native binding mode of 1hiv complex. CPK - native binding mode, lines (red and grey) the top ranked mode by energy of binding. The RMSD from the native are ~3.1 A for all atoms. The major difference between native and calculated modes are the orientation of one aromatic double-ring at the top of molecule NOA, the RMSD = 1.1. A over all atoms except the later aromatic system.
B. 1hiv complex of HIV1 protease with inhibitor NOA CPK - native mode, red and grey lines - are calculated modes.

7) 1hvr complex of HIV1 protease with inhibitor XK2

Fig. 7. A. Calculated binding mode of XK2, red lines, CPK - native binding mode of XK2 ligand. RMSD = 0.95 A for all atom.
B. Calculated docking mode for the ligand XK2 in complex with HIV1 protease, CPK - the native binding mode of the XK2 ligand.

8) 1hvp complex of 1HIV protease with VAC molecule inhibitor

Fig. 8. A. Calculated best binding mode of VAC is in red lines, CPK - native VAC inhibitor in the 1hvp complex; the RMSD = 0.99 A.
B. 4hvp complex, red lines is the calculated mode, CPK - the native binding mode of VAC inhibitor.

Table 1. Results of MdDock method for a set of complexes
complex Ntors RMSD, A DEgap
1) 1bty trypsin/benz 0 0.5 9.7
2) 1dwb a-thrombin/benz 0 0.5 13.3
3) 1stp streptavidine/biotin 5 0.96 29.5
4) 3tpi trypsinogen/Ile-Vla 6 0.42 10.6
5) 1dwc a-thrombin/MIT 8 0.2 10.8
6) 1hiv HIV1 protease/NOA 16 1.1/3.1 2.6
7) 1hvr HIV1 protease/XK263 8 0.95 39.1
8) 4phv HIV1 protease/VAC 15 0.9 3.4

Ntors - number of flexible torsion angles.
DEgap - energy gap between lowest energy binding mod and the next energy mode.

The developed method of blind docking show a good accuracy in prediction of the native bindig modes of flexible ligands. At the test set of 8 ligands the method shows 100% accuracy, i.e. the native binding mode are found as the mode with highest binding energy calculated in the all-atom force field.


Tamar Schlick. Molecular Modeling and simulation. Springer-Verlag, New York, 2000.
Cornell W.D., Cieplak P., Bayly C.I., Gould I.R., Mertz K.M., Ferguson D., Spellmeyer D.C., Fox T., Caldwell J.W., Kollam P.A. A second generation force field for the simulation of proteins, nucleic acids and organic molecules. J.Am.Chem.Soc. 1995: 117, p.5179-5197
Lazaridis T., Karplus M. Proteins: Structu, Funct., and Gen. 1999: 35, p.133-152