[Fwd: Re: SMD on center of mass]

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Fri Nov 03 2006 - 14:05:47 CST

(with apologies, didn't reply to the group before)


attached mail follows:


Hi Gianluca,
the code does indeed behave as documented. You should set your SMD force
to be your intended spring constant (NOT k/N); the force actually applied to each
atom will be scaled by its share of the group mass.

The reason for the confusion is that you've found the section of the
code where NAMD calculates what total force the group needs, but this
isn't where the forces are actually applied to atoms.

   modifyGroupForces()[0] = f;

toggles a boolean flag and then sets one of the group force vectors to
be f. The place this is actually applied (if I'm correct -- Jim is the
final authority on this) is around line 139 of ComputeGlobal.C (in the
version of the source tree I'm looking at; it's around the comment
"calculate forces for atoms in groups"), which
iterates over each group and over each atom in each group. As you'll see
there, the forces applied to each atom are scaled such that the overall
force is k(x-x0)**2; the force is NOT applied raw to each atom. I know
this may seem indirect, but it's necessary to properly and efficiently
collect and distribute the forces for paralell calculations.

I hope this is helpful; please let me know if I've left something
unclear.

Peter

On Thu, Nov 02, 2006 at 07:27:41PM -0800, Gianluca Interlandi wrote:
> I have found part of the answer myself by looking into the code
> (GlobalMasterSMD.C):
>
> GlobalMasterSMD::GlobalMasterSMD(BigReal spring_constant, BigReal
> velocity,
> const Vector direction, int output_frequency,
> int first_timestep, const char *filename) {
>
> ...
>
> BigReal diff = (curcm - cm)*moveDir;
> Force f = k*(moveVel*currentTime - diff)*moveDir;
> modifyGroupForces()[0] = f;
>
>
> The force "f" is calculated based on the displacement of the center of
> mass. However, from the line "modifyGroupForces()[0] = f;" it seems that
> the calculated force is applied to each SMD targeted atom. Does this mean
> that the more targeted atoms I have the higher will be the force applied
> to the whole group of SMD atoms? Conversely, in the documentation for TMD
> I read that the spring constant "k" is scaled down by the number N of
> targeted atoms. This does not seem the case for SMD. Am I right? What does
> the line "modifyGroupForces()[0] = f;" exactly do? Does it apply "f" to
> each atom or is the force split equally among the targeted atoms? In the
> first case, does it make sense to set SMDk = k/N, where "k" is my
> "intended" spring constant and "N" the number of targeted atoms?
>
> I would appreciate any help to get this clear.
>
> Gianluca
>
>
> On Thu, 2 Nov 2006, Gianluca Interlandi wrote:
>
> > I'm sorry. This is the second e-mail I'm posting about this issue. Again a
> > question about what I read in the documentation:
> >
> > > The center of mass of the SMD atoms will be harmonically constrained with
> > > force constant"k" (SMDk) to move with velocity "v" (SMDVel) in the
> > > direction "n" (SMDDir).
> >
> > SMD thus results in the following potential being applied to the system:
> >
> > U(r_1, r_2, ..., t) = 1/2 k [vt - (R(t) - R_0) . n ]^2.
> >
> > What is the SYSTEM? Is it the whole protein or just the SMD atoms?
> >
> > Thanks a lot,
> >
> > Gianluca
> >
> > -----------------------------------------------------
> > Dr. Gianluca Interlandi gianluca_at_u.washington.edu
> > +1 (206) 685 4435
> > +1 (206) 714 4303
> > http://biocroma.unizh.ch/gianluca/
> >
> > Postdoc at the Department of Bioengineering
> > at the University of Washington, Seattle WA U.S.A.
> > -----------------------------------------------------
> >
>
> -----------------------------------------------------
> Dr. Gianluca Interlandi gianluca_at_u.washington.edu
> +1 (206) 685 4435
> +1 (206) 714 4303
> http://biocroma.unizh.ch/gianluca/
>
> Postdoc at the Department of Bioengineering
> at the University of Washington, Seattle WA U.S.A.
> -----------------------------------------------------

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:42:47 CST