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