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];
  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.  We assume the standard
      // unit vector was originally along the z axis below.
      
      phi = currAtom->getEulerPhi() * M_PI / 180.0;
      theta = currAtom->getEulerTheta() * M_PI / 180.0;
      psi = currAtom->getEulerPsi()* M_PI / 180.0;
      
      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 );
      ux = ux / u;
      uy = uy / u;
      uz = uz / u;
      
      dAtom->setSUx( ux );
      dAtom->setSUy( uy );
      dAtom->setSUz( uz );
      
      dAtom->setA( A );
      
      dAtom->setJx( 0.0 );
      dAtom->setJy( 0.0 );
      dAtom->setJz( 0.0 );
      
    }
    else{
      atomArray[j] = new GeneralAtom( 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:	986
Committed:	Mon Jan 26 22:01:23 2004 UTC (21 years, 3 months ago) by gezelter
File size:	5665 byte(s)
Log Message:	Changes to BASS reader to use Euler angles for orientation instead of unit vectors required changes in MoLocator
#	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			double r[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			// in space using the three Euler angles. We assume the standard
54			// unit vector was originally along the z axis below.
55
56			phi = currAtom->getEulerPhi() * M_PI / 180.0;
57			theta = currAtom->getEulerTheta() * M_PI / 180.0;
58			psi = currAtom->getEulerPsi()* M_PI / 180.0;
59
60			Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
61			Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
62			Axz = sin(theta) * sin(psi);
63
64			Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
65			Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
66			Ayz = sin(theta) * cos(psi);
67
68			Azx = sin(phi) * sin(theta);
69			Azy = -cos(phi) * sin(theta);
70			Azz = cos(theta);
71
72			sux = 0.0;
73			suy = 0.0;
74			suz = 1.0;
75
76			ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
77			uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
78			uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
79	chuckv	678
80			uSqr = (ux * ux) + (uy * uy) + (uz * uz);
81
82			u = sqrt( uSqr );
83			ux = ux / u;
84			uy = uy / u;
85			uz = uz / u;
86
87			dAtom->setSUx( ux );
88			dAtom->setSUy( uy );
89			dAtom->setSUz( uz );
90
91			dAtom->setA( A );
92
93			dAtom->setJx( 0.0 );
94			dAtom->setJy( 0.0 );
95			dAtom->setJz( 0.0 );
96
97			}
98			else{
99			atomArray[j] = new GeneralAtom( j, myConfig);
100	chuckv	700	atomArray[j]->setCoords();
101	chuckv	678	}
102
103			atomArray[j]->setType( currAtom->getType() );
104
105			for(k=0; k<3; k++) r[k] = myCoords[(i*3)+k];
106
107			rotMe( r, A );
108
109			for(k=0; k<3; k++) r[k] += pos[k];
110
111	chuckv	700	atomArray[j]->setPos( r );
112
113			atomArray[j]->setVel( vel );;
114	chuckv	678	}
115			}
116
117			void MoLocator::calcRefCoords( void ){
118
119			int i,j,k;
120			AtomStamp* currAtom;
121			double centerX, centerY, centerZ;
122			double smallX, smallY, smallZ;
123			double bigX, bigY, bigZ;
124			double dx, dy, dz;
125			double dsqr;
126
127
128			centerX = 0.0;
129			centerY = 0.0;
130			centerZ = 0.0;
131
132			for(i=0; i<nAtoms; i++){
133
134			currAtom = myStamp->getAtom(i);
135			if( !currAtom->havePosition() ){
136			sprintf( painCave.errMsg,
137			"MoLocator error.\n"
138			" Component %s, atom %s does not have a position specified.\n"
139			" This means MoLocator cannot initalize it's position.\n",
140			myStamp->getID(),
141			currAtom->getType() );
142			painCave.isFatal = 1;
143			simError();
144			}
145
146
147			centerX += currAtom->getPosX();
148			centerY += currAtom->getPosY();
149			centerZ += currAtom->getPosZ();
150			}
151
152			centerX /= nAtoms;
153			centerY /= nAtoms;
154			centerZ /= nAtoms;
155
156			myCoords = new double[nAtoms*3];
157
158			j = 0;
159			for(i=0; i<nAtoms; i++){
160
161			currAtom = myStamp->getAtom(i);
162			j = i*3;
163
164			myCoords[j] = currAtom->getPosX() - centerX;
165			myCoords[j+1] = currAtom->getPosY() - centerY;
166			myCoords[j+2] = currAtom->getPosZ() - centerZ;
167			}
168
169			smallX = myCoords[0];
170			smallY = myCoords[1];
171			smallZ = myCoords[2];
172
173			bigX = myCoords[0];
174			bigY = myCoords[1];
175			bigZ = myCoords[2];
176
177			j=0;
178			for(i=1; i<nAtoms; i++){
179			j= i*3;
180
181			if( myCoords[j] < smallX ) smallX = myCoords[j];
182			if( myCoords[j+1] < smallY ) smallY = myCoords[j+1];
183			if( myCoords[j+2] < smallZ ) smallZ = myCoords[j+2];
184
185			if( myCoords[j] > bigX ) bigX = myCoords[j];
186			if( myCoords[j+1] > bigY ) bigY = myCoords[j+1];
187			if( myCoords[j+2] > bigZ ) bigZ = myCoords[j+2];
188			}
189
190
191			dx = bigX - smallX;
192			dy = bigY - smallY;
193			dz = bigZ - smallZ;
194
195			dsqr = (dx * dx) + (dy * dy) + (dz * dz);
196			maxLength = sqrt( dsqr );
197			}
198
199			void MoLocator::rotMe( double r[3], double A[3][3] ){
200
201			double rt[3];
202			int i,j;
203
204			for(i=0; i<3; i++) rt[i] = r[i];
205
206			for(i=0; i<3; i++){
207			r[i] = 0.0;
208			for(j=0; j<3; j++){
209			r[i] += A[i][j] * rt[j];
210			}
211			}
212			}
213	mmeineke	821
214			void getRandomRot( double rot[3][3] ){
215
216			double theta, phi, psi;
217			double cosTheta;
218
219			// select random phi, psi, and cosTheta
220
221			phi = 2.0 * M_PI * drand48();
222			psi = 2.0 * M_PI * drand48();
223			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
224
225			theta = acos( cosTheta );
226
227			getEulerRot( theta, phi, psi, rot );
228			}
229
230
231			void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){
232
233			rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
234			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
235			rot[0][2] = sin(theta) * sin(psi);
236
237			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
238			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
239			rot[1][2] = sin(theta) * cos(psi);
240
241			rot[2][0] = sin(phi) * sin(theta);
242			rot[2][1] = -cos(phi) * sin(theta);
243			rot[2][2] = cos(theta);
244			}
245