Merge pull request #50 from rvhonorato/flexref

Implement the semi flexible refinement module (`flexref`)

examples/rigid_body_docking/rigid_body_docking.toml → examples/simple_docking/simple_docking.toml

@@ -3,7 +3,7 @@
 
 [input]
 # in which order the steps must be executed?
-order = ["topology", "rigidbody"]
+order = ["topology", "rigidbody", "flexref"]
 
 # directory in which the scoring will be done
 project_dir = "run1"
@@ -24,7 +24,12 @@ autohis = true
 ########################
 [stage.rigidbody]
 ambig = 'ambig.tbl'
-sampling = 100
+sampling = 2
+
+########################
+[stage.flexref]
+ambig = 'ambig.tbl'
+cool1_steps = 500
 
 # ====================================================================
 

src/haddock/cns/util.py

@@ -3,8 +3,6 @@
 from haddock.pdbutil import PDBFactory
 from haddock.mathutil import RandomNumberGenerator
 from haddock.defaults import Default
-from haddock.ontology import Format
-# from haddock.cns.util import load_workflow_params, prepare_input
 
 RND = RandomNumberGenerator()
 
@@ -227,7 +225,7 @@ def prepare_multiple_input(pdb_input_list, psf_input_list):
     for pdb in pdb_input_list:
         input_str += f'coor @@{pdb}{linesep}'
 
-    ncomponents = len(pdb_input_list)
+    ncomponents = len(psf_input_list)
     input_str += f'eval ($ncomponents={ncomponents}){linesep}'
 
     seed = RND.randint(100, 999)

src/haddock/modules/__init__.py

@@ -63,8 +63,7 @@ def _load_previous_io(self):
 
     def previous_path(self):
         if self.order > 1:
-            return (self.path.resolve().parent.absolute() /
-                    f"{MODULE_PATH_NAME}{self.order-1}")
+            return (self.path.resolve().parent.absolute() / self.stream['input']['order'][self.order-1])
         if self.order == 1:
             return self.path.resolve().parent.absolute() / TOPOLOGY_PATH
 

src/haddock/modules/flexref/cns/calc_free-ene.cns

@@ -0,0 +1,24 @@
+! calc_free-ene.cns
+!    Calculate the total energy of the separated components
+!
+! Script uses the store1 and store2 atom object, and COMP, REF and RMSD object
+! ***********************************************************************
+! * Copyright 2003-2018 Alexandre Bonvin, Utrecht University.           *
+! * All rights reserved.                                                *
+! * This code is part of the HADDOCK software and governed by its       *
+! * license. Please see the LICENSE file that should have been included *
+! * as part of this package.                                            *
+! ***********************************************************************
+!
+@RUN:scale_intra_only.cns
+
+evaluate ($eintfree = 0.0)
+@RUN:flex_segment_back.cns
+fix sele=(((attr store5 = 0) or resn ANI or resn DAN or resn XAN or resn SHA) and not name H* and not (resn WAT or resn HOH or resn TIP*)) end
+
+minimize powell nstep=100 drop=10.0 nprint=25 end
+fix sele=(not all) end
+energy end
+evaluate ($eintfree = $bond + $angl + $impr + $dihe + $vdw + $elec)
+
+display FREE MOLECULES INTERNAL ENERGY = $eintfree

src/haddock/modules/flexref/cns/centroids_create.cns

@@ -0,0 +1,56 @@
+!module (ncomponents; Toppar;)
+    {Create and initialize Dummy residues for each segid
+
+    Parameters
+    ----------
+    ncomponents : int
+        Number of components/segids
+
+    Toppar : data-structure
+        The Toppar data-structure
+
+    Side-effects
+    ------------
+    A DUM residue is appended to each segid and placed at the origin
+    }
+
+    ! define topology of DUM residue
+    topology
+        mass DD 100.00
+        residue DUM group
+            atom DUM type=DD charge=0.000 end
+            atom MAP type=DD charge=0.000 end
+        end
+    end
+    ! define non-bonded parameters
+    parameter
+        nonbonded DD 0.001 0.001 0.001 0.001 end
+    end
+
+    ! place dummy atom of dummy residue in the center of each chain
+    evaluate($nchain1 = 0)
+    while ($nchain1 < &ncomponents) loop nloopdum
+        evaluate($nchain1 = $nchain1 + 1)
+
+        ! create the dummy residue
+        segment
+            name="TMP1"
+            chain 
+                sequence "DUM" end
+            end
+        end
+
+        ! set the residue id of dummy to last residue + 1
+        show max(decode(resid)) (segid &Toppar.prot_segid_$nchain1)
+        evaluate($maxresid = $result)
+
+        ! change the segid of the dummy residue
+        do (segid = &Toppar.prot_segid_$nchain1) (segid TMP1)
+
+        ! set residue id and coordinate of the dummy residue
+        do (x = 0) (segid &Toppar.prot_segid_$nchain1 and resn DUM)
+        do (y = 0) (segid &Toppar.prot_segid_$nchain1 and resn DUM)
+        do (z = 0) (segid &Toppar.prot_segid_$nchain1 and resn DUM)
+        do (resid = encode($maxresid + 1)) (segid &Toppar.prot_segid_$nchain1 and resn DUM)
+
+    end loop nloopdum

src/haddock/modules/flexref/cns/centroids_initialize.cns

@@ -0,0 +1,17 @@
+{Places DUM-residue on the center of each segid}
+
+eval($nchain1 = 0)
+while ($nchain1 < $data.ncomponents) loop nloop1
+    eval($nchain1 = $nchain1 + 1)
+
+    show ave(x) (segid $Toppar.prot_segid_$nchain1 and (not name H*) and (not resn DUM))
+    eval($center.x = $RESULT)
+    show ave(y) (segid $Toppar.prot_segid_$nchain1 and (not name H*) and (not resn DUM))
+    eval($center.y = $RESULT)
+    show ave(z) (segid $Toppar.prot_segid_$nchain1 and (not name H*) and (not resn DUM))
+    eval($center.z = $RESULT)
+
+    do (x = $center.x) (segid $Toppar.prot_segid_$nchain1 and resn DUM)
+    do (y = $center.y) (segid $Toppar.prot_segid_$nchain1 and resn DUM)
+    do (z = $center.z) (segid $Toppar.prot_segid_$nchain1 and resn DUM)
+end loop nloop1

src/haddock/modules/flexref/cns/centroids_set_restraints.cns

@@ -0,0 +1,70 @@
+{Sets the centroid based distance restraints}
+
+if ($iteration = 0) then
+
+  noe
+    ! define a new class for the centroid based distance restraints
+    class       centroid 
+    averaging   centroid sum 
+    potential   centroid soft
+    scale       centroid $Data.centroids.kscale
+    sqconstant  centroid 1.0
+    sqexponent  centroid 2
+    soexponent  centroid 1
+    rswitch     centroid 1.0
+    sqoffset    centroid 0.0
+    asymptote   centroid 2.0
+    msoexponent centroid 1
+    masymptote  centroid -0.1
+    mrswitch    centroid 1.0
+
+    ! set the restraints between each DUM and MAP atom for each chain. in case
+    ! of ambiguous restraints, each DUM atom is assigned with all MAP atoms
+    evaluate($nchain1 = 0)
+    while ($nchain1 < $Data.ncomponents) loop distloop
+        evaluate($nchain1 = $nchain1 + 1)
+
+        if ($Data.centroids.ambi_$nchain1 = true) then
+            assign (name MAP)
+                   ((name DUM) and (segid $Toppar.prot_segid_$nchain1)) 0.0 0.0 0.0
+        else
+            assign ((name MAP) and (segid $Toppar.prot_segid_$nchain1))
+                   ((name DUM) and (segid $Toppar.prot_segid_$nchain1)) 0.0 0.0 0.0
+        end if
+    end loop distloop
+  end
+end if
+
+if ($iteration > 0) then
+
+! define centroid based distance restraints for expand/refine protocol
+
+  noe
+    ! define a new class for the centroid based distance restraints
+    class       centroid 
+    averaging   centroid center
+    potential   centroid soft
+    scale       centroid $Data.centroids.kscale
+    sqconstant  centroid 1.0
+    sqexponent  centroid 2
+    soexponent  centroid 1
+    rswitch     centroid 1.0
+    sqoffset    centroid 0.0
+    asymptote   centroid 2.0
+    msoexponent centroid 1
+    masymptote  centroid -0.1
+    mrswitch    centroid 1.0
+
+    ! set the restraints between the heavy atoms of each chain and the
+    ! initially stored center of mass of each chain
+    evaluate($nchain1 = 0)
+    while ($nchain1 < $Data.ncomponents) loop distloop
+        evaluate($nchain1 = $nchain1 + 1)
+
+        assign (not name H* and not name DUM and segid $Toppar.prot_segid_$nchain1)
+               (name DUM and segid $Toppar.prot_segid_$nchain1) 0.0 0.0 0.0
+    end loop distloop
+
+  end
+
+end if

src/haddock/modules/flexref/cns/charge-beads-interactions.cns

@@ -0,0 +1,33 @@
+for $id1 in id ( (resn SER or resn THR or resn ASN or resn GLN) and name BB ) loop c1beads
+
+   show (segid) (id $id1)
+   evaluate ($cgsegid=$result)
+   show (resid) (id $id1)
+   evaluate ($cgresid=$result)
+
+   igroup
+     ! turn off all vdw interactions between charged beads and all beads
+     interaction (segid $cgsegid and resid $cgresid and name SCD1) (all) weight vdw 0.0 end
+     interaction (segid $cgsegid and resid $cgresid and name SCD2) (all) weight vdw 0.0 end
+
+     ! turn off all vdw and elec interactions between the charged beads within one residue
+     interaction (segid $cgsegid and resid $cgresid and name SCD1)
+                 (segid $cgsegid and resid $cgresid and name SCD2) weight * 1.0 vdw 0.0 elec 0.0 end
+   end
+
+end loop c1beads
+
+for $id1 in id ( (resn LYS or resn ARG or resn ASP or resn GLU) and name BB ) loop c2beads
+
+   show (segid) (id $id1)
+   evaluate ($cgsegid=$result)
+   show (resid) (id $id1)
+   evaluate ($cgresid=$result)
+
+   igroup
+     ! turn off all vdw interactions between charged beads and all beads
+     interaction (segid $cgsegid and resid $cgresid and name SCD1) (all) weight vdw 0.0 end
+   end
+
+end loop c2beads
+

src/haddock/modules/flexref/cns/cm-restraints.cns

@@ -0,0 +1,142 @@
+! cm-restraints.cns
+!    Define center-of-mass restraints between molecules
+!
+! ***********************************************************************
+! * Copyright 2003-2018 Alexandre Bonvin, Utrecht University.           *
+! * All rights reserved.                                                *
+! * This code is part of the HADDOCK software and governed by its       *
+! * license. Please see the LICENSE file that should have been included *
+! * as part of this package.                                            *
+! ***********************************************************************
+!
+!define center of mass restraints between all molecules
+!using distance restraints between CA, BB or N1 atoms with center averaging
+  set echo=on message=on end
+
+  evaluate ($ncount = 0)
+
+  ! store original coordinates 
+  do (refx = x) (all)
+  do (refy = y) (all)
+  do (refz = z) (all)
+
+  while ($ncount < $data.ncomponents) loop nloop1
+    evaluate ($ncount = $ncount +1)
+    evaluate ($dim_$ncount = 0.0)
+
+    !orient molecule
+    coor orient sele=(segid $Toppar.prot_segid_$ncount) end
+
+    ! find dimensions
+    show max (x) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    evaluate ($xdim = $result)
+    show max (y) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    evaluate ($ydim = $result)
+    show max (z) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    evaluate ($zdim = $result)
+    show min (x) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    evaluate ($xdim = $xdim - $result)
+    show min (y) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    evaluate ($ydim = $ydim - $result)
+    show min (z) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    evaluate ($zdim = $zdim - $result)
+
+    evaluate ($corr = max($xdim,$ydim))
+    evaluate ($corr = max($corr,$zdim))
+
+    if  ($data.cmtight eq false) then
+      ! use average of all three dimensions + 10A
+      evaluate ($dim_$ncount = ($xdim + $ydim + $zdim)/6.0)
+    else
+      ! use only the average of the smallest two dimensions
+      evaluate ($dim_$ncount = ($xdim + $ydim + $zdim - $corr)/4.0)
+    end if
+
+    if ($toppar.dna_$ncount = true) then
+      ! Check first if not protein-DNA complex
+      do (store6 = 0) (all)
+      do (store6 = 1) (segid $Toppar.prot_segid_$ncount and (name CA or name BB))
+      show sum (store6) (all)
+      if ($result eq 0) then
+        ! we are dealing with a DNA - set dimension to 0
+        evaluate ($dim_$ncount = 0.0)
+      end if
+    end if
+
+    do (store6 = 0) (all)
+    do (store6 = 1) (segid $Toppar.prot_segid_$ncount and (name CA or name BB or name N1))
+    show sum (store6) (all)
+    if ($result eq 0) then
+      ! we are dealing with a ligand set dimension to 0
+      evaluate ($dim_$ncount = 0.0)
+    end if
+
+  end loop nloop1
+
+  ! restore original coordinates
+  do (x = refx) (all)
+  do (y = refy) (all)
+  do (z = refz) (all)
+
+  eval($nchain = 0)
+  do (store9 = 0) (all)
+  do (store9 = 1) (name CA or name BB or name N1)
+  while ($nchain < $data.ncomponents) loop nloop0
+    eval($nchain = $nchain + 1)
+    show sum (store9) (segid $Toppar.prot_segid_$nchain)
+    if ($result < 3) then
+      evaluate ($selat$nchain = 0 )
+    else
+      evaluate ($selat$nchain = 1 )
+    end if
+  end loop nloop0
+
+  eval($nchain1 = 0)
+  noe
+    class contact
+    while ($nchain1 < $data.ncomponents) loop nloop1
+      eval($nchain1 = $nchain1 + 1)
+      eval($nchain2 = $nchain1 )
+      if ($Toppar.shape_$nchain1 eq false) then
+        while ($nchain2 < $data.ncomponents) loop nloop2
+          eval($nchain2 = $nchain2 + 1)
+          if ($Toppar.shape_$nchain2 eq false) then
+            if  ($data.cmtight eq false) then
+              eval($cm_dist = $dim_$nchain1 + $dim_$nchain2)
+            else
+              eval($cm_dist = ($dim_$nchain1 + $dim_$nchain2)/2 )
+            end if
+            if ($selat$nchain1 = 1) then
+              if ($selat$nchain2 = 1) then
+                assign (segid $Toppar.prot_segid_$nchain1 and ( name CA or name BB or name N1 ))
+                   (segid $Toppar.prot_segid_$nchain2 and ( name CA or name BB or name N1 )) $cm_dist $cm_dist 1.0
+              else
+                assign (segid $Toppar.prot_segid_$nchain1 and ( name CA or name BB or name N1 ))
+                   (segid $Toppar.prot_segid_$nchain2) $cm_dist $cm_dist 1.0
+              end if
+            else
+              if ($selat$nchain2 = 1) then
+                assign (segid $Toppar.prot_segid_$nchain1)
+                   (segid $Toppar.prot_segid_$nchain2 and ( name CA or name BB or name N1 )) $cm_dist $cm_dist 1.0
+              else
+                assign (segid $Toppar.prot_segid_$nchain1)
+                   (segid $Toppar.prot_segid_$nchain2) $cm_dist $cm_dist 1.0
+              end if
+            end if
+          end if
+        end loop nloop2
+      end if
+    end loop nloop1
+
+    averaging  contact center
+    scale      contact $Data.kcont
+    sqconstant contact 1.0
+    sqexponent contact 2
+    soexponent contact 1
+    rswitch    contact 1.0
+    sqoffset   contact 0.0
+    asymptote  contact 2.0
+    msoexponent contact 1
+    masymptote  contact -0.1
+    mrswitch    contact 1.0
+  end