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 shapes(name);
  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){
  ofstream shapes(name, ios::app);
  
  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) > 0.01 * biggest) nfuncs++;
      if (fabs(sm) > 0.01 * 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) > 0.01 * biggest) 
        shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
      if (fabs(sm) > 0.01 * 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:	1309
Committed:	Fri Jun 25 20:10:53 2004 UTC (20 years ago) by chrisfen
File size:	7213 byte(s)
Log Message:	Coefficients are now multiplied by the normalization factor - no longer do we need normalization in the visualizer
#	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			ofstream shapes(name);
115	chrisfen	1298	shapes << "begin ShapeInfo\n";
116			shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n";
117	chrisfen	1287	shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t"
118	chrisfen	1292	<< momInert[1][1] << "\t" << momInert[2][2] << "\n";
119	chrisfen	1298	shapes << "end ShapeInfo\n";
120	chrisfen	1287	}
121
122			void SphereHarm::printToShapesFile(char name[200], int index){
123			ofstream shapes(name, ios::app);
124
125			biggest = 0.0;
126			nfuncs = 0;
127			for ( l = 0 ; l < bw ; l++ ) {
128			for (m = 0; m < l+1; m++) {
129			dummy1 = seanindex(m, l, bw);
130			dummy2 = seanindex(-m, l, bw);
131
132			if (m == 0) {
133	chrisfen	1309	cm = normFactor(l,m)*rcoeffs[dummy1];
134			sm = normFactor(l,m)*icoeffs[dummy1];
135	chrisfen	1287	} else {
136	chrisfen	1309	cm = normFactor(l,m)(pow(-1.0,(double)m)rcoeffs[dummy1]
137			+ rcoeffs[dummy2]);
138			sm = normFactor(l,m)(pow(-1.0,(double)m)icoeffs[dummy1]
139			- icoeffs[dummy2]);
140	chrisfen	1287	}
141
142			if (fabs(cm) > biggest) biggest = fabs(cm);
143			if (fabs(sm) > biggest) biggest = fabs(sm);
144			}
145			}
146			for ( l = 0 ; l < bw ; l++ ) {
147			for (m = 0; m < l+1; m++) {
148			dummy1 = seanindex(m, l, bw);
149			dummy2 = seanindex(-m, l, bw);
150
151			if (m == 0) {
152	chrisfen	1309	cm = normFactor(l,m)*rcoeffs[dummy1];
153			sm = normFactor(l,m)*icoeffs[dummy1];
154	chrisfen	1287	} else {
155	chrisfen	1309	cm = normFactor(l,m)(pow(-1.0,(double)m)rcoeffs[dummy1]
156			+ rcoeffs[dummy2]);
157			sm = normFactor(l,m)(pow(-1.0,(double)m)icoeffs[dummy1]
158			- icoeffs[dummy2]);
159	chrisfen	1287	}
160
161			if (fabs(cm) > 0.01 * biggest) nfuncs++;
162			if (fabs(sm) > 0.01 * biggest) nfuncs++;
163			}
164			}
165
166			switch(index){
167			case 0:{
168			shapes << "\nbegin ContactFunctions\n";
169			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
170			}; break;
171			case 1:{
172			shapes << "\nbegin RangeFunctions\n";
173			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
174			}; break;
175			case 2:{
176			shapes << "\nbegin StrengthFunctions\n";
177			shapes << "#l\tm\tsin or cos\tcoeff (kcal/mol)\n";
178			}; break;
179			}
180
181			for ( l = 0 ; l < bw ; l++ ) {
182			for (m = 0; m < l+1; m++) {
183			dummy1 = seanindex(m, l, bw);
184			dummy2 = seanindex(-m, l, bw);
185
186			if (m == 0) {
187	chrisfen	1309	cm = normFactor(l,m)*rcoeffs[dummy1];
188			sm = normFactor(l,m)*icoeffs[dummy1];
189	chrisfen	1287	} else {
190	chrisfen	1309	cm = normFactor(l,m)(pow(-1.0,(double)m)rcoeffs[dummy1]
191			+ rcoeffs[dummy2]);
192			sm = normFactor(l,m)(pow(-1.0,(double)m)icoeffs[dummy1]
193			- icoeffs[dummy2]);
194	chrisfen	1287	}
195
196			if (fabs(cm) > 0.01 * biggest)
197	chrisfen	1292	shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
198	chrisfen	1287	if (fabs(sm) > 0.01 * biggest)
199	chrisfen	1292	shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n";
200	chrisfen	1287	}
201			}
202			switch(index){
203			case 0:{
204	chrisfen	1292	shapes << "end ContactFunctions\n";
205	chrisfen	1287	}; break;
206			case 1:{
207	chrisfen	1292	shapes << "end RangeFunctions\n";
208	chrisfen	1287	}; break;
209			case 2:{
210	chrisfen	1292	shapes << "end StrengthFunctions\n";
211	chrisfen	1287	}; break;
212			}
213			}
214
215	chrisfen	1309	double SphereHarm::normFactor(int L, int M){
216			// normalization factor:
217			if (L+M > 170){
218			printf("Warning: A coefficient was omitted because l + m > 170.\n"
219			"\tThe double buffer overflows with factorial calculations\n"
220			"\tof 170 and higher. You should consider using a smaller\n"
221			"\tbandwidth if you aren't okay with the loss of the %i, %i\n"
222			"\tspherical harmonic.\n", L, M);
223			return 0.0;
224			}
225			else
226			return sqrt( (2L+1)/(4.0M_PI)*factorialFunc((double)(L-M))
227			/ factorialFunc(double(L+M)) );
228			}
229
230			double SphereHarm::factorialFunc(double n) {
231			if (n < 0.0) return NAN;
232			else {
233			if (n < 2.0) return 1.0;
234			else
235			return n*factorialFunc(n-1.0);
236			}
237			}