OOPSE/props/AllCorr.cpp

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

#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( double 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();
  }


  theList.sort();
  
  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;

    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 atomNames[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]->initCoor(maxLength, theNames );
    }
  }
  
  if( haveStaticCorrs ){

    // empty for now
  }

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


void AllCorr::calcCorrelations(){

  SimInfor* 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[k]);

    // 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( &(infoArray[k]), 0 );
    }
    
    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->wrapVector( ri );
              
              currInfo->atoms[i]->getPos( rj );
              currInfo->wrapVector( rj );             

              for(l=0;l<3;l++)
                rij[l] = rj[l] - ri[l];
              
              dSqr=0;
              for(l=0;l<3;l++)
                dSqr += rj[l] * rj[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[3] Rij, double dist, 
                             double[3] uHatI, double[3] uHatJ ){
  
  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<nPair;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]->writeCorrs( outPrefix );
    }

    if( haveStaticCorrs ){

      // empty for now
    }

  }
  else{

    // empty for now 
  }

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