--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/01/19 21:17:39	965
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/01/26 21:45:03	983
@@ -174,7 +174,6 @@ void SimSetup::makeMolecules(void){
   bend_set* theBends;
   torsion_set* theTorsions;
 
-
   //init the forceField paramters
 
   the_ff->readParams();
@@ -182,6 +181,9 @@ void SimSetup::makeMolecules(void){
 
   // init the atoms
 
+  double phi, theta, psi;
+  double sux, suy, suz;
+  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
   double ux, uy, uz, u, uSqr;
 
   for (k = 0; k < nInfo; k++){
@@ -218,10 +220,34 @@ void SimSetup::makeMolecules(void){
           info[k].n_oriented++;
           molInfo.myAtoms[j] = dAtom;
 
-          ux = currentAtom->getOrntX();
-          uy = currentAtom->getOrntY();
-          uz = currentAtom->getOrntZ();
+          // Directional Atoms have standard unit vectors which are oriented
+          // in space using the three Euler angles.  We assume the standard
+          // unit vector was originally along the z axis below.
+
+          phi = currentAtom->getEulerPhi();
+          theta = currentAtom->getEulerTheta();
+          psi = currentAtom->getEulerPsi();
+            
+          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
+          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
+          Axz = sin(theta) * sin(psi);
+          
+          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
+          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
+          Ayz = sin(theta) * cos(psi);
+          
+          Azx = sin(phi) * sin(theta);
+          Azy = -cos(phi) * sin(theta);
+          Azz = cos(theta);
+
+          sux = 0.0;
+          suy = 0.0;
+          suz = 1.0;
 
+          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
+          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
+          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
+
           uSqr = (ux * ux) + (uy * uy) + (uz * uz);
 
           u = sqrt(uSqr);