RE: tcl force

From: JC Gumbart (gumbart_at_ks.uiuc.edu)
Date: Sun Oct 18 2009 - 08:57:10 CDT

>From my reading of the manual, you do have to manually add the energy, in
addition to the force.

 

I think your tcl script is fine, except it's more complicated than
necessary. The force is just F = -k * (z-z0); you don't need to set
conditions on r0.

 

 

From: owner-namd-l_at_ks.uiuc.edu [mailto:owner-namd-l_at_ks.uiuc.edu] On Behalf
Of Seungho Choe
Sent: Saturday, October 17, 2009 1:08 PM
To: JC Gumbart
Cc: namd-l_at_ks.uiuc.edu
Subject: Re: namd-l: tcl force

 

Hi JC and NAMD users,

I'm trying to compare both methods: tcl force & SMD.
Could you please look at my tcl script again and tell me whether it is
correct ?
I wonder whether I'm missing some factors in the script. Thanks a lot.

Best,
Seungho

JC Gumbart wrote:

Yes, you can. But it will only work in 1D as you've proposed, e.g., not n =
(1,1,1).
 
Regarding your initial question, see the addenergy term in the manual:
http://www.ks.uiuc.edu/Research/namd/2.7b1/ug/node43.html
 
 
-----Original Message-----
From: owner-namd-l_at_ks.uiuc.edu [mailto:owner-namd-l_at_ks.uiuc.edu] On Behalf
Of Seungho Choe
Sent: Saturday, October 17, 2009 12:21 AM
To: Seungho Choe
Cc: namd-l_at_ks.uiuc.edu
Subject: Re: namd-l: tcl force
 
Hi,
 
I wonder whether I can use SMD with v=0, n= (0,0,1) to restrain the
center of mass of a molecule in z direction ?
It seems that in this case I don't need to worry about the energy term
that I mentioned in the previous my post. Could someone clarify this ?
Apology if this issue was already discussed before.
 
Best,
Seungho
 
 
 
 
Seungho Choe wrote:
  

Dear all,
 
I'm trying to restrain the center of mass (only z position) of a
molecule. Below is my tcl script. What I'm wondering is whether a new
energy term from the added force is applied correctly to the total
energy. Do I need to explicitly add the energy term in the script ? Is
there any way to check this ? Any comments are welcome. Thank you in
advance.
 
Best,
Seungho
 
************************************************************************
set numatoms 63230
 
set atoms {}
for {set i 63207} { $i <= $numatoms} { incr i} {
   lappend atoms $i
}
 
print "atoms:$atoms"
 
foreach atom $atoms {
   addatom $atom
}
 
set groupid [addgroup $atoms]
 
 
# Set spring constant
set k 7.0
set z0 -3
 
proc calcforces {} {
 
global atoms groupid numatoms k z0
loadcoords p
 
#print "$p($groupid)"
 
set z1 [lindex $p($groupid) 2]
 
print "z comp:$z1"
 
set r01 [expr $z1-$z0]
#print "r01:$r01"
 
if {$r01 > 0} {
   set n01 {0 0 -1}
} elseif {$r01 < 0} {
   set n01 {0 0 +1}
}
 
# Calculate the "force"
set force [expr abs($k*$r01)]
 
set forcez [vecscale $force $n01]
 
print "forcez:$forcez"
 
addforce $groupid $forcez
 
}
 
 
    

 
  

 

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:51:35 CST