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();

  myConfig=NULL;
  
  myCoords = NULL;
  
  calcRefCoords();
}

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

void MoLocator::setConfig(SimState * theConfig){

        myConfig = theConfig;
        haveConfig = true;
}

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;


if( !haveConfig ){
        sprintf(painCave.errMsg,
                "attempt to placeMol without setting the SimState in the MoLocator.\n",);
        painCave.isFatal = 1;
        simError(void);
}
  
  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(void);
          
      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(void);
    }
    
    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:	678
Committed:	Mon Aug 11 22:12:31 2003 UTC (20 years, 11 months ago) by chuckv
File size:	4059 byte(s)
Log Message:	Arranged sysbuilder into a subdirectory. Fixed some of sysbuilder to work with new atom allocation in libmdtools.
#	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			myConfig=NULL;
17
18			myCoords = NULL;
19
20			calcRefCoords();
21			}
22
23			MoLocator::~MoLocator(){
24
25			if( myCoords != NULL ) delete[] myCoords;
26			}
27
28			void MoLocator::setConfig(SimState * theConfig){
29
30			myConfig = theConfig;
31			haveConfig = true;
32			}
33
34			void MoLocator::placeMol( double pos[3], double A[3][3], Atom** atomArray,
35			int atomIndex ){
36
37			int i,j,k;
38			double r[3];
39			double ux, uy, uz, u, uSqr;
40
41			AtomStamp* currAtom;
42			DirectionalAtom* dAtom;
43
44
45			if( !haveConfig ){
46			sprintf(painCave.errMsg,
47			"attempt to placeMol without setting the SimState in the MoLocator.\n",);
48			painCave.isFatal = 1;
49			simError(void);
50			}
51
52			for(i=0; i<nAtoms; i++){
53
54			currAtom = myStamp->getAtom( i );
55			j = atomIndex+i;
56
57			if( currAtom->haveOrientation()){
58
59			dAtom = new DirectionalAtom( j, myConfig);
60			atomArray[j] = dAtom;
61			atomArray[j]->setCoords(void);
62
63			ux = currAtom->getOrntX();
64			uy = currAtom->getOrntY();
65			uz = currAtom->getOrntZ();
66
67			uSqr = (ux * ux) + (uy * uy) + (uz * uz);
68
69			u = sqrt( uSqr );
70			ux = ux / u;
71			uy = uy / u;
72			uz = uz / u;
73
74			dAtom->setSUx( ux );
75			dAtom->setSUy( uy );
76			dAtom->setSUz( uz );
77
78			dAtom->setA( A );
79
80			dAtom->setJx( 0.0 );
81			dAtom->setJy( 0.0 );
82			dAtom->setJz( 0.0 );
83
84			}
85			else{
86			atomArray[j] = new GeneralAtom( j, myConfig);
87			atomArray[j]->setCoords(void);
88			}
89
90			atomArray[j]->setType( currAtom->getType() );
91
92			for(k=0; k<3; k++) r[k] = myCoords[(i*3)+k];
93
94			rotMe( r, A );
95
96			for(k=0; k<3; k++) r[k] += pos[k];
97
98			atomArray[j]->setX( r[0] );
99			atomArray[j]->setY( r[1] );
100			atomArray[j]->setZ( r[2] );
101
102			atomArray[j]->set_vx( 0.0 );
103			atomArray[j]->set_vy( 0.0 );
104			atomArray[j]->set_vz( 0.0 );
105			}
106			}
107
108			void MoLocator::calcRefCoords( void ){
109
110			int i,j,k;
111			AtomStamp* currAtom;
112			double centerX, centerY, centerZ;
113			double smallX, smallY, smallZ;
114			double bigX, bigY, bigZ;
115			double dx, dy, dz;
116			double dsqr;
117
118
119			centerX = 0.0;
120			centerY = 0.0;
121			centerZ = 0.0;
122
123			for(i=0; i<nAtoms; i++){
124
125			currAtom = myStamp->getAtom(i);
126			if( !currAtom->havePosition() ){
127			sprintf( painCave.errMsg,
128			"MoLocator error.\n"
129			" Component %s, atom %s does not have a position specified.\n"
130			" This means MoLocator cannot initalize it's position.\n",
131			myStamp->getID(),
132			currAtom->getType() );
133			painCave.isFatal = 1;
134			simError();
135			}
136
137
138			centerX += currAtom->getPosX();
139			centerY += currAtom->getPosY();
140			centerZ += currAtom->getPosZ();
141			}
142
143			centerX /= nAtoms;
144			centerY /= nAtoms;
145			centerZ /= nAtoms;
146
147			myCoords = new double[nAtoms*3];
148
149			j = 0;
150			for(i=0; i<nAtoms; i++){
151
152			currAtom = myStamp->getAtom(i);
153			j = i*3;
154
155			myCoords[j] = currAtom->getPosX() - centerX;
156			myCoords[j+1] = currAtom->getPosY() - centerY;
157			myCoords[j+2] = currAtom->getPosZ() - centerZ;
158			}
159
160			smallX = myCoords[0];
161			smallY = myCoords[1];
162			smallZ = myCoords[2];
163
164			bigX = myCoords[0];
165			bigY = myCoords[1];
166			bigZ = myCoords[2];
167
168			j=0;
169			for(i=1; i<nAtoms; i++){
170			j= i*3;
171
172			if( myCoords[j] < smallX ) smallX = myCoords[j];
173			if( myCoords[j+1] < smallY ) smallY = myCoords[j+1];
174			if( myCoords[j+2] < smallZ ) smallZ = myCoords[j+2];
175
176			if( myCoords[j] > bigX ) bigX = myCoords[j];
177			if( myCoords[j+1] > bigY ) bigY = myCoords[j+1];
178			if( myCoords[j+2] > bigZ ) bigZ = myCoords[j+2];
179			}
180
181
182			dx = bigX - smallX;
183			dy = bigY - smallY;
184			dz = bigZ - smallZ;
185
186			dsqr = (dx * dx) + (dy * dy) + (dz * dz);
187			maxLength = sqrt( dsqr );
188			}
189
190			void MoLocator::rotMe( double r[3], double A[3][3] ){
191
192			double rt[3];
193			int i,j;
194
195			for(i=0; i<3; i++) rt[i] = r[i];
196
197			for(i=0; i<3; i++){
198			r[i] = 0.0;
199			for(j=0; j<3; j++){
200			r[i] += A[i][j] * rt[j];
201			}
202			}
203			}