00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "QCPartSubDomain.hpp"
00022 #include <algorithm>
00023
00024
00025
00026 const int QCPartSubDomain::DIM;
00027
00028
00029
00033 QCPartSubDomain::QCPartSubDomain (void) :
00034 id(0),
00035 nbAtoms(0),
00036 nbAtomsCore(0),
00037 nbAtomsShell1(0),
00038 nbAtomsShell2(0),
00039 xMin(),
00040 xMax(),
00041 neighborsIdxs(),
00042 atomsInKernel(),
00043 atomIsInKernel(),
00044 atomsInShell1(NULL),
00045 atomsInShell2(NULL),
00046 nbSharedAtoms(),
00047 nbOverlp(0)
00048 {}
00049
00050
00051
00052
00056 QCPartSubDomain::~QCPartSubDomain (void) {
00057
00058 if (atomsInShell1) {
00059 delete [] atomsInShell1;
00060 atomsInShell1 = NULL;
00061 }
00062
00063 if (atomsInShell2) {
00064 delete [] atomsInShell2;
00065 atomsInShell2 = NULL;
00066 }
00067 }
00068
00069
00070
00071
00075 void
00076 QCPartSubDomain::resetData (void) {
00077
00078 nbAtoms = nbAtomsCore = nbAtomsShell1 = nbAtomsShell2 = 0;
00079
00080 atomIsInKernel.clear();
00081 atomsInKernel.clear();
00082
00083 for(unsigned int i = 0 ; i < neighborsIdxs.size(); ++i){
00084 atomsInShell1[i].clear() ;
00085 atomsInShell2[i].clear() ;
00086 }
00087
00088 weight[QC_CORE] = 0;
00089 weight[QC_SHELL1] = 0;
00090 weight[QC_SHELL2] = 0;
00091 }
00092
00093
00094
00095
00099 bool
00100 QCPartSubDomain::initNbOverlapDomain(void) {
00101
00102 bool ret = true;
00103 int size = neighborsIdxs.size() ;
00104
00105 if (size == 0) {
00106 return false ;
00107
00108 } else {
00109
00110 if (atomsInShell1 == NULL) {
00111 atomsInShell1 = new vector<int>[size];
00112 } else {
00113 ret = false;
00114 }
00115
00116 if (atomsInShell2 == NULL) {
00117 atomsInShell2 = new vector<int>[size];
00118 } else {
00119 ret = false;
00120 }
00121 }
00122
00123 return ret;
00124 }
00125
00126
00127
00128
00129
00133 void
00134 QCPartSubDomain::addAtomsInKernel (const QCAtoms& atoms, bool ) {
00135
00136 QCFloat size[DIM], dist[DIM], *xyz;
00137
00138 size[COORDX] = xMax[COORDX] - xMin[COORDX];
00139 size[COORDY] = xMax[COORDY] - xMin[COORDY];
00140 size[COORDZ] = xMax[COORDZ] - xMin[COORDZ];
00141
00142 for (int i = 0; i < atoms.getNbAtoms(); ++i) {
00143
00144 xyz = atoms.getCoordsAt(i);
00145
00146 dist[COORDX] = xyz[COORDX] - xMin[COORDX];
00147 dist[COORDY] = xyz[COORDY] - xMin[COORDY];
00148 dist[COORDZ] = xyz[COORDZ] - xMin[COORDZ];
00149
00150 if ( ( dist[COORDX]>=0 && dist[COORDX]<=size[COORDX] ) &&
00151 ( dist[COORDY]>=0 && dist[COORDY]<=size[COORDY] ) &&
00152 ( dist[COORDZ]>=0 && dist[COORDZ]<=size[COORDZ] ) ) {
00153
00154 atomsInKernel.push_back(i);
00155 ++nbAtoms;
00156 ++nbAtomsCore;
00157 }
00158 }
00159
00160 atomIsInKernel.resize(atomsInKernel.size());
00161 for(int i = 0 ; i < static_cast<int>(atomsInKernel.size()) ; ++i){
00162 atomIsInKernel[i] = 1;
00163 }
00164 }
00168 void
00169 QCPartSubDomain::addAtomsInBuffers (const QCAtoms& atoms, QCPartSubDomain& domain,
00170 int localOverlapSubDomain,
00171 QCFloat buffer1, QCFloat buffer2, bool ) {
00172
00173
00174 QCFloat coeff[DIM], bord[DIM], R, RHalf;
00175 int * lNum;
00176 QCAtomIn atomNSD1, atomNSD2;
00177
00178
00179 lNum = domain.getLocalNum() ;
00180 for(int i = 0; i < DIM ; ++i){
00181 coeff[i] = lNum[i] - localNum[i] ;
00182 if (coeff[i] > 0) {
00183 bord[i] = xMax[i];
00184 } else {
00185 bord[i] = xMin[i];
00186 }
00187 }
00188 RHalf = 0.5 * buffer1;
00189 R = RHalf + buffer2;
00190
00191 vector<int> & index1 = this->getAtomsInKernel() ;
00192 vector<int> & index2 = domain.getAtomsInKernel() ;
00193 vector<int> & atomIsInKernel2 = domain.getAtomIsInKernel() ;
00194 vector<int> indexToRemove ;
00195
00196 atomNSD1.numDomain = id;
00197 atomNSD1.typeZone = QC_SHELL1;
00198 atomNSD1.localNum = -1;
00199 atomNSD2.numDomain = domain.getId();
00200 atomNSD2.typeZone = QC_SHELL1;
00201 atomNSD2.localNum = -1;
00202
00203 int localOverlapSubDomain1, localOverlapSubDomain2 ;
00204 localOverlapSubDomain1 = localOverlapSubDomain ;
00205 localOverlapSubDomain2 = domain.getLocalNumberOfSubdomain(atomNSD1.numDomain);
00206
00207 if (index1.size() == 0 || index2.size() == 0) {
00208 return;
00209 }
00210
00211 int numAtom;
00212 QCFloat * xyz, dist;
00213 for (int at = 0 ; at < static_cast<int>(index2.size()) ; ++at){
00214 numAtom = index2[at];
00215 xyz = atoms.getCoordsAt(numAtom);
00216
00217 dist = FQCAbs((bord[COORDX] - xyz[COORDX])*coeff[COORDX]);
00218 dist = FQCMax(dist, FQCAbs((bord[COORDY] - xyz[COORDY])*coeff[COORDY]));
00219 dist = FQCMax(dist, FQCAbs((bord[COORDZ] - xyz[COORDZ])*coeff[COORDZ]));
00220
00221 if (dist > R)
00222 {continue;}
00223
00224 if (dist < RHalf){
00225 atomIsInKernel2[at] = 0;
00226 indexToRemove.push_back(at);
00227 this->addIndexInOverlap(localOverlapSubDomain1,numAtom);
00228 domain.addIndexInOverlap(localOverlapSubDomain2,numAtom);
00229 }
00230 else {
00231 atomNSD1.typeZone = QC_SHELL2;
00232 this->addIndexInCut(localOverlapSubDomain1,numAtom);
00233 atomNSD1.typeZone = QC_SHELL1;
00234 }
00235 }
00236 sort(indexToRemove.begin(),indexToRemove.end());
00237
00238
00239 vector<int>::iterator startIterator = index2.begin();
00240 index2.erase(startIterator, startIterator + indexToRemove.size());
00241
00242 indexToRemove.clear();
00243
00244 for (int at = 0 ; at < static_cast<int>(index1.size()) ; ++at){
00245 numAtom = index1[at];
00246 xyz = atoms.getCoordsAt(numAtom);
00247
00248 dist = FQCAbs((bord[COORDX] - xyz[COORDX])*coeff[COORDX]);
00249 dist = FQCMax(dist, FQCAbs((bord[COORDY] - xyz[COORDY])*coeff[COORDY]));
00250 dist = FQCMax(dist, FQCAbs((bord[COORDZ] - xyz[COORDZ])*coeff[COORDZ]));
00251
00252 if (dist > R) {
00253 continue; }
00254
00255 if (dist < RHalf){
00256 this->addIndexInOverlap(localOverlapSubDomain1,numAtom);
00257 atomIsInKernel[at] = 0;
00258 indexToRemove.push_back(at);
00259 domain.addIndexInOverlap(localOverlapSubDomain2,numAtom);
00260
00261 }
00262 else {
00263 atomNSD2.typeZone = QC_SHELL2;
00264 domain.addIndexInCut(localOverlapSubDomain2,numAtom);
00265 atomNSD2.typeZone = QC_SHELL1;
00266 }
00267 }
00268
00269 startIterator = index1.end() - indexToRemove.size();
00270 index1.erase(startIterator, index1.end());
00271
00272 indexToRemove.clear() ;
00273 }
00274
00278 void QCPartSubDomain::write (void) {
00279
00280 QC_TRACE_PART("BEGIN QCPartSubDomain::write()");
00281 std::cout << "Domain : "<< id << " ("<< localNum[0]<<", "<< localNum[1]<<", "<< localNum[2]<<") "<<std::endl
00282 << "Numbers : all atoms = "<< nbAtoms<< " core = "<< nbAtomsCore<< " buffer1 = "<< nbAtomsShell1
00283 << " buffer2 = "<< nbAtomsShell2 <<std::endl
00284 << " nbOverlp : " << nbOverlp <<std::endl
00285 << " Size : " << xMin << " "<< xMax <<std::endl;
00286 unsigned int i ;
00287 std::cout << "Atoms in Kernel : " << std::endl << " ";
00288 for( i = 0 ; i < atomsInKernel.size(); ++i) {
00289 std::cout << "("<<i<<" , "<<atomsInKernel[i]<< ") ";
00290 }
00291 std::cout << std::endl ;
00292 int num = atomsInKernel.size() ;
00293 for( i = 0 ; i < neighborsIdxs.size(); ++i) {
00294 std::cout << " Atoms in SHELL1 and in domain " << neighborsIdxs[i] << std::endl << " ";
00295 vector<int> & indexS1 = atomsInShell1[i] ;
00296 for( unsigned int j = 0 ; j < indexS1.size(); ++j) {
00297 std::cout << "("<<num + j<<" , "<<indexS1[j]<< ") ";
00298 }
00299 num += indexS1.size() ;
00300 std::cout << std::endl ;
00301 }
00302 for( i = 0 ; i < neighborsIdxs.size(); ++i) {
00303 std::cout << " Atoms in SHELL2 and in domain " << neighborsIdxs[i] << std::endl << " ";
00304 vector<int> & indexS1 = atomsInShell2[i] ;
00305 for( unsigned int j = 0 ; j < indexS1.size(); ++j) {
00306 std::cout << "("<<num + j<<" , "<<indexS1[j]<< ") ";
00307 }
00308 num += indexS1.size() ;
00309 std::cout << std::endl ;
00310 }
00311
00312 std::cout << "Shared Atoms : " << std::endl << " ";
00313 for( i = 0 ; i < nbSharedAtoms.size(); ++i) {
00314 std::cout << "("<<i<<" , "<<nbSharedAtoms[i]<< ") ";
00315 }
00316 std::cout << std::endl ;
00317
00318 QC_TRACE_PART("END QCPartSubDomain::write ()");
00319 }
00320 void QCPartSubDomain::computeNumbers (void) {
00321 QC_TRACE_PART("BEGIN QCPartSubDomain::computeNumbers ()");
00322
00323 unsigned int i;
00324 nbAtomsCore = nbAtomsShell1 = nbAtomsShell2 = 0;
00325
00326 for( i = 0 ; i < neighborsIdxs.size(); ++i) {
00327 nbAtomsShell1 += atomsInShell1[i].size() ;
00328 nbAtomsShell2 += atomsInShell2[i].size() ;
00329 }
00330 nbAtomsCore = atomsInKernel.size() ;
00331
00332 nbAtoms = nbAtomsCore + nbAtomsShell1 + nbAtomsShell2;
00333 nbOverlp = nbAtomsShell1 + nbAtomsShell2;
00334
00335 #ifndef QC_NO_DEBUG
00336 cout << "Domain " << id
00337 << " (ATOMS): core = " << nbAtomsCore
00338 << " overlap = " << nbAtomsShell1
00339 << " cut = " << nbAtomsShell2
00340 << " total = " << nbAtoms
00341 << " nbOverlp = "<< nbOverlp
00342 << endl;
00343 #endif
00344 QC_TRACE_PART("END QCPartSubDomain::computeNumbers ()");
00345 }