OOPSE/utils/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 ){
  
  int i,j,k;
  double r[3];
  double ux, uy, uz, u, uSqr;
  
  AtomStamp* currAtom;
  DirectionalAtom* dAtom;
  
  for(i=0; i<nAtoms; i++){
    
    currAtom = myStamp->getAtom( i );
    j = atomIndex+i;

    if( currAtom->haveOrientation()){
      
      dAtom = new DirectionalAtom( j );
      atomArray[j] = dAtom;
      
      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 );
    }
    
    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]->setX( r[0] );
    atomArray[j]->setY( r[1] );
    atomArray[j]->setZ( r[2] );

    atomArray[j]->set_vx( 0.0 );
    atomArray[j]->set_vy( 0.0 );
    atomArray[j]->set_vz( 0.0 );
  }
}

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];
    }
  }
}
Revision:	504
Committed:	Thu Apr 17 21:54:18 2003 UTC (21 years, 5 months ago) by mmeineke
File size:	3696 byte(s)
Log Message:	fixed up sysBuild to where it should now build our systems
#	User	Rev	Content
1	mmeineke	504	#include <iostream>
2
3	mmeineke	503	#include <cstdlib>
4			#include <cmath>
5
6	mmeineke	504	#include "simError.h"
7
8	mmeineke	501	#include "MoLocator.hpp"
9	mmeineke	503
10
11	mmeineke	504	MoLocator::MoLocator( MoleculeStamp* theStamp ){
12	mmeineke	503
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	mmeineke	504	void MoLocator::placeMol( double pos[3], double A[3][3], Atom** atomArray,
27			int atomIndex ){
28	mmeineke	503
29			int i,j,k;
30			double r[3];
31	mmeineke	504	double ux, uy, uz, u, uSqr;
32	mmeineke	503
33			AtomStamp* currAtom;
34			DirectionalAtom* dAtom;
35
36	mmeineke	504	for(i=0; i<nAtoms; i++){
37
38			currAtom = myStamp->getAtom( i );
39			j = atomIndex+i;
40
41			if( currAtom->haveOrientation()){
42
43			dAtom = new DirectionalAtom( j );
44			atomArray[j] = dAtom;
45
46			ux = currAtom->getOrntX();
47			uy = currAtom->getOrntY();
48			uz = currAtom->getOrntZ();
49
50			uSqr = (ux * ux) + (uy * uy) + (uz * uz);
51
52			u = sqrt( uSqr );
53			ux = ux / u;
54			uy = uy / u;
55			uz = uz / u;
56
57			dAtom->setSUx( ux );
58			dAtom->setSUy( uy );
59			dAtom->setSUz( uz );
60
61			dAtom->setA( A );
62
63			dAtom->setJx( 0.0 );
64			dAtom->setJy( 0.0 );
65			dAtom->setJz( 0.0 );
66
67			}
68			else{
69			atomArray[j] = new GeneralAtom( j );
70			}
71
72			atomArray[j]->setType( currAtom->getType() );
73
74			for(k=0; k<3; k++) r[k] = myCoords[(i*3)+k];
75
76			rotMe( r, A );
77
78			for(k=0; k<3; k++) r[k] += pos[k];
79
80			atomArray[j]->setX( r[0] );
81			atomArray[j]->setY( r[1] );
82			atomArray[j]->setZ( r[2] );
83
84			atomArray[j]->set_vx( 0.0 );
85			atomArray[j]->set_vy( 0.0 );
86			atomArray[j]->set_vz( 0.0 );
87			}
88	mmeineke	503	}
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	mmeineke	504	for(i=0; i<nAtoms; i++){
105	mmeineke	503
106			currAtom = myStamp->getAtom(i);
107	mmeineke	504	if( !currAtom->havePosition() ){
108			sprintf( painCave.errMsg,
109			"MoLocator error.\n"
110			" Component %s, atom %s does not have a position specified.\n"
111			" This means MoLocator cannot initalize it's position.\n",
112			myStamp->getID(),
113			currAtom->getType() );
114			painCave.isFatal = 1;
115			simError();
116			}
117
118	mmeineke	503
119			centerX += currAtom->getPosX();
120			centerY += currAtom->getPosY();
121			centerZ += currAtom->getPosZ();
122			}
123
124			centerX /= nAtoms;
125			centerY /= nAtoms;
126			centerZ /= nAtoms;
127
128			myCoords = new double[nAtoms*3];
129
130			j = 0;
131	mmeineke	504	for(i=0; i<nAtoms; i++){
132	mmeineke	503
133			currAtom = myStamp->getAtom(i);
134			j = i*3;
135
136			myCoords[j] = currAtom->getPosX() - centerX;
137			myCoords[j+1] = currAtom->getPosY() - centerY;
138			myCoords[j+2] = currAtom->getPosZ() - centerZ;
139			}
140
141			smallX = myCoords[0];
142			smallY = myCoords[1];
143			smallZ = myCoords[2];
144
145			bigX = myCoords[0];
146			bigY = myCoords[1];
147			bigZ = myCoords[2];
148
149			j=0;
150			for(i=1; i<nAtoms; i++){
151			j= i*3;
152
153			if( myCoords[j] < smallX ) smallX = myCoords[j];
154			if( myCoords[j+1] < smallY ) smallY = myCoords[j+1];
155			if( myCoords[j+2] < smallZ ) smallZ = myCoords[j+2];
156
157			if( myCoords[j] > bigX ) bigX = myCoords[j];
158			if( myCoords[j+1] > bigY ) bigY = myCoords[j+1];
159			if( myCoords[j+2] > bigZ ) bigZ = myCoords[j+2];
160			}
161
162
163			dx = bigX - smallX;
164			dy = bigY - smallY;
165			dz = bigZ - smallZ;
166
167			dsqr = (dx * dx) + (dy * dy) + (dz * dz);
168			maxLength = sqrt( dsqr );
169			}
170
171			void MoLocator::rotMe( double r[3], double A[3][3] ){
172
173			double rt[3];
174			int i,j;
175
176	mmeineke	504	for(i=0; i<3; i++) rt[i] = r[i];
177	mmeineke	503
178			for(i=0; i<3; i++){
179			r[i] = 0.0;
180			for(j=0; j<3; j++){
181			r[i] += A[i][j] * rt[j];
182			}
183			}
184			}