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], string particle,
                                      double mass, double momInert[3][3]){
  ofstream shapes(name);
  shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t" 
         << momInert[1][1] << "\t" << momInert[2][2] << "\n\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 = rcoeffs[dummy1];
        sm = icoeffs[dummy1];
      } else {
        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
        sm = 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 = rcoeffs[dummy1];
        sm = icoeffs[dummy1];
      } else {
        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
        sm = 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 = rcoeffs[dummy1];
        sm = icoeffs[dummy1];
      } else {
        cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
        sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
      }
        
      if (fabs(cm) > 0.01 * biggest) 
        shapes << l << "\t" << m << "\tcos\t" << cm << "\n";
      if (fabs(sm) > 0.01 * biggest) 
        shapes << l << "\t" << m << "\tsin\t" << sm << "\n";
    }
  }
  switch(index){
    case 0:{
      shapes << "\nend ContactFunctions\n";
    }; break;
    case 1:{
      shapes << "\nend RangeFunctions\n";
    }; break;
    case 2:{
      shapes << "\nend StrengthFunctions\n";
    }; break;
  }
}

Revision:	1287
Committed:	Wed Jun 23 20:18:48 2004 UTC (20 years ago) by chrisfen
File size:	5979 byte(s)
Log Message:	Major progress towards inclusion of spherical harmonic transform capability - still having some build issues...
#	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			void SphereHarm::printShapesFileStart(char name[200], string particle,
113			double mass, double momInert[3][3]){
114			ofstream shapes(name);
115			shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t"
116			<< momInert[1][1] << "\t" << momInert[2][2] << "\n\n";
117			}
118
119			void SphereHarm::printToShapesFile(char name[200], int index){
120			ofstream shapes(name, ios::app);
121
122			biggest = 0.0;
123			nfuncs = 0;
124			for ( l = 0 ; l < bw ; l++ ) {
125			for (m = 0; m < l+1; m++) {
126			dummy1 = seanindex(m, l, bw);
127			dummy2 = seanindex(-m, l, bw);
128
129			if (m == 0) {
130			cm = rcoeffs[dummy1];
131			sm = icoeffs[dummy1];
132			} else {
133			cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
134			sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
135			}
136
137			if (fabs(cm) > biggest) biggest = fabs(cm);
138			if (fabs(sm) > biggest) biggest = fabs(sm);
139			}
140			}
141			for ( l = 0 ; l < bw ; l++ ) {
142			for (m = 0; m < l+1; m++) {
143			dummy1 = seanindex(m, l, bw);
144			dummy2 = seanindex(-m, l, bw);
145
146			if (m == 0) {
147			cm = rcoeffs[dummy1];
148			sm = icoeffs[dummy1];
149			} else {
150			cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
151			sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
152			}
153
154			if (fabs(cm) > 0.01 * biggest) nfuncs++;
155			if (fabs(sm) > 0.01 * biggest) nfuncs++;
156			}
157			}
158
159			switch(index){
160			case 0:{
161			shapes << "\nbegin ContactFunctions\n";
162			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
163			}; break;
164			case 1:{
165			shapes << "\nbegin RangeFunctions\n";
166			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
167			}; break;
168			case 2:{
169			shapes << "\nbegin StrengthFunctions\n";
170			shapes << "#l\tm\tsin or cos\tcoeff (kcal/mol)\n";
171			}; break;
172			}
173
174			for ( l = 0 ; l < bw ; l++ ) {
175			for (m = 0; m < l+1; m++) {
176			dummy1 = seanindex(m, l, bw);
177			dummy2 = seanindex(-m, l, bw);
178
179			if (m == 0) {
180			cm = rcoeffs[dummy1];
181			sm = icoeffs[dummy1];
182			} else {
183			cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
184			sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
185			}
186
187			if (fabs(cm) > 0.01 * biggest)
188			shapes << l << "\t" << m << "\tcos\t" << cm << "\n";
189			if (fabs(sm) > 0.01 * biggest)
190			shapes << l << "\t" << m << "\tsin\t" << sm << "\n";
191			}
192			}
193			switch(index){
194			case 0:{
195			shapes << "\nend ContactFunctions\n";
196			}; break;
197			case 1:{
198			shapes << "\nend RangeFunctions\n";
199			}; break;
200			case 2:{
201			shapes << "\nend StrengthFunctions\n";
202			}; break;
203			}
204			}
205