Hello Discussion Dwellers,

I'm writing a script that places a sorbant molecule 10 Angstroms above a surface in a simulation box that is periodic x,y and non-periodic in z. Ideally, I would like to give the molecule a shove in the direction of the surface at the beginning of the simulation. I adopted a simplistic approach to do this by assigning a velocity to each atom in the sorbant molecule a la:

       foreach my \$atom (@{\$sorbant->Atoms})

       {

             \$atom->Velocity = Point ( X => 0, Y => 0, Z => -0.001 );

        }

and then running a brief test simulation:

        my \$forcite = Modules->Forcite;

        \$forcite->ChangeSettings([Ensemble3D=> "NVT",WriteForces=> "no",InitialVelocities => "Current",NumberOfSteps => 100,TrajectoryFrequency => 20,CurrentForcefield => "COMPASS26",TimeStep => 0.5,Thermostat => "Nose"]);

         my \$results= Modules->Forcite->Dynamics->Run(\$start_structure);

This routine doesn't apparently work, because the 'velocity' field in properties explorer shows that the velocity of the specified atoms is not what I assigned them to be- although the simulation seems to run without problem. (I did verify that the foreach loop above is looping through the correct atoms.)

I'm not totally sure what the velocity field in propert explorer tells me, so I added another foreach loop that prints out the velocities of these atoms:

foreach my \$atom (@{\$sorbant->Atoms})

  {

   my \$vel = \$atom->Velocity;

   print "\$vel is this\n";

  }

And the output appeared to print a non-decimal number which I assume is the memory location for the actual velocity:

ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this
ISDPoint3d [0x1a037538] is this

I'm not sure if this last bit of information is uesful.

Anyway, here are my two questions: (1) Any idea as to what I'm not doing correctly in my script? (2) Is it sufficient to assign a reasonable velocity to all atoms (i.e., consistent with the Boltzmann dist at the simulation temperature of 298K) and hope that the SHAKE algorithm holds the molecule together? If not, any quick way to assign net velocity to a set of atoms?

Thanks,

Vanessa

PS I've attached the full script just in case it helps with diagnosis.