00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef QC_ATOMS_H_
00022 #define QC_ATOMS_H_
00023
00024
00025 #include <float.h>
00026 #include <string>
00027 #include <vector>
00028
00029 #include "QCPoint.hpp"
00030
00031
00032 template <class TSolver>
00033 class QCMDSystem;
00034
00040
00046
00047
00048
00049 class QCAtoms
00050 {
00051
00052 template <class TSolver> friend class QCMDSystem;
00053 template <class TPSystem> friend class QCLinearPartitioner;
00054
00056 public:
00057
00058
00059
00060
00062 QCAtoms (void):
00063 atoms(NULL), disps(NULL), atomType(NULL), fragTab(NULL), nbAtoms(0),
00064 nbFrags(0),_fragName(NULL),_nbLinks(0),_connect(NULL),_nbAA(0)
00065 {}
00066
00069 QCAtoms (int nb) : nbAtoms(nb) {
00070 atoms = new QCFloat [nb * DIMENSION];
00071 atomType = new int [nb];
00072 fragTab = new int [nb];
00073
00074 _fragName = new std::string[nb];
00075 _connect = new std::vector<int>[nb] ;
00076 nbFrags = nb ;_nbAA = 0 ;
00077 for (int i=0; i<nb; ++i) {
00078 fragTab[i] = i;
00079 }
00080 }
00081
00082
00083
00084
00085
00086 virtual ~QCAtoms (void) {
00087 QC_TRACE_END("BEGIN QCAtoms<TPSolver>::~QCAtoms");
00088 if (atoms) {
00089 delete [] atoms;
00090 }
00091 if (disps) {
00092 delete [] disps;
00093 }
00094 if (atomType) {
00095 delete [] atomType;
00096 }
00097 if (fragTab) {
00098 delete [] fragTab;
00099 }
00100 if (_fragName) {
00101 delete [] _fragName;
00102 }
00103 if (_connect) {
00104 delete [] _connect; _nbLinks=0;
00105 }
00106 QC_TRACE_END("BEGIN QCAtoms<TPSolver>::~QCAtoms");
00107 }
00108
00109
00110
00114 void allocateAtoms (int nb) {
00115 #ifndef QC_NO_DEBUG
00116 assert(atoms == NULL);
00117 assert(atomType == NULL);
00118 assert(fragTab == NULL);
00119 assert(_fragName == NULL);
00120 assert(_connect == NULL);
00121 #endif
00122
00123 atoms = new QCFloat [nb * DIMENSION];
00124 atomType = new int [nb];
00125 fragTab = new int [nb];
00126 _fragName = new std::string [nb];
00127 _connect = new std::vector<int> [nb];
00128 nbAtoms = nb;
00129 nbFrags = nb;
00130 for (int i=0; i<nb; ++i) {
00131 fragTab[i] = i;
00132 }
00133 }
00134
00135
00136
00140 void allocateDisps (void) {
00141 #ifndef QC_NO_DEBUG
00142 assert(disps == NULL);
00143 assert(nbAtoms != 0);
00144 #endif
00145 disps = new QCFloat [nbAtoms*DIMENSION];
00146 resetDisps();
00147 }
00148
00149
00150
00154 void resetDisps (void) {
00155 for (int i=0; i<nbAtoms*DIMENSION; ++i) {
00156 disps[i] = QC_ZERO;
00157 }
00158 }
00159
00160
00161
00162
00166 int getNbAtoms (void) const { return nbAtoms; }
00167
00168
00172 int getNbFrags (void) const { return nbFrags; }
00176 int getNbAA (void) const { return _nbAA; }
00177 void setNbAA (const int n) { _nbAA=n; }
00181 int getNbLinks (void) const { return _nbLinks; }
00182 void setNbLinks (const int nbLink) { _nbLinks = nbLink; }
00183
00187 std::string & fragName (const int i) { return _fragName[i]; }
00188
00189
00193 std::vector<int> & getConnectivity (const int i) { return _connect[i]; }
00194
00195
00202 void getPointAt (int index, QCPoint3D& p) const {
00203 #ifndef QC_NO_DEBUG
00204 assert (index < nbAtoms);
00205 #endif
00206 p[COORDX] = atoms[index * DIMENSION + COORDX];
00207 p[COORDY] = atoms[index * DIMENSION + COORDY];
00208 p[COORDZ] = atoms[index * DIMENSION + COORDZ];
00209 }
00210
00211
00212
00218 void getDispAt (int index, QCPoint3D& p) const {
00219 #ifndef QC_NO_DEBUG
00220 assert (index < nbAtoms);
00221 #endif
00222 p[COORDX] = disps[index * DIMENSION + COORDX];
00223 p[COORDY] = disps[index * DIMENSION + COORDY];
00224 p[COORDZ] = disps[index * DIMENSION + COORDZ];
00225 }
00226
00227
00228
00234 int getType (int index) const {
00235 #ifndef QC_NO_DEBUG
00236 assert (index < nbAtoms);
00237 #endif
00238 int ret = -1;
00239 if (atomType) {
00240 ret = atomType[index];
00241 }
00242 return ret;
00243 }
00244
00245
00246
00250 int getFragment(int index) const {
00251 #ifndef QC_NO_DEBUG
00252 assert (index < nbAtoms);
00253 #endif
00254 int ret = -1;
00255 if (fragTab) {
00256 ret = fragTab[index];
00257 }
00258 return ret;
00259 }
00260
00261
00262
00268 QCFloat * getCoordsAt (int index) const {
00269 #ifndef QC_NO_DEBUG
00270 assert (index < nbAtoms);
00271 #endif
00272 return &atoms[index * DIMENSION + COORDX];
00273 }
00274
00278 void buildGravityCenter(QCPoint3D& G) const {
00279
00280 G[COORDX] = G[COORDY] = G[COORDZ] = 0.0 ;
00281
00282 int index = 0 ;
00283 for(int i = 0 ; i < nbAtoms ; ++i ){
00284 G[COORDX] += atoms[index + COORDX] ;
00285 G[COORDY] += atoms[index + COORDY] ;
00286 G[COORDZ] += atoms[index + COORDZ] ;
00287 index += DIMENSION ;
00288 }
00289
00290 G[COORDX] /= nbAtoms ;
00291 G[COORDY] /= nbAtoms ;
00292 G[COORDZ] /= nbAtoms ;
00293 }
00294
00295
00296
00297
00298
00299
00300 void incrXDisp (int idx, QCFloat val) { disps[idx*DIMENSION+COORDX] += val; }
00301 void incrYDisp (int idx, QCFloat val) { disps[idx*DIMENSION+COORDY] += val; }
00302 void incrZDisp (int idx, QCFloat val) { disps[idx*DIMENSION+COORDZ] += val; }
00303
00304 void divDisps (QCFloat val) {
00305 for (int i=0; i<nbAtoms*DIMENSION; ++i) {
00306 disps[i] /= val;
00307 }
00308 }
00309
00310 void displayDisps (ostream& out);
00311
00312
00313
00317 void swap (const int index1, const int index2) {
00318
00319 int tmp ;
00320 if(atomType != NULL) {
00321 tmp = atomType[index1] ;
00322 atomType[index1] = atomType[index2] ;
00323 atomType[index2] = tmp ;
00324 }
00325 if(fragTab != NULL) {
00326 tmp = fragTab[index1] ;
00327 fragTab[index1] = fragTab[index2] ;
00328 fragTab[index2] = tmp ;
00329 }
00330
00331 QCFloat p ;
00332 p = atoms[index1 * DIMENSION + COORDX];
00333 atoms[index1 * DIMENSION + COORDX] = atoms[index2 * DIMENSION + COORDX];
00334 atoms[index2 * DIMENSION + COORDX] = p;
00335 p = atoms[index1 * DIMENSION + COORDY];
00336 atoms[index1 * DIMENSION + COORDY] = atoms[index2 * DIMENSION + COORDY];
00337 atoms[index2 * DIMENSION + COORDY] = p;
00338 p = atoms[index1 * DIMENSION + COORDZ];
00339 atoms[index1 * DIMENSION + COORDZ] = atoms[index2 * DIMENSION + COORDZ];
00340 atoms[index2 * DIMENSION + COORDZ] = p;
00341
00342 if(disps != NULL) {
00343 p = disps[index1] ;
00344 disps[index1] = disps[index2] ;
00345 disps[index2] = p ;
00346 }
00347 if(_fragName != NULL) {
00348 std::string s(_fragName[index1]) ;
00349 _fragName[index1] = _fragName[index2] ;
00350 _fragName[index2] = s ;
00351 }
00352 if(_connect != NULL) {
00353 std::vector<int> v(_connect[index1]) ;
00354 _connect[index1] = _connect[index2] ;
00355 _connect[index2] = _connect[index1] ;
00356 }
00357
00358 }
00364
00365 void findMinMax (QCPoint3D& atomMin, QCPoint3D& atomMax) const {
00366
00367 atomMax[COORDX] = atomMax[COORDY] = atomMax[COORDZ] = -DBL_MAX;
00368 atomMin[COORDX] = atomMin[COORDY] = atomMin[COORDZ] = DBL_MAX;
00369
00370 int index = 0 ;
00371 for(int i = 0 ; i < nbAtoms ; ++i ){
00372 atomMax[COORDX] = max(atomMax[COORDX], atoms[index + COORDX]);
00373 atomMax[COORDY] = max(atomMax[COORDY], atoms[index + COORDY]);
00374 atomMax[COORDZ] = max(atomMax[COORDZ], atoms[index + COORDZ]);
00375
00376 atomMin[COORDX] = min(atomMin[COORDX], atoms[index + COORDX]);
00377 atomMin[COORDY] = min(atomMin[COORDY], atoms[index + COORDY]);
00378 atomMin[COORDZ] = min(atomMin[COORDZ], atoms[index + COORDZ]);
00379
00380 index += DIMENSION;
00381 }
00382 }
00383
00384
00385
00389 const QCAtoms& getAtoms (void) const { return *this; }
00390
00391
00392
00393
00395 protected:
00396
00397
00403 void setPointAt (int index, const QCFloat * array) {
00404 #ifndef QC_NO_DEBUG
00405 assert (index < nbAtoms);
00406 #endif
00407 atoms[index * DIMENSION + COORDX] = array [COORDX];
00408 atoms[index * DIMENSION + COORDY] = array [COORDY];
00409 atoms[index * DIMENSION + COORDZ] = array [COORDZ];
00410 }
00411
00412
00413
00419 void setPointAt (int index, const QCPoint3D& p) {
00420 #ifndef QC_NO_DEBUG
00421 assert (index < nbAtoms);
00422 #endif
00423 atoms[index * DIMENSION + COORDX] = p [COORDX];
00424 atoms[index * DIMENSION + COORDY] = p [COORDY];
00425 atoms[index * DIMENSION + COORDZ] = p [COORDZ];
00426 }
00427
00428
00429
00435 int setType (int index, int type) {
00436 #ifndef QC_NO_DEBUG
00437 assert (index < nbAtoms);
00438 #endif
00439 return atomType[index] = type;
00440 }
00441
00442
00443
00449 int setFragment (int index, int fragNumber) {
00450 #ifndef QC_NO_DEBUG
00451 assert (index < nbAtoms);
00452 #endif
00453 return fragTab[index] = fragNumber;
00454 }
00455
00459 const QCAtoms& asBase (void) const { return *this; }
00460
00461
00463 private:
00464
00465
00466
00468 public:
00469
00470 static const int DIMENSION = 3;
00471
00472 protected:
00476 QCFloat * QCRestrict atoms;
00477
00478
00482
00483
00484 QCFloat * QCRestrict disps;
00485
00487
00491 int * atomType;
00495 int * QCRestrict fragTab ;
00499 int nbAtoms;
00503 int nbFrags ;
00507 std::string *_fragName ;
00508
00512 int _nbLinks;
00513 std::vector<int> *_connect ;
00517 int _nbAA ;
00518 private:
00519
00520 };
00521
00522
00523
00530 ostream&
00531 operator << (ostream& os, const QCAtoms& atoms);
00532
00533
00534 #endif // QC_ATOMS_H_