# RE: tcl force

From: JC Gumbart (gumbart_at_ks.uiuc.edu)
Date: Sat Oct 17 2009 - 00:56:45 CDT

Yes, you can. But it will only work in 1D as you've proposed, e.g., not n =
(1,1,1).

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
>
> Best,
> Seungho
>
> ************************************************************************
> set numatoms 63230
>
> set atoms {}
> for {set i 63207} { \$i <= \$numatoms} { incr i} {
> lappend atoms \$i
> }
>
> print "atoms:\$atoms"
>
> foreach atom \$atoms {
> }
>
>
>
> # Set spring constant
> set k 7.0
> set z0 -3
>
> proc calcforces {} {
>
> global atoms groupid numatoms k z0
>
> #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"
>