trunk/SHAPES/SphereHarm.cpp

#include "SphereHarm.hpp"

SphereHarm::SphereHarm( int bandWidth ){
  bw = bandWidth;
  
  /*** ASSUMING WILL SEMINAIVE ALL ORDERS ***/
  cutoff = bw;
  size = 2*bw;
  
  /* allocate memory */
  rdata = (double *) fftw_malloc(sizeof(double) * (size * size));
  idata = (double *) fftw_malloc(sizeof(double) * (size * size));
  rcoeffs = (double *) fftw_malloc(sizeof(double) * (bw * bw));
  icoeffs = (double *) fftw_malloc(sizeof(double) * (bw * bw));
  weights = (double *) fftw_malloc(sizeof(double) * 4 * bw);
  seminaive_naive_tablespace =
      (double *) fftw_malloc(sizeof(double) *
                           (Reduced_Naive_TableSize(bw,cutoff) +
                            Reduced_SpharmonicTableSize(bw,cutoff)));
  workspace = (double *) fftw_malloc(sizeof(double) * 
                                ((8 * (bw*bw)) + 
                                 (7 * bw)));
  
  
  /****
    At this point, check to see if all the memory has been
    allocated. If it has not, there's no point in going further.
    ****/
  
  if ( (rdata == NULL) || (idata == NULL) ||
       (rcoeffs == NULL) || (icoeffs == NULL) ||
       (seminaive_naive_tablespace == NULL) ||
       (workspace == NULL) )
  {
    perror("Error in allocating memory");
    exit( 1 ) ;
  }  
  
  //precompute the Legendres
  fprintf(stdout,"Precomputing the Legendres...\n");
  seminaive_naive_table = SemiNaive_Naive_Pml_Table( bw, cutoff,
                                                     seminaive_naive_tablespace,
                                                     workspace );
  
  //construct fftw plans using the GURU interface  
  /* forward DCT */
  dctPlan = fftw_plan_r2r_1d( 2*bw, weights, rdata,
                              FFTW_REDFT10, FFTW_ESTIMATE ) ;
  
  /*
   fftw "preamble" ;
   note that this plan places the output in a transposed array
   */
  rank = 1 ;
  dims[0].n = 2*bw ;
  dims[0].is = 1 ;
  dims[0].os = 2*bw ;
  howmany_rank = 1 ;
  howmany_dims[0].n = 2*bw ;
  howmany_dims[0].is = 2*bw ;
  howmany_dims[0].os = 1 ;
  
  /* forward fft */
  fftPlan = fftw_plan_guru_split_dft( rank, dims,
                                      howmany_rank, howmany_dims,
                                      rdata, idata,
                                      workspace, workspace+(4*bw*bw),
                                      FFTW_ESTIMATE );
  
  //make the weights
  makeweights( bw, weights );
}

SphereHarm::~SphereHarm(){
  //free up memory
  fftw_destroy_plan( fftPlan );
  fftw_destroy_plan( dctPlan );
  
  fftw_free(workspace);
  fftw_free(seminaive_naive_table);
  fftw_free(seminaive_naive_tablespace);
  fftw_free(weights);
  fftw_free(icoeffs);
  fftw_free(rcoeffs);
  fftw_free(idata);
  fftw_free(rdata);
}

void SphereHarm::doTransforms(vector<double> gridData){
  int i;
  
  //load the data
  for (i=0; i<size*size; i++){
    rdata[i] = gridData[i];
    //our data is all real, so load the imaginary part with zeros
    idata[i] = 0.0;
  }
  
  //do the forward spherical transform
  FST_semi_memo(rdata, idata,
                rcoeffs, icoeffs,
                bw,
                seminaive_naive_table,
                workspace,
                0,
                cutoff,
                &dctPlan,
                &fftPlan,
                weights );
}

void SphereHarm::printShapesFileStart(char name[200], char particle[80],
                                      double mass, double momInert[3][3]){
  ofstream myShapes(name);
  printShapesStreamStart(myShapes, particle, mass, momInert);
}

void SphereHarm::printShapesStreamStart(ostream& shapes, char particle[80],
                                      double mass, double momInert[3][3]){
  shapes << "begin ShapeInfo\n";
  shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n";
  shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t" 
         << momInert[1][1] << "\t" << momInert[2][2] << "\n";
  shapes << "end ShapeInfo\n";
}


void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){
  ofstream shapes(name, ios::app);

  printToShapesStream(shapes, index, tolVal);
}

void SphereHarm::printToShapesStream(ostream& shapes, int index, double tolVal) {
  
  biggest = 0.0;
  nfuncs = 0;
  for ( l = 0 ; l < bw ; l++ ) {     
    for (m = 0; m < l+1; m++) {
      dummy1 = seanindex(m, l, bw);
      dummy2 = seanindex(-m, l, bw);
        
      if (m == 0) {
        cm = normFactor(l,m)*rcoeffs[dummy1];
        sm = normFactor(l,m)*icoeffs[dummy1];
      } else {
        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] 
                              + rcoeffs[dummy2]);
        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] 
                              - icoeffs[dummy2]);
      }
        
      if (fabs(cm) > biggest) biggest = fabs(cm);
      if (fabs(sm) > biggest) biggest = fabs(sm);
    }
  }
  for ( l = 0 ; l < bw ; l++ ) {
    for (m = 0; m < l+1; m++) {
      dummy1 = seanindex(m, l, bw);
      dummy2 = seanindex(-m, l, bw);
        
      if (m == 0) {
        cm = normFactor(l,m)*rcoeffs[dummy1];
        sm = normFactor(l,m)*icoeffs[dummy1];
      } else {
        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] 
                              + rcoeffs[dummy2]);
        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] 
                              - icoeffs[dummy2]);
      }
        
      if (fabs(cm) > tolVal * biggest) nfuncs++;
      if (fabs(sm) > tolVal * biggest) nfuncs++;
    }
  }
    
  switch(index){
    case 0:{
      shapes << "\nbegin ContactFunctions\n";
      shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
    }; break;
    case 1:{
      shapes << "\nbegin RangeFunctions\n";
      shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
    }; break;
    case 2:{
      shapes << "\nbegin StrengthFunctions\n";
      shapes << "#l\tm\tsin or cos\tcoeff (kcal/mol)\n";
    }; break;
  }
    
  for ( l = 0 ; l < bw ; l++ ) {
    for (m = 0; m < l+1; m++) {
      dummy1 = seanindex(m, l, bw);
      dummy2 = seanindex(-m, l, bw);
        
      if (m == 0) {
        cm = normFactor(l,m)*rcoeffs[dummy1];
        sm = normFactor(l,m)*icoeffs[dummy1];
      } else {
        cm = normFactor(l,m)*(pow(-1.0,(double)m)*rcoeffs[dummy1] 
                              + rcoeffs[dummy2]);
        sm = normFactor(l,m)*(pow(-1.0,(double)m)*icoeffs[dummy1] 
                              - icoeffs[dummy2]);
      }
        
      if (fabs(cm) > tolVal * biggest) 
        shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
      if (fabs(sm) > tolVal * biggest) 
        shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n";
    }
  }
  switch(index){
    case 0:{
      shapes << "end ContactFunctions\n";
    }; break;
    case 1:{
      shapes << "end RangeFunctions\n";
    }; break;
    case 2:{
      shapes << "end StrengthFunctions\n";
    }; break;
  }
}

double SphereHarm::normFactor(int L, int M){
  // normalization factor:
  if (L+M > 170){
    printf("Warning: A coefficient was omitted because l + m > 170.\n"
           "\tThe double buffer overflows with factorial calculations\n"
           "\tof 170 and higher.  You should consider using a smaller\n"
           "\tbandwidth if you aren't okay with the loss of the %i, %i\n"
           "\tspherical harmonic.\n", L, M);
    return 0.0;
  }
  else
    return sqrt( (2*L+1)/(4.0*M_PI)*factorialFunc((double)(L-M)) 
                 / factorialFunc(double(L+M)) );
}

double SphereHarm::factorialFunc(double n) {
  if (n < 0.0) return NAN;
  else {
    if (n < 2.0) return 1.0;
    else
      return n*factorialFunc(n-1.0);
  }
}
Revision:	1360
Committed:	Tue Jul 20 20:02:15 2004 UTC (21 years, 3 months ago) by gezelter
File size:	7588 byte(s)
Log Message:	Changes for CGI and for gengetopts --unamed-opt="PDBFILE"
#	User	Rev	Content
1	chrisfen	1287	#include "SphereHarm.hpp"
2
3			SphereHarm::SphereHarm( int bandWidth ){
4			bw = bandWidth;
5
6			/* ASSUMING WILL SEMINAIVE ALL ORDERS */
7			cutoff = bw;
8			size = 2*bw;
9
10			/* allocate memory */
11			rdata = (double ) fftw_malloc(sizeof(double) (size * size));
12			idata = (double ) fftw_malloc(sizeof(double) (size * size));
13			rcoeffs = (double ) fftw_malloc(sizeof(double) (bw * bw));
14			icoeffs = (double ) fftw_malloc(sizeof(double) (bw * bw));
15			weights = (double ) fftw_malloc(sizeof(double) 4 * bw);
16			seminaive_naive_tablespace =
17			(double ) fftw_malloc(sizeof(double)
18			(Reduced_Naive_TableSize(bw,cutoff) +
19			Reduced_SpharmonicTableSize(bw,cutoff)));
20			workspace = (double ) fftw_malloc(sizeof(double)
21			((8 * (bw*bw)) +
22			(7 * bw)));
23
24
25			/****
26			At this point, check to see if all the memory has been
27			allocated. If it has not, there's no point in going further.
28			****/
29
30			if ( (rdata == NULL) \|\| (idata == NULL) \|\|
31			(rcoeffs == NULL) \|\| (icoeffs == NULL) \|\|
32			(seminaive_naive_tablespace == NULL) \|\|
33			(workspace == NULL) )
34			{
35			perror("Error in allocating memory");
36			exit( 1 ) ;
37			}
38
39			//precompute the Legendres
40			fprintf(stdout,"Precomputing the Legendres...\n");
41			seminaive_naive_table = SemiNaive_Naive_Pml_Table( bw, cutoff,
42			seminaive_naive_tablespace,
43			workspace );
44
45			//construct fftw plans using the GURU interface
46			/* forward DCT */
47			dctPlan = fftw_plan_r2r_1d( 2*bw, weights, rdata,
48			FFTW_REDFT10, FFTW_ESTIMATE ) ;
49
50			/*
51			fftw "preamble" ;
52			note that this plan places the output in a transposed array
53			*/
54			rank = 1 ;
55			dims[0].n = 2*bw ;
56			dims[0].is = 1 ;
57			dims[0].os = 2*bw ;
58			howmany_rank = 1 ;
59			howmany_dims[0].n = 2*bw ;
60			howmany_dims[0].is = 2*bw ;
61			howmany_dims[0].os = 1 ;
62
63			/* forward fft */
64			fftPlan = fftw_plan_guru_split_dft( rank, dims,
65			howmany_rank, howmany_dims,
66			rdata, idata,
67			workspace, workspace+(4bwbw),
68			FFTW_ESTIMATE );
69
70			//make the weights
71			makeweights( bw, weights );
72			}
73
74			SphereHarm::~SphereHarm(){
75			//free up memory
76			fftw_destroy_plan( fftPlan );
77			fftw_destroy_plan( dctPlan );
78
79			fftw_free(workspace);
80			fftw_free(seminaive_naive_table);
81			fftw_free(seminaive_naive_tablespace);
82			fftw_free(weights);
83			fftw_free(icoeffs);
84			fftw_free(rcoeffs);
85			fftw_free(idata);
86			fftw_free(rdata);
87			}
88
89			void SphereHarm::doTransforms(vector<double> gridData){
90			int i;
91
92			//load the data
93			for (i=0; i<size*size; i++){
94			rdata[i] = gridData[i];
95			//our data is all real, so load the imaginary part with zeros
96			idata[i] = 0.0;
97			}
98
99			//do the forward spherical transform
100			FST_semi_memo(rdata, idata,
101			rcoeffs, icoeffs,
102			bw,
103			seminaive_naive_table,
104			workspace,
105			0,
106			cutoff,
107			&dctPlan,
108			&fftPlan,
109			weights );
110			}
111
112	chrisfen	1292	void SphereHarm::printShapesFileStart(char name[200], char particle[80],
113	chrisfen	1287	double mass, double momInert[3][3]){
114	gezelter	1360	ofstream myShapes(name);
115			printShapesStreamStart(myShapes, particle, mass, momInert);
116			}
117
118			void SphereHarm::printShapesStreamStart(ostream& shapes, char particle[80],
119			double mass, double momInert[3][3]){
120	chrisfen	1298	shapes << "begin ShapeInfo\n";
121			shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n";
122	chrisfen	1287	shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t"
123	chrisfen	1292	<< momInert[1][1] << "\t" << momInert[2][2] << "\n";
124	chrisfen	1298	shapes << "end ShapeInfo\n";
125	chrisfen	1287	}
126
127	gezelter	1360
128
129	chrisfen	1313	void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){
130	chrisfen	1287	ofstream shapes(name, ios::app);
131	gezelter	1360
132			printToShapesStream(shapes, index, tolVal);
133			}
134
135			void SphereHarm::printToShapesStream(ostream& shapes, int index, double tolVal) {
136	chrisfen	1287
137			biggest = 0.0;
138			nfuncs = 0;
139			for ( l = 0 ; l < bw ; l++ ) {
140			for (m = 0; m < l+1; m++) {
141			dummy1 = seanindex(m, l, bw);
142			dummy2 = seanindex(-m, l, bw);
143
144			if (m == 0) {
145	chrisfen	1309	cm = normFactor(l,m)*rcoeffs[dummy1];
146			sm = normFactor(l,m)*icoeffs[dummy1];
147	chrisfen	1287	} else {
148	chrisfen	1309	cm = normFactor(l,m)(pow(-1.0,(double)m)rcoeffs[dummy1]
149			+ rcoeffs[dummy2]);
150			sm = normFactor(l,m)(pow(-1.0,(double)m)icoeffs[dummy1]
151			- icoeffs[dummy2]);
152	chrisfen	1287	}
153
154			if (fabs(cm) > biggest) biggest = fabs(cm);
155			if (fabs(sm) > biggest) biggest = fabs(sm);
156			}
157			}
158			for ( l = 0 ; l < bw ; l++ ) {
159			for (m = 0; m < l+1; m++) {
160			dummy1 = seanindex(m, l, bw);
161			dummy2 = seanindex(-m, l, bw);
162
163			if (m == 0) {
164	chrisfen	1309	cm = normFactor(l,m)*rcoeffs[dummy1];
165			sm = normFactor(l,m)*icoeffs[dummy1];
166	chrisfen	1287	} else {
167	chrisfen	1309	cm = normFactor(l,m)(pow(-1.0,(double)m)rcoeffs[dummy1]
168			+ rcoeffs[dummy2]);
169			sm = normFactor(l,m)(pow(-1.0,(double)m)icoeffs[dummy1]
170			- icoeffs[dummy2]);
171	chrisfen	1287	}
172
173	chrisfen	1313	if (fabs(cm) > tolVal * biggest) nfuncs++;
174			if (fabs(sm) > tolVal * biggest) nfuncs++;
175	chrisfen	1287	}
176			}
177
178			switch(index){
179			case 0:{
180			shapes << "\nbegin ContactFunctions\n";
181			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
182			}; break;
183			case 1:{
184			shapes << "\nbegin RangeFunctions\n";
185			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
186			}; break;
187			case 2:{
188			shapes << "\nbegin StrengthFunctions\n";
189			shapes << "#l\tm\tsin or cos\tcoeff (kcal/mol)\n";
190			}; break;
191			}
192
193			for ( l = 0 ; l < bw ; l++ ) {
194			for (m = 0; m < l+1; m++) {
195			dummy1 = seanindex(m, l, bw);
196			dummy2 = seanindex(-m, l, bw);
197
198			if (m == 0) {
199	chrisfen	1309	cm = normFactor(l,m)*rcoeffs[dummy1];
200			sm = normFactor(l,m)*icoeffs[dummy1];
201	chrisfen	1287	} else {
202	chrisfen	1309	cm = normFactor(l,m)(pow(-1.0,(double)m)rcoeffs[dummy1]
203			+ rcoeffs[dummy2]);
204			sm = normFactor(l,m)(pow(-1.0,(double)m)icoeffs[dummy1]
205			- icoeffs[dummy2]);
206	chrisfen	1287	}
207
208	chrisfen	1313	if (fabs(cm) > tolVal * biggest)
209	chrisfen	1292	shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
210	chrisfen	1313	if (fabs(sm) > tolVal * biggest)
211	chrisfen	1292	shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n";
212	chrisfen	1287	}
213			}
214			switch(index){
215			case 0:{
216	chrisfen	1292	shapes << "end ContactFunctions\n";
217	chrisfen	1287	}; break;
218			case 1:{
219	chrisfen	1292	shapes << "end RangeFunctions\n";
220	chrisfen	1287	}; break;
221			case 2:{
222	chrisfen	1292	shapes << "end StrengthFunctions\n";
223	chrisfen	1287	}; break;
224			}
225			}
226
227	chrisfen	1309	double SphereHarm::normFactor(int L, int M){
228			// normalization factor:
229			if (L+M > 170){
230			printf("Warning: A coefficient was omitted because l + m > 170.\n"
231			"\tThe double buffer overflows with factorial calculations\n"
232			"\tof 170 and higher. You should consider using a smaller\n"
233			"\tbandwidth if you aren't okay with the loss of the %i, %i\n"
234			"\tspherical harmonic.\n", L, M);
235			return 0.0;
236			}
237			else
238			return sqrt( (2L+1)/(4.0M_PI)*factorialFunc((double)(L-M))
239			/ factorialFunc(double(L+M)) );
240			}
241
242			double SphereHarm::factorialFunc(double n) {
243			if (n < 0.0) return NAN;
244			else {
245			if (n < 2.0) return 1.0;
246			else
247			return n*factorialFunc(n-1.0);
248			}
249			}