Help regarding PMF with SMD

From: Katherine Parra Pulido (kparra_at_mail.usf.edu)
Date: Thu Aug 01 2013 - 10:06:03 CDT

Hello I'm trying to calculate Potential Mean Force between a complex of a
receptor and its ligand, using Steered Molecular Dynamics. The problem is
that I've been getting this error message:
TCL: Running for 1000000 steps
------------- Processor 0 Exiting: Called CmiAbort ------------
Reason: FATAL ERROR: floating-point value too large to represent
    while executing
"expr $k*($c1x-$r1x)"
    (procedure "calcforces" line 21)
    invoked from within
"calcforces"

*****************************************************************************************************************************
Here is the conf file:

#############################################################
## JOB DESCRIPTION ##
#############################################################

# SMD simulation of 4-BP pulling from DHP
# Constant temperature

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################

structure dhp.psf
coordinates dhp.pdb
outputName dhp0

set temperature 300

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters par_all27_prot_lipid.inp
temperature $temperature

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 12.
switching on
switchdist 10.
pairlistdist 14.

# Integrator Parameters
timestep 2.0 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10

# Constant Temperature Control
langevin on ;# do langevin dynamics
langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
langevinTemp $temperature
langevinHydrogen no ;# don't couple langevin bath to hydrogens

# Output
binaryoutput no
dcdfreq 500 ;# 500steps = every 1ps
outputEnergies 100

# Particle mesh Ewald

PME yes
PMEGridSpacing 1.0
PMEGridSizeX 128
PMEGridSizeY 128
PMEGridSizeZ 216

#############################################################
## PBC PARAMETERS ##
#############################################################

#Periodic Boundary Conditions [if a water box is involved]

cellBasisVector1 93.3 0.0 0.0 ;# x-value for box dimension

cellBasisVector2 0.0 93.3 0.0 ;# y-value for box dimension

cellBasisVector3 0.0 0.0 138.9 ;# z-value for box dimension

cellOrigin 0.0 0.0 1.4 ;# x-center, y-center, z-center

############################################################

#############################################################
## EXTRA PARAMETERS ##
#############################################################

# Tcl interface
tclForces on
tclForcesScript dhp0.tcl
run 1000000 ;# 2 ns

*************************************************************************************************************************************
And tcl script:
# $Id: smd.tcl,v 1.2 2005/02/18 18:07:11 mbach Exp $

# Atoms selected for force application

set targets1 {}
set targets2 {}
set inStream [open dhp.pdb r]
foreach line [split [read $inStream] "\n"] {

 set type [string trim [string range $line 0 5]]
 set name [string trim [string range $line 12 15]]
 set resid [string trim [string range $line 22 25]]
 set beta [string trim [string range $line 60 65]]
 set segname [string trim [string range $line 72 75]]

 if { $beta == 1.00 } {
   lappend targets1 "$segname $resid $name"
                         }
 if { $beta == 2.00 } {
   lappend targets2 "$segname $resid $name"

}
}
close $inStream

set atoms1 {}
foreach target1 $targets1 {
  foreach {segname resid atom} $target1 { break
}

  set atomindex1 [atomid $segname $resid $atom]
  lappend atoms1 $atomindex1
                           }

set atoms2 {}
foreach target2 $targets2 {
  foreach {segname resid atom} $target2 { break
}

  set atomindex2 [atomid $segname $resid
$atom]
  lappend atoms2
$atomindex2

}

set a1 [addgroup $atoms1]
set a2 [addgroup $atoms2]

# set the output frequency, initialize the time counter
set Tclfreq 50
set t 0

# contraint points

set c1x -0.25
set c1y 0.30
set c1z -3.82

set c2x -1.30
set c2y 2.53
set c2z 58.84

# force constant (kcal/mol/A^2)
set k 10.0

# pulling velocity (A/timestep)
set v 0.00002

set outfilename dhp0.out
open $outfilename w

proc calcforces {} {

  global Tclfreq t k v a1 a2 c1x c1y c1z c2x c2y c2z outfilename

  # get coordinates

  loadcoords coordinate

  set r1 $coordinate($a1)
  set r1x [lindex $r1 0]
  set r1y [lindex $r1 1]
  set r1z [lindex $r1 2]

  set r2 $coordinate($a2)
  set r2x [lindex $r2 0]
  set r2y [lindex $r2 1]
  set r2z [lindex $r2 2]

  # calculate forces

  set f1x [expr $k*($c1x-$r1x)]
  set f1y [expr $k*($c1y-$r1y)]
  set f1z [expr $k*($c1z-$r1z)]
  lappend f1 $f1x $f1y $f1z

  set f2x [expr $k*($c2x-$r2x+$v*$t)]
  set f2y [expr $k*($c2y-$r2y+$v*$t)]
  set f2z [expr $k*($c2z-$r2z)]
  lappend f2 $f2x $f2y $f2z

  # apply forces

  addforce $a1 $f1
  addforce $a2 $f2

  # output

  set foo [expr $t % $Tclfreq]
  if { $foo == 0 } {
      set outfile [open $outfilename a]
      set time [expr $t*2/1000.0]
      puts $outfile "$time $r2x $r2y $r2z $f2x $f2y $f2z"
      close $outfile
  }
  incr t
  return
}

Hope you can help, this is my time using SMD.

KP

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:23:31 CST