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;

    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[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->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:	810
Committed:	Fri Oct 17 21:19:07 2003 UTC (20 years, 8 months ago) by mmeineke
File size:	6701 byte(s)
Log Message:	added the staticProps directory to the build list for both configure and configure.in fixed a number of bugs in the staticProps code. gofr is now working.
#	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			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			}
118
119			pairCorrSet = true;
120			}
121
122			void AllCorr::initCorrelations( char* theOutPrefix, bool theIsSeparateOut,
123			bool theHavePairCorrs,
124			bool theHaveStaticCorrs ){
125			int i;
126	mmeineke	810	atomName* theNames = new atomName[nAtoms];
127	mmeineke	803
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	mmeineke	810	pairCorrs[i]->initCorr(maxLength, theNames );
148	mmeineke	803	}
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	mmeineke	810	SimInfo* currInfo;
164	mmeineke	803	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	mmeineke	810
209			currInfo->atoms[j]->getPos( rj );
210	mmeineke	803
211			for(l=0;l<3;l++)
212			rij[l] = rj[l] - ri[l];
213	mmeineke	810
214			currInfo->wrapVector( rij );
215
216	mmeineke	803	dSqr=0;
217			for(l=0;l<3;l++)
218	mmeineke	810	dSqr += rij[l] * rij[l];
219	mmeineke	803
220			rDist = sqrt( dSqr );
221
222			if( currInfo->atoms[i]->isDirectional() ){
223
224	mmeineke	810	dAtom = (DirectionalAtom*)currInfo->atoms[i];
225	mmeineke	803	dAtom->getU( uHatI );
226
227			if( currInfo->atoms[j]->isDirectional() ){
228
229	mmeineke	810	dAtom = (DirectionalAtom*)currInfo->atoms[j];
230	mmeineke	803	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	mmeineke	810	dAtom = (DirectionalAtom*)currInfo->atoms[j];
242	mmeineke	803	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	mmeineke	810	void AllCorr::pairCorrelate( double Rij[3], double dist,
300			double uHatI[3], double uHatJ[3] ){
301	mmeineke	803
302			int i;
303
304			for(i=0;i<nPairs;i++)
305	mmeineke	810	pairCorrs[i]->correlate( Rij, dist, uHatI, uHatJ );
306	mmeineke	803	}
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	mmeineke	810	for(i=0;i<nPairs;i++)
331	mmeineke	803	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	mmeineke	810	pairCorrs[i]->writeCorr( outPrefix );
349	mmeineke	803	}
350
351			if( haveStaticCorrs ){
352
353			// empty for now
354			}
355
356			}
357			else{
358
359			// empty for now
360			}
361
362			}