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 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();
          
      ux = currAtom->getOrntX();
      uy = currAtom->getOrntY();
      uz = currAtom->getOrntZ();
      
      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);  
}


void getUnitRot( double u[3], double rot[3][3] ){

  double theta, phi, psi;

  theta = acos(u[2]);
  phi = atan(u[1] / u[0]);
  psi = 0.0;

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

Revision:	821
Committed:	Fri Oct 24 22:17:45 2003 UTC (21 years, 6 months ago) by mmeineke
File size:	4786 byte(s)
Log Message:	put QuickBass, MoLocator, and latticeBuilder into a Builder Library overhauled latticeBilayer into its own program. Removed sysBuild from the Makefile
#	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			double ux, uy, uz, u, uSqr;
32
33			AtomStamp* currAtom;
34			DirectionalAtom* dAtom;
35	chuckv	700	double vel[3];
36			for(i=0;i<3;i++)vel[i]=0.0;
37	chuckv	678
38			for(i=0; i<nAtoms; i++){
39
40			currAtom = myStamp->getAtom( i );
41			j = atomIndex+i;
42
43			if( currAtom->haveOrientation()){
44
45			dAtom = new DirectionalAtom( j, myConfig);
46			atomArray[j] = dAtom;
47	chuckv	700	atomArray[j]->setCoords();
48	chuckv	678
49			ux = currAtom->getOrntX();
50			uy = currAtom->getOrntY();
51			uz = currAtom->getOrntZ();
52
53			uSqr = (ux * ux) + (uy * uy) + (uz * uz);
54
55			u = sqrt( uSqr );
56			ux = ux / u;
57			uy = uy / u;
58			uz = uz / u;
59
60			dAtom->setSUx( ux );
61			dAtom->setSUy( uy );
62			dAtom->setSUz( uz );
63
64			dAtom->setA( A );
65
66			dAtom->setJx( 0.0 );
67			dAtom->setJy( 0.0 );
68			dAtom->setJz( 0.0 );
69
70			}
71			else{
72			atomArray[j] = new GeneralAtom( j, myConfig);
73	chuckv	700	atomArray[j]->setCoords();
74	chuckv	678	}
75
76			atomArray[j]->setType( currAtom->getType() );
77
78			for(k=0; k<3; k++) r[k] = myCoords[(i*3)+k];
79
80			rotMe( r, A );
81
82			for(k=0; k<3; k++) r[k] += pos[k];
83
84	chuckv	700	atomArray[j]->setPos( r );
85
86			atomArray[j]->setVel( vel );;
87	chuckv	678	}
88			}
89
90			void MoLocator::calcRefCoords( void ){
91
92			int i,j,k;
93			AtomStamp* currAtom;
94			double centerX, centerY, centerZ;
95			double smallX, smallY, smallZ;
96			double bigX, bigY, bigZ;
97			double dx, dy, dz;
98			double dsqr;
99
100
101			centerX = 0.0;
102			centerY = 0.0;
103			centerZ = 0.0;
104
105			for(i=0; i<nAtoms; i++){
106
107			currAtom = myStamp->getAtom(i);
108			if( !currAtom->havePosition() ){
109			sprintf( painCave.errMsg,
110			"MoLocator error.\n"
111			" Component %s, atom %s does not have a position specified.\n"
112			" This means MoLocator cannot initalize it's position.\n",
113			myStamp->getID(),
114			currAtom->getType() );
115			painCave.isFatal = 1;
116			simError();
117			}
118
119
120			centerX += currAtom->getPosX();
121			centerY += currAtom->getPosY();
122			centerZ += currAtom->getPosZ();
123			}
124
125			centerX /= nAtoms;
126			centerY /= nAtoms;
127			centerZ /= nAtoms;
128
129			myCoords = new double[nAtoms*3];
130
131			j = 0;
132			for(i=0; i<nAtoms; i++){
133
134			currAtom = myStamp->getAtom(i);
135			j = i*3;
136
137			myCoords[j] = currAtom->getPosX() - centerX;
138			myCoords[j+1] = currAtom->getPosY() - centerY;
139			myCoords[j+2] = currAtom->getPosZ() - centerZ;
140			}
141
142			smallX = myCoords[0];
143			smallY = myCoords[1];
144			smallZ = myCoords[2];
145
146			bigX = myCoords[0];
147			bigY = myCoords[1];
148			bigZ = myCoords[2];
149
150			j=0;
151			for(i=1; i<nAtoms; i++){
152			j= i*3;
153
154			if( myCoords[j] < smallX ) smallX = myCoords[j];
155			if( myCoords[j+1] < smallY ) smallY = myCoords[j+1];
156			if( myCoords[j+2] < smallZ ) smallZ = myCoords[j+2];
157
158			if( myCoords[j] > bigX ) bigX = myCoords[j];
159			if( myCoords[j+1] > bigY ) bigY = myCoords[j+1];
160			if( myCoords[j+2] > bigZ ) bigZ = myCoords[j+2];
161			}
162
163
164			dx = bigX - smallX;
165			dy = bigY - smallY;
166			dz = bigZ - smallZ;
167
168			dsqr = (dx * dx) + (dy * dy) + (dz * dz);
169			maxLength = sqrt( dsqr );
170			}
171
172			void MoLocator::rotMe( double r[3], double A[3][3] ){
173
174			double rt[3];
175			int i,j;
176
177			for(i=0; i<3; i++) rt[i] = r[i];
178
179			for(i=0; i<3; i++){
180			r[i] = 0.0;
181			for(j=0; j<3; j++){
182			r[i] += A[i][j] * rt[j];
183			}
184			}
185			}
186	mmeineke	821
187			void getRandomRot( double rot[3][3] ){
188
189			double theta, phi, psi;
190			double cosTheta;
191
192			// select random phi, psi, and cosTheta
193
194			phi = 2.0 * M_PI * drand48();
195			psi = 2.0 * M_PI * drand48();
196			cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
197
198			theta = acos( cosTheta );
199
200			getEulerRot( theta, phi, psi, rot );
201			}
202
203
204			void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){
205
206			rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
207			rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
208			rot[0][2] = sin(theta) * sin(psi);
209
210			rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
211			rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
212			rot[1][2] = sin(theta) * cos(psi);
213
214			rot[2][0] = sin(phi) * sin(theta);
215			rot[2][1] = -cos(phi) * sin(theta);
216			rot[2][2] = cos(theta);
217			}
218
219
220			void getUnitRot( double u[3], double rot[3][3] ){
221
222			double theta, phi, psi;
223
224			theta = acos(u[2]);
225			phi = atan(u[1] / u[0]);
226			psi = 0.0;
227
228			getEulerRot( theta, phi, psi, rot );
229			}
230
231