# The IDs of the atoms defining the restraint
 
  set aid15  [atomid NBDA  309 OY ]
  set aid16  [atomid NBDA  312 N  ]
  set aid17  [atomid NBDB  309 OY ]
  set aid18  [atomid NBDB  312 N  ]
  set aid21  [atomid IC1A   91 OY]
  set aid22  [atomid IC1A   94 N ]
  set aid23  [atomid IC1B   91 OY]
  set aid24  [atomid IC1B   94 N]
  set aid25  [atomid IC2A  188 OY]
  set aid26  [atomid IC2A  191 N ]
  set aid27  [atomid IC2B  188 OY]
  set aid28  [atomid IC2B  191 N]
  set aid29  [atomid IC1A   127 O ]
  set aid30  [atomid IC1A   130 NT]
  set aid31  [atomid IC1B   127 O ]
  set aid32  [atomid IC1B   130 NT]
  set aid33  [atomid IC2A   234 O ]
  set aid34  [atomid IC2A   237 NT]
  set aid35  [atomid IC2B   234 O ]
  set aid36  [atomid IC2B   237 NT]



 
# The "spring constant" for the harmonic constraint
  set k  10
  addatom $aid15
  addatom $aid16
  addatom $aid17
  addatom $aid18
  addatom $aid21
  addatom $aid22
  addatom $aid23
  addatom $aid24
  addatom $aid25
  addatom $aid26
  addatom $aid27
  addatom $aid28

  addatom $aid29
  addatom $aid30
  addatom $aid31
  addatom $aid32
  addatom $aid33
  addatom $aid34
  addatom $aid35
  addatom $aid36


  proc calcforces {} {
  
    global aid15 aid16 aid17 aid18 aid21 aid22 aid23 aid24 aid25 aid26 aid27 aid28 aid29 aid30
    
    global aid31 aid32 aid33 aid34 aid35 aid36  k


    loadcoords p
    
    # Calculate the current distance
    set MG13 [getbond $p($aid15) $p($aid16)]
    set MG14 [getbond $p($aid17) $p($aid18)]
    set MG15 [getbond $p($aid21) $p($aid22)]
    set MG16 [getbond $p($aid23) $p($aid24)]
    set MG17 [getbond $p($aid25) $p($aid26)]
    set MG18 [getbond $p($aid27) $p($aid28)]
    set MG19 [getbond $p($aid29) $p($aid30)]
    set MG20 [getbond $p($aid31) $p($aid32)]
    set MG21 [getbond $p($aid33) $p($aid34)]
    set MG22 [getbond $p($aid35) $p($aid36)]

 # Calculate the current force/distance
    set MG13 [expr ($MG13-3.0)*$k/$MG13]
    set MG14 [expr ($MG14-3.0)*$k/$MG14]
    set MG15 [expr ($MG15-3.0)*$k/$MG15]
    set MG16 [expr ($MG16-3.0)*$k/$MG16]
    set MG17 [expr ($MG17-3.0)*$k/$MG17]
    set MG18 [expr ($MG18-3.0)*$k/$MG18]
    set MG19 [expr ($MG19-3.0)*$k/$MG19]
    set MG20 [expr ($MG20-3.0)*$k/$MG20]
    set MG21 [expr ($MG21-3.0)*$k/$MG21]
    set MG22 [expr ($MG22-3.0)*$k/$MG22]


  set g14  [vecsub $p($aid16) $p($aid15)]
  set g15  [vecsub $p($aid18) $p($aid17)]
  set g16  [vecsub $p($aid22) $p($aid21)]
  set g17  [vecsub $p($aid24) $p($aid23)]
  set g18  [vecsub $p($aid26) $p($aid25)]
  set g19  [vecsub $p($aid28) $p($aid27)]
  set g20  [vecsub $p($aid29) $p($aid30)]
  set g21  [vecsub $p($aid31) $p($aid32)]
  set g22  [vecsub $p($aid33) $p($aid34)]
  set g23  [vecsub $p($aid35) $p($aid36)]


    
    # The force to be applied on each atom is proportional to its
    # corresponding gradient

    addforce $aid15 [vecscale $g14 $MG13]
    addforce $aid17 [vecscale $g15 $MG14]
    addforce $aid21 [vecscale $g16 $MG15]
    addforce $aid23 [vecscale $g17 $MG16]    
    addforce $aid25 [vecscale $g18 $MG17]
    addforce $aid27 [vecscale $g19 $MG18]    

    addforce $aid30 [vecscale $g20 $MG19]
    addforce $aid32 [vecscale $g21 $MG20]
    addforce $aid34 [vecscale $g22 $MG21]
    addforce $aid36 [vecscale $g23 $MG22]

  }