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 (19 years, 11 months ago) by gezelter
File size:	7588 byte(s)
Log Message:	Changes for CGI and for gengetopts --unamed-opt="PDBFILE"
#	Content
1	#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	void SphereHarm::printShapesFileStart(char name[200], char particle[80],
113	double mass, double momInert[3][3]){
114	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	shapes << "begin ShapeInfo\n";
121	shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n";
122	shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t"
123	<< momInert[1][1] << "\t" << momInert[2][2] << "\n";
124	shapes << "end ShapeInfo\n";
125	}
126
127
128
129	void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){
130	ofstream shapes(name, ios::app);
131
132	printToShapesStream(shapes, index, tolVal);
133	}
134
135	void SphereHarm::printToShapesStream(ostream& shapes, int index, double tolVal) {
136
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	cm = normFactor(l,m)*rcoeffs[dummy1];
146	sm = normFactor(l,m)*icoeffs[dummy1];
147	} else {
148	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	}
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	cm = normFactor(l,m)*rcoeffs[dummy1];
165	sm = normFactor(l,m)*icoeffs[dummy1];
166	} else {
167	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	}
172
173	if (fabs(cm) > tolVal * biggest) nfuncs++;
174	if (fabs(sm) > tolVal * biggest) nfuncs++;
175	}
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	cm = normFactor(l,m)*rcoeffs[dummy1];
200	sm = normFactor(l,m)*icoeffs[dummy1];
201	} else {
202	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	}
207
208	if (fabs(cm) > tolVal * biggest)
209	shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n";
210	if (fabs(sm) > tolVal * biggest)
211	shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n";
212	}
213	}
214	switch(index){
215	case 0:{
216	shapes << "end ContactFunctions\n";
217	}; break;
218	case 1:{
219	shapes << "end RangeFunctions\n";
220	}; break;
221	case 2:{
222	shapes << "end StrengthFunctions\n";
223	}; break;
224	}
225	}
226
227	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	}