Re: tcl force

From: Seungho Choe (sec53_at_pitt.edu)
Date: Sat Oct 17 2009 - 02:02:47 CDT

Hi JC,

According to the manual, I don't need to explicitly add the energy term.
Am I correct ?
If I apply a force by using "addforce" to the center of mass, the force
will be distributed to each atom.
What I saw in one of previous posts is that all atoms will move in the
same direction.
However, it will be more reasonable to me if the atoms move in random
directions.

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
>>
>> 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"
>>