utils/sysbuilder/MoLocator.cpp

#include <iostream>

#include <cstdlib>
#include <cmath>

#include "simError.h"

#include "MoLocator.hpp"


MoLocator::MoLocator( MoleculeStamp* theStamp ){

  myStamp = theStamp;
  nAtoms = myStamp->getNAtoms();

  myCoords = NULL;
  
  calcRefCoords();
}

MoLocator::~MoLocator(){
  
  if( myCoords != NULL ) delete[] myCoords;
}

void MoLocator::placeMol( double pos[3], double A[3][3], Atom** atomArray,
                          int atomIndex, SimState* myConfig ){

  int i,j,k;
  double r[3], ji[3];
  double phi, theta, psi;
  double sux, suy, suz;
  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
  double ux, uy, uz, u, uSqr;
  
  AtomStamp* currAtom;
  DirectionalAtom* dAtom;
  double vel[3];
  for(i=0;i<3;i++)vel[i]=0.0;

  for(i=0; i<nAtoms; i++){
    
    currAtom = myStamp->getAtom( i );
    j = atomIndex+i;

    if( currAtom->haveOrientation()){
      
      dAtom = new DirectionalAtom( j, myConfig);
      atomArray[j] = dAtom;
      atomArray[j]->setCoords();

      // Directional Atoms have standard unit vectors which are oriented
      // in space using the three Euler angles.
      
      phi = currAtom->getEulerPhi() * M_PI / 180.0;
      theta = currAtom->getEulerTheta() * M_PI / 180.0;
      psi = currAtom->getEulerPsi()* M_PI / 180.0;

      dAtom->setUnitFrameFromEuler(phi, theta, psi);      
      dAtom->setA( A );

      ji[0] = 0.0;
      ji[1] = 0.0;
      ji[2] = 0.0;
      dAtom->setJ( ji );
      
    }
    else{
      atomArray[j] = new Atom( j, myConfig);
      atomArray[j]->setCoords();
    }
    
    atomArray[j]->setType( currAtom->getType() );
    
    for(k=0; k<3; k++) r[k] = myCoords[(i*3)+k];

    rotMe( r, A );

    for(k=0; k<3; k++) r[k] += pos[k];

    atomArray[j]->setPos( r );
        
    atomArray[j]->setVel( vel );;
  }
}

void MoLocator::calcRefCoords( void ){

  int i,j,k;
  AtomStamp* currAtom;
  double centerX, centerY, centerZ;
  double smallX, smallY, smallZ;
  double bigX, bigY, bigZ;
  double dx, dy, dz;
  double dsqr;

  
  centerX = 0.0;
  centerY = 0.0;
  centerZ = 0.0;

  for(i=0; i<nAtoms; i++){
    
    currAtom = myStamp->getAtom(i);
    if( !currAtom->havePosition() ){
      sprintf( painCave.errMsg,
               "MoLocator error.\n"
               "  Component %s, atom %s does not have a position specified.\n"
               "  This means MoLocator cannot initalize it's position.\n",
               myStamp->getID(),
               currAtom->getType() );
      painCave.isFatal = 1;
      simError();
    }

    
    centerX += currAtom->getPosX();
    centerY += currAtom->getPosY();
    centerZ += currAtom->getPosZ();
  }

  centerX /= nAtoms;
  centerY /= nAtoms;
  centerZ /= nAtoms;
  
  myCoords = new double[nAtoms*3];

  j = 0;
  for(i=0; i<nAtoms; i++){
    
    currAtom = myStamp->getAtom(i);
    j = i*3;
    
    myCoords[j]   = currAtom->getPosX() - centerX;
    myCoords[j+1] = currAtom->getPosY() - centerY;
    myCoords[j+2] = currAtom->getPosZ() - centerZ;
  }
  
  smallX = myCoords[0];
  smallY = myCoords[1];
  smallZ = myCoords[2];

  bigX = myCoords[0];
  bigY = myCoords[1];
  bigZ = myCoords[2];

  j=0;
  for(i=1; i<nAtoms; i++){
    j= i*3;
    
    if( myCoords[j]   < smallX ) smallX = myCoords[j];
    if( myCoords[j+1] < smallY ) smallY = myCoords[j+1];
    if( myCoords[j+2] < smallZ ) smallZ = myCoords[j+2];

    if( myCoords[j]   > bigX ) bigX = myCoords[j];
    if( myCoords[j+1] > bigY ) bigY = myCoords[j+1];
    if( myCoords[j+2] > bigZ ) bigZ = myCoords[j+2];
  }


  dx = bigX - smallX;
  dy = bigY - smallY;
  dz = bigZ - smallZ;

  dsqr = (dx * dx) + (dy * dy) + (dz * dz);
  maxLength = sqrt( dsqr );
}

void MoLocator::rotMe( double r[3], double A[3][3] ){

  double rt[3];
  int i,j;

  for(i=0; i<3; i++) rt[i] = r[i];

  for(i=0; i<3; i++){
    r[i] = 0.0;
    for(j=0; j<3; j++){
      r[i] += A[i][j] * rt[j];
    }
  }
}

void getRandomRot( double rot[3][3] ){

  double theta, phi, psi;
  double cosTheta;

  // select random phi, psi, and cosTheta

  phi = 2.0 * M_PI * drand48();
  psi = 2.0 * M_PI * drand48();
  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1

  theta = acos( cosTheta );

  getEulerRot( theta, phi, psi, rot );
}


void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){

  rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
  rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
  rot[0][2] = sin(theta) * sin(psi);
  
  rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
  rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
  rot[1][2] = sin(theta) * cos(psi);

  rot[2][0] = sin(phi) * sin(theta);
  rot[2][1] = -cos(phi) * sin(theta);
  rot[2][2] = cos(theta);  
}

Revision:	1099
Committed:	Mon Apr 12 20:33:12 2004 UTC (21 years ago) by gezelter
File size:	4704 byte(s)
Log Message:	Changes for RigidBody dynamics
#	User	Rev	Content
1	chuckv	678	#include <iostream>
2
3			#include <cstdlib>
4			#include <cmath>
5
6			#include "simError.h"
7
8			#include "MoLocator.hpp"
9
10
11			MoLocator::MoLocator( MoleculeStamp* theStamp ){
12
13			myStamp = theStamp;
14			nAtoms = myStamp->getNAtoms();
15
16			myCoords = NULL;
17
18			calcRefCoords();
19			}
20
21			MoLocator::~MoLocator(){
22
23			if( myCoords != NULL ) delete[] myCoords;
24			}
25
26			void MoLocator::placeMol( double pos[3], double A[3][3], Atom** atomArray,
27	chuckv	700	int atomIndex, SimState* myConfig ){
28	chuckv	678
29			int i,j,k;
30	gezelter	1099	double r[3], ji[3];
31	gezelter	986	double phi, theta, psi;
32			double sux, suy, suz;
33			double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
34	chuckv	678	double ux, uy, uz, u, uSqr;
35
36			AtomStamp* currAtom;
37			DirectionalAtom* dAtom;
38	chuckv	700	double vel[3];
39			for(i=0;i<3;i++)vel[i]=0.0;
40	chuckv	678
41			for(i=0; i<nAtoms; i++){
42
43			currAtom = myStamp->getAtom( i );
44			j = atomIndex+i;
45
46			if( currAtom->haveOrientation()){
47
48			dAtom = new DirectionalAtom( j, myConfig);
49			atomArray[j] = dAtom;
50	chuckv	700	atomArray[j]->setCoords();
51	gezelter	986
52			// Directional Atoms have standard unit vectors which are oriented
53	gezelter	1099	// in space using the three Euler angles.
54	gezelter	986
55			phi = currAtom->getEulerPhi() * M_PI / 180.0;
56			theta = currAtom->getEulerTheta() * M_PI / 180.0;
57			psi = currAtom->getEulerPsi()* M_PI / 180.0;
58	gezelter	1099
59			dAtom->setUnitFrameFromEuler(phi, theta, psi);
60	chuckv	678	dAtom->setA( A );
61	gezelter	1099
62			ji[0] = 0.0;
63			ji[1] = 0.0;
64			ji[2] = 0.0;
65			dAtom->setJ( ji );
66	chuckv	678
67			}
68			else{
69	gezelter	1099	atomArray[j] = new Atom( j, myConfig);
70	chuckv	700	atomArray[j]->setCoords();
71	chuckv	678	}
72
73			atomArray[j]->setType( currAtom->getType() );
74
75			for(k=0; k<3; k++) r[k] = myCoords[(i*3)+k];
76
77			rotMe( r, A );
78
79			for(k=0; k<3; k++) r[k] += pos[k];
80
81	chuckv	700	atomArray[j]->setPos( r );
82
83			atomArray[j]->setVel( vel );;
84	chuckv	678	}
85			}
86
87			void MoLocator::calcRefCoords( void ){
88
89			int i,j,k;
90			AtomStamp* currAtom;
91			double centerX, centerY, centerZ;
92			double smallX, smallY, smallZ;
93			double bigX, bigY, bigZ;
94			double dx, dy, dz;
95			double dsqr;
96
97
98			centerX = 0.0;
99			centerY = 0.0;
100			centerZ = 0.0;
101
102			for(i=0; i<nAtoms; i++){
103
104			currAtom = myStamp->getAtom(i);
105			if( !currAtom->havePosition() ){
106			sprintf( painCave.errMsg,
107			"MoLocator error.\n"
108			" Component %s, atom %s does not have a position specified.\n"
109			" This means MoLocator cannot initalize it's position.\n",
110			myStamp->getID(),
111			currAtom->getType() );
112			painCave.isFatal = 1;
113			simError();
114			}
115
116
117			centerX += currAtom->getPosX();
118			centerY += currAtom->getPosY();
119			centerZ += currAtom->getPosZ();
120			}
121
122			centerX /= nAtoms;
123			centerY /= nAtoms;
124			centerZ /= nAtoms;
125
126			myCoords = new double[nAtoms*3];
127
128			j = 0;
129			for(i=0; i<nAtoms; i++){
130
131			currAtom = myStamp->getAtom(i);
132			j = i*3;
133
134			myCoords[j] = currAtom->getPosX() - centerX;
135			myCoords[j+1] = currAtom->getPosY() - centerY;
136			myCoords[j+2] = currAtom->getPosZ() - centerZ;
137			}
138
139			smallX = myCoords[0];
140			smallY = myCoords[1];
141			smallZ = myCoords[2];
142
143			bigX = myCoords[0];
144			bigY = myCoords[1];
145			bigZ = myCoords[2];
146
147			j=0;
148			for(i=1; i<nAtoms; i++){
149			j= i*3;
150
151			if( myCoords[j] < smallX ) smallX = myCoords[j];
152			if( myCoords[j+1] < smallY ) smallY = myCoords[j+1];
153			if( myCoords[j+2] < smallZ ) smallZ = myCoords[j+2];
154
155			if( myCoords[j] > bigX ) bigX = myCoords[j];
156			if( myCoords[j+1] > bigY ) bigY = myCoords[j+1];
157			if( myCoords[j+2] > bigZ ) bigZ = myCoords[j+2];
158			}
159
160
161			dx = bigX - smallX;
162			dy = bigY - smallY;
163			dz = bigZ - smallZ;
164
165			dsqr = (dx * dx) + (dy * dy) + (dz * dz);
166			maxLength = sqrt( dsqr );
167			}
168
169			void MoLocator::rotMe( double r[3], double A[3][3] ){
170
171			double rt[3];
172			int i,j;
173
174			for(i=0; i<3; i++) rt[i] = r[i];
175
176			for(i=0; i<3; i++){
177			r[i] = 0.0;
178			for(j=0; j<3; j++){
179			r[i] += A[i][j] * rt[j];
180			}
181			}
182			}
183	mmeineke	821
184			void getRandomRot( double rot[3][3] ){
185
186			double theta, phi, psi;
187			double cosTheta;
188
189			// select random phi, psi, and cosTheta
190
191			phi = 2.0 * M_PI * drand48();
192			psi = 2.0 * M_PI * drand48();
193			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
194
195			theta = acos( cosTheta );
196
197			getEulerRot( theta, phi, psi, rot );
198			}
199
200
201			void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){
202
203			rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
204			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
205			rot[0][2] = sin(theta) * sin(psi);
206
207			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
208			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
209			rot[1][2] = sin(theta) * cos(psi);
210
211			rot[2][0] = sin(phi) * sin(theta);
212			rot[2][1] = -cos(phi) * sin(theta);
213			rot[2][2] = cos(theta);
214			}
215