OOPSE/staticProps/AllCorr.cpp

#include <stdlib.h>
#include <math.h>
#include <vector>

#include "simError.h"
#include "AllCorr.hpp"


AllCorr::AllCorr(){

  pairCorrs   = NULL;
  haveFrames  = false;
  pairCorrSet = false;
  haveLength  = false;
  haveNbins = false;
}

AllCorr::~AllCorr(){
  int i;

  if( pairCorrs != NULL ){
    for(i=0; i<nPairs; i++)
      delete pairCorrs[i];
    delete[] pairCorrs;
  }
}

void AllCorr::setMaxLength( double theLength ){
  maxLength = theLength;
  haveLength = true;
}

void AllCorr::setNbins( int theNbins ){
  nBins = theNbins;
  haveNbins = true;
}

void AllCorr::setFrames( SimInfo* theInfoArray, int theNframes, 
                         DumpReader* theReader ){
  int i;
  

  infoArray = theInfoArray;
  nFrames = theNframes;
  reader = theReader;
  
  nAtoms = infoArray[0].n_atoms;

  
  // allocate the first frame's atom arrays

  (infoArray[0].getConfiguration())->createArrays(infoArray[0].n_atoms);
  for (i = 0; i < infoArray[0].n_atoms; i++)
    infoArray[0].atoms[i]->setCoords();

  // read in the first frame's positions

  reader->readFrame( &(infoArray[0]), 0 );

  if( !haveLength ){
    maxLength = infoArray[0].getMaxCutoff();
    haveLength = true;    
  }

  haveFrames = true;
}

void AllCorr::setPairCorrList( vector <PairCorrList> &theList ){
  
  int i;
  char pair1[PCL_NAME_LENGTH];
  char pair2[PCL_NAME_LENGTH];
  
  if( !haveFrames ){
    sprintf( painCave.errMsg,
             "\n"
             "AllCorr Error: attempt to create the pair correlations without"
             " first setting the Frames\n" );
    painCave.isFatal = 1;
    simError();
  }
  
  if( !haveNbins ){
    sprintf( painCave.errMsg,
             "\n"
             "AllCorr Error: attempt to create the pair correlations without"
             " first setting nBins\n" );
    painCave.isFatal = 1;
    simError();
  }


  nPairs = (int)theList.size();
  pairCorrs = new PairCorrType*[nPairs];
  
  for(i=0;i<nPairs;i++){

    theList[i].getAtom1( pair1 );
    theList[i].getAtom2( pair2 );
    
    switch( theList[i].getType() ){

    case gofr:
      pairCorrs[i] = new GofR(pair1, pair2, nAtoms, nBins);
      break;

    case gofrTheta:
      pairCorrs[i] = new GofRtheta(pair1, pair2, nAtoms, nBins);
      break;

    case gofrOmega:
      pairCorrs[i] = new GofRomega(pair1, pair2, nAtoms, nBins);
      break;

    default:
      
      sprintf( painCave.errMsg,
               "AllCorr Error: Unrecognized pair corr type in "
               " setPairCorrList()->theList[%d]\n",
               i );
      painCave.isFatal = 1;
      simError();
      break;
    }
  }

  pairCorrSet = true;
}
  
void AllCorr::initCorrelations( char* theOutPrefix, bool theIsSeparateOut, 
                                bool theHavePairCorrs,
                                bool theHaveStaticCorrs ){
  int i;
  atomName* theNames = new atomName[nAtoms];
  
  havePairCorrs = theHavePairCorrs;
  haveStaticCorrs = theHaveStaticCorrs;

  if( havePairCorrs){
   
    if(!pairCorrSet){
      sprintf( painCave.errMsg,
               "AllCorr error: Attempt to initialize pair correlations "
               " without first setting the pair correlation list.\n" );
      painCave.isFatal = 1;
      simError();
    }
 
    for(i=0;i<nAtoms;i++){
      
      strcpy(theNames[i].name, infoArray[0].atoms[i]->getType() );
    }
    
    for(i=0;i<nPairs;i++){
      pairCorrs[i]->initCorr(maxLength, theNames );
    }
  }
  
  if( haveStaticCorrs ){

    // empty for now
  }

  strcpy(outPrefix, theOutPrefix);
  isSeparateOut = theIsSeparateOut;  
}


void AllCorr::calcCorrelations(){

  SimInfo* currInfo;
  int i,j,k,l;
  double rij[3], rDist, uHatI[3], uHatJ[3];
  double* grndOut = NULL;
  double dSqr, ri[3], rj[3];
  DirectionalAtom* dAtom;
  double frameVolume;
  
 
  for(k=0; k<nFrames; k++){
    
    currInfo = &(infoArray[0]); // only use 1 frame to save memory

    // skip the first frame, as it is already initialized
    if( k != 0 ){
      
      // allocate the frame's atom arrays

      (currInfo->getConfiguration())->createArrays(currInfo->n_atoms);
      for (i = 0; i < currInfo->n_atoms; i++)
        currInfo->atoms[i]->setCoords();
      
      // read in the frame's positions
      
      reader->readFrame( currInfo, k );
    }
    
    frameVolume = currInfo->boxVol;
    this->setFrameVolume( frameVolume );

    if(havePairCorrs){

      // do the pair loop

      for(i=0;i<(nAtoms-1);i++){
        
        if( this->matchI(i) ){
          
          for(j=i+1;j<nAtoms;j++){
            
            if( this->matchJ(j) ){
              
              // calc rdist and rij;
              
              currInfo->atoms[i]->getPos( ri );
                      
              currInfo->atoms[j]->getPos( rj );
              
              for(l=0;l<3;l++)
                rij[l] = rj[l] - ri[l];

              currInfo->wrapVector( rij );

              dSqr=0;
              for(l=0;l<3;l++)
                dSqr += rij[l] * rij[l];
              
              rDist = sqrt( dSqr );
              
              if( currInfo->atoms[i]->isDirectional() ){
                
                dAtom = (DirectionalAtom*)currInfo->atoms[i];
                dAtom->getU( uHatI );
                
                if( currInfo->atoms[j]->isDirectional() ){
                  
                  dAtom = (DirectionalAtom*)currInfo->atoms[j];
                  dAtom->getU( uHatJ );
                  
                  this->pairCorrelate( rij, rDist, uHatI, uHatJ );
                }
                else{
                  
                  this->pairCorrelate( rij, rDist, uHatI, grndOut );
                }
              }
              else if( currInfo->atoms[j]->isDirectional() ){
                
                dAtom = (DirectionalAtom*)currInfo->atoms[j];
                dAtom->getU( uHatJ );
                
                this->pairCorrelate( rij, rDist, grndOut, uHatJ );
              }
              else{
                this->pairCorrelate( rij, rDist, grndOut, grndOut );
              }
            }
          }
        }
      }
    }

    if( haveStaticCorrs ){

      // empty for now 

    }

    // accumulate the averages for this frame
    
    this->accumulateFrame();

    // de-allocate the frame
    
    (currInfo->getConfiguration())->destroyArrays();
  }
                
  this->writeCorrs();  
}

bool AllCorr::matchI(int i){
  
  int k;
  bool result;
  
  result = false;
  
  for(k=0;k<nPairs;k++)
    result = (result || pairCorrs[k]->matchI(i) );

  return result;  
}

bool AllCorr::matchJ(int j){
  
  int k;
  bool result;
  
  result = false;
  
  for(k=0;k<nPairs;k++)
    result = (result || pairCorrs[k]->matchJ(j) );

  return result;  
}

void AllCorr::pairCorrelate( double Rij[3], double dist, 
                             double uHatI[3], double uHatJ[3] ){
  
  int i;
  
  for(i=0;i<nPairs;i++)
    pairCorrs[i]->correlate( Rij, dist, uHatI, uHatJ );
}

void AllCorr::setFrameVolume( double theVolume ){
  
  int i;

  if( havePairCorrs ){

    for(i=0;i<nPairs;i++)
      pairCorrs[i]->setFrameVolume( theVolume );
  }

  if( haveStaticCorrs ){

    // empty for now 
  }
}

void AllCorr::accumulateFrame( void ){
  
  int i;
  
  if( havePairCorrs ){
    
    for(i=0;i<nPairs;i++)
      pairCorrs[i]->accumulateFrame();
  }

  if( haveStaticCorrs ){

    // empty for now
  }
}

void AllCorr::writeCorrs( void ){
  int i;
  
  if( isSeparateOut ){

    if( havePairCorrs ){

      for(i=0;i<nPairs;i++)
        pairCorrs[i]->writeCorr( outPrefix );
    }

    if( haveStaticCorrs ){

      // empty for now
    }

  }
  else{

    // empty for now 
  }

}
Revision:	1043
Committed:	Mon Feb 9 20:38:38 2004 UTC (21 years, 8 months ago) by mmeineke
File size:	6927 byte(s)
Log Message:	started a bug fix on the massive memory overusage by OOPSE
#	User	Rev	Content
1	mmeineke	803	#include <stdlib.h>
2	mmeineke	810	#include <math.h>
3	mmeineke	803	#include <vector>
4
5	mmeineke	810	#include "simError.h"
6	mmeineke	803	#include "AllCorr.hpp"
7
8
9			AllCorr::AllCorr(){
10
11			pairCorrs = NULL;
12			haveFrames = false;
13			pairCorrSet = false;
14			haveLength = false;
15			haveNbins = false;
16			}
17
18			AllCorr::~AllCorr(){
19			int i;
20
21			if( pairCorrs != NULL ){
22			for(i=0; i<nPairs; i++)
23			delete pairCorrs[i];
24			delete[] pairCorrs;
25			}
26			}
27
28			void AllCorr::setMaxLength( double theLength ){
29			maxLength = theLength;
30			haveLength = true;
31			}
32
33	mmeineke	810	void AllCorr::setNbins( int theNbins ){
34	mmeineke	803	nBins = theNbins;
35			haveNbins = true;
36			}
37
38			void AllCorr::setFrames( SimInfo* theInfoArray, int theNframes,
39			DumpReader* theReader ){
40			int i;
41
42
43			infoArray = theInfoArray;
44			nFrames = theNframes;
45			reader = theReader;
46
47			nAtoms = infoArray[0].n_atoms;
48
49
50			// allocate the first frame's atom arrays
51
52			(infoArray[0].getConfiguration())->createArrays(infoArray[0].n_atoms);
53			for (i = 0; i < infoArray[0].n_atoms; i++)
54			infoArray[0].atoms[i]->setCoords();
55
56			// read in the first frame's positions
57
58			reader->readFrame( &(infoArray[0]), 0 );
59
60			if( !haveLength ){
61			maxLength = infoArray[0].getMaxCutoff();
62			haveLength = true;
63			}
64
65			haveFrames = true;
66			}
67
68			void AllCorr::setPairCorrList( vector <PairCorrList> &theList ){
69
70			int i;
71			char pair1[PCL_NAME_LENGTH];
72			char pair2[PCL_NAME_LENGTH];
73
74			if( !haveFrames ){
75			sprintf( painCave.errMsg,
76			"\n"
77			"AllCorr Error: attempt to create the pair correlations without"
78			" first setting the Frames\n" );
79			painCave.isFatal = 1;
80			simError();
81			}
82
83			if( !haveNbins ){
84			sprintf( painCave.errMsg,
85			"\n"
86			"AllCorr Error: attempt to create the pair correlations without"
87			" first setting nBins\n" );
88			painCave.isFatal = 1;
89			simError();
90			}
91
92
93			nPairs = (int)theList.size();
94			pairCorrs = new PairCorrType*[nPairs];
95
96			for(i=0;i<nPairs;i++){
97
98			theList[i].getAtom1( pair1 );
99			theList[i].getAtom2( pair2 );
100
101			switch( theList[i].getType() ){
102
103			case gofr:
104			pairCorrs[i] = new GofR(pair1, pair2, nAtoms, nBins);
105			break;
106
107	mmeineke	885	case gofrTheta:
108			pairCorrs[i] = new GofRtheta(pair1, pair2, nAtoms, nBins);
109			break;
110
111			case gofrOmega:
112			pairCorrs[i] = new GofRomega(pair1, pair2, nAtoms, nBins);
113			break;
114
115	mmeineke	803	default:
116
117			sprintf( painCave.errMsg,
118			"AllCorr Error: Unrecognized pair corr type in "
119			" setPairCorrList()->theList[%d]\n",
120			i );
121			painCave.isFatal = 1;
122			simError();
123			break;
124			}
125			}
126
127			pairCorrSet = true;
128			}
129
130			void AllCorr::initCorrelations( char* theOutPrefix, bool theIsSeparateOut,
131			bool theHavePairCorrs,
132			bool theHaveStaticCorrs ){
133			int i;
134	mmeineke	810	atomName* theNames = new atomName[nAtoms];
135	mmeineke	803
136			havePairCorrs = theHavePairCorrs;
137			haveStaticCorrs = theHaveStaticCorrs;
138
139			if( havePairCorrs){
140
141			if(!pairCorrSet){
142			sprintf( painCave.errMsg,
143			"AllCorr error: Attempt to initialize pair correlations "
144			" without first setting the pair correlation list.\n" );
145			painCave.isFatal = 1;
146			simError();
147			}
148
149			for(i=0;i<nAtoms;i++){
150
151			strcpy(theNames[i].name, infoArray[0].atoms[i]->getType() );
152			}
153
154			for(i=0;i<nPairs;i++){
155	mmeineke	810	pairCorrs[i]->initCorr(maxLength, theNames );
156	mmeineke	803	}
157			}
158
159			if( haveStaticCorrs ){
160
161			// empty for now
162			}
163
164			strcpy(outPrefix, theOutPrefix);
165			isSeparateOut = theIsSeparateOut;
166			}
167
168
169			void AllCorr::calcCorrelations(){
170
171	mmeineke	810	SimInfo* currInfo;
172	mmeineke	803	int i,j,k,l;
173			double rij[3], rDist, uHatI[3], uHatJ[3];
174			double* grndOut = NULL;
175			double dSqr, ri[3], rj[3];
176			DirectionalAtom* dAtom;
177			double frameVolume;
178
179
180			for(k=0; k<nFrames; k++){
181
182	mmeineke	1043	currInfo = &(infoArray[0]); // only use 1 frame to save memory
183	mmeineke	803
184			// skip the first frame, as it is already initialized
185			if( k != 0 ){
186
187			// allocate the frame's atom arrays
188
189			(currInfo->getConfiguration())->createArrays(currInfo->n_atoms);
190			for (i = 0; i < currInfo->n_atoms; i++)
191			currInfo->atoms[i]->setCoords();
192
193			// read in the frame's positions
194
195	mmeineke	1043	reader->readFrame( currInfo, k );
196	mmeineke	803	}
197
198			frameVolume = currInfo->boxVol;
199			this->setFrameVolume( frameVolume );
200
201			if(havePairCorrs){
202
203			// do the pair loop
204
205			for(i=0;i<(nAtoms-1);i++){
206
207			if( this->matchI(i) ){
208
209			for(j=i+1;j<nAtoms;j++){
210
211			if( this->matchJ(j) ){
212
213			// calc rdist and rij;
214
215			currInfo->atoms[i]->getPos( ri );
216	mmeineke	810
217			currInfo->atoms[j]->getPos( rj );
218	mmeineke	803
219			for(l=0;l<3;l++)
220			rij[l] = rj[l] - ri[l];
221	mmeineke	810
222			currInfo->wrapVector( rij );
223
224	mmeineke	803	dSqr=0;
225			for(l=0;l<3;l++)
226	mmeineke	810	dSqr += rij[l] * rij[l];
227	mmeineke	803
228			rDist = sqrt( dSqr );
229
230			if( currInfo->atoms[i]->isDirectional() ){
231
232	mmeineke	810	dAtom = (DirectionalAtom*)currInfo->atoms[i];
233	mmeineke	803	dAtom->getU( uHatI );
234
235			if( currInfo->atoms[j]->isDirectional() ){
236
237	mmeineke	810	dAtom = (DirectionalAtom*)currInfo->atoms[j];
238	mmeineke	803	dAtom->getU( uHatJ );
239
240			this->pairCorrelate( rij, rDist, uHatI, uHatJ );
241			}
242			else{
243
244			this->pairCorrelate( rij, rDist, uHatI, grndOut );
245			}
246			}
247			else if( currInfo->atoms[j]->isDirectional() ){
248
249	mmeineke	810	dAtom = (DirectionalAtom*)currInfo->atoms[j];
250	mmeineke	803	dAtom->getU( uHatJ );
251
252			this->pairCorrelate( rij, rDist, grndOut, uHatJ );
253			}
254			else{
255			this->pairCorrelate( rij, rDist, grndOut, grndOut );
256			}
257			}
258			}
259			}
260			}
261			}
262
263			if( haveStaticCorrs ){
264
265			// empty for now
266
267			}
268
269			// accumulate the averages for this frame
270
271			this->accumulateFrame();
272
273			// de-allocate the frame
274
275			(currInfo->getConfiguration())->destroyArrays();
276			}
277
278			this->writeCorrs();
279			}
280
281			bool AllCorr::matchI(int i){
282
283			int k;
284			bool result;
285
286			result = false;
287
288			for(k=0;k<nPairs;k++)
289			result = (result \|\| pairCorrs[k]->matchI(i) );
290
291			return result;
292			}
293
294			bool AllCorr::matchJ(int j){
295
296			int k;
297			bool result;
298
299			result = false;
300
301			for(k=0;k<nPairs;k++)
302			result = (result \|\| pairCorrs[k]->matchJ(j) );
303
304			return result;
305			}
306
307	mmeineke	810	void AllCorr::pairCorrelate( double Rij[3], double dist,
308			double uHatI[3], double uHatJ[3] ){
309	mmeineke	803
310			int i;
311
312			for(i=0;i<nPairs;i++)
313	mmeineke	810	pairCorrs[i]->correlate( Rij, dist, uHatI, uHatJ );
314	mmeineke	803	}
315
316			void AllCorr::setFrameVolume( double theVolume ){
317
318			int i;
319
320			if( havePairCorrs ){
321
322			for(i=0;i<nPairs;i++)
323			pairCorrs[i]->setFrameVolume( theVolume );
324			}
325
326			if( haveStaticCorrs ){
327
328			// empty for now
329			}
330			}
331
332			void AllCorr::accumulateFrame( void ){
333
334			int i;
335
336			if( havePairCorrs ){
337
338	mmeineke	810	for(i=0;i<nPairs;i++)
339	mmeineke	803	pairCorrs[i]->accumulateFrame();
340			}
341
342			if( haveStaticCorrs ){
343
344			// empty for now
345			}
346			}
347
348			void AllCorr::writeCorrs( void ){
349			int i;
350
351			if( isSeparateOut ){
352
353			if( havePairCorrs ){
354
355			for(i=0;i<nPairs;i++)
356	mmeineke	810	pairCorrs[i]->writeCorr( outPrefix );
357	mmeineke	803	}
358
359			if( haveStaticCorrs ){
360
361			// empty for now
362			}
363
364			}
365			else{
366
367			// empty for now
368			}
369
370			}