trunk/SHAPES/SphereHarm.cpp

#include "SphereHarm.hpp"
#include "cospmls.h"
#include "makeweights.h"
#include "FST_semi_memo.h"

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:	1288
Committed:	Wed Jun 23 20:28:39 2004 UTC (20 years ago) by chrisfen
File size:	6052 byte(s)
Log Message:	touch up changes
#	User	Rev	Content
1	chrisfen	1287	#include "SphereHarm.hpp"
2	chrisfen	1288	#include "cospmls.h"
3			#include "makeweights.h"
4			#include "FST_semi_memo.h"
5	chrisfen	1287
6			SphereHarm::SphereHarm( int bandWidth ){
7			bw = bandWidth;
8
9			/* ASSUMING WILL SEMINAIVE ALL ORDERS */
10			cutoff = bw;
11			size = 2*bw;
12
13			/* allocate memory */
14			rdata = (double ) fftw_malloc(sizeof(double) (size * size));
15			idata = (double ) fftw_malloc(sizeof(double) (size * size));
16			rcoeffs = (double ) fftw_malloc(sizeof(double) (bw * bw));
17			icoeffs = (double ) fftw_malloc(sizeof(double) (bw * bw));
18			weights = (double ) fftw_malloc(sizeof(double) 4 * bw);
19			seminaive_naive_tablespace =
20			(double ) fftw_malloc(sizeof(double)
21			(Reduced_Naive_TableSize(bw,cutoff) +
22			Reduced_SpharmonicTableSize(bw,cutoff)));
23			workspace = (double ) fftw_malloc(sizeof(double)
24			((8 * (bw*bw)) +
25			(7 * bw)));
26
27
28			/****
29			At this point, check to see if all the memory has been
30			allocated. If it has not, there's no point in going further.
31			****/
32
33			if ( (rdata == NULL) \|\| (idata == NULL) \|\|
34			(rcoeffs == NULL) \|\| (icoeffs == NULL) \|\|
35			(seminaive_naive_tablespace == NULL) \|\|
36			(workspace == NULL) )
37			{
38			perror("Error in allocating memory");
39			exit( 1 ) ;
40			}
41
42			//precompute the Legendres
43			fprintf(stdout,"Precomputing the Legendres...\n");
44			seminaive_naive_table = SemiNaive_Naive_Pml_Table( bw, cutoff,
45			seminaive_naive_tablespace,
46			workspace );
47
48			//construct fftw plans using the GURU interface
49			/* forward DCT */
50			dctPlan = fftw_plan_r2r_1d( 2*bw, weights, rdata,
51			FFTW_REDFT10, FFTW_ESTIMATE ) ;
52
53			/*
54			fftw "preamble" ;
55			note that this plan places the output in a transposed array
56			*/
57			rank = 1 ;
58			dims[0].n = 2*bw ;
59			dims[0].is = 1 ;
60			dims[0].os = 2*bw ;
61			howmany_rank = 1 ;
62			howmany_dims[0].n = 2*bw ;
63			howmany_dims[0].is = 2*bw ;
64			howmany_dims[0].os = 1 ;
65
66			/* forward fft */
67			fftPlan = fftw_plan_guru_split_dft( rank, dims,
68			howmany_rank, howmany_dims,
69			rdata, idata,
70			workspace, workspace+(4bwbw),
71			FFTW_ESTIMATE );
72
73			//make the weights
74			makeweights( bw, weights );
75			}
76
77			SphereHarm::~SphereHarm(){
78			//free up memory
79			fftw_destroy_plan( fftPlan );
80			fftw_destroy_plan( dctPlan );
81
82			fftw_free(workspace);
83			fftw_free(seminaive_naive_table);
84			fftw_free(seminaive_naive_tablespace);
85			fftw_free(weights);
86			fftw_free(icoeffs);
87			fftw_free(rcoeffs);
88			fftw_free(idata);
89			fftw_free(rdata);
90			}
91
92			void SphereHarm::doTransforms(vector<double> gridData){
93			int i;
94
95			//load the data
96			for (i=0; i<size*size; i++){
97			rdata[i] = gridData[i];
98			//our data is all real, so load the imaginary part with zeros
99			idata[i] = 0.0;
100			}
101
102			//do the forward spherical transform
103			FST_semi_memo(rdata, idata,
104			rcoeffs, icoeffs,
105			bw,
106			seminaive_naive_table,
107			workspace,
108			0,
109			cutoff,
110			&dctPlan,
111			&fftPlan,
112			weights );
113			}
114
115			void SphereHarm::printShapesFileStart(char name[200], string particle,
116			double mass, double momInert[3][3]){
117			ofstream shapes(name);
118			shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t"
119			<< momInert[1][1] << "\t" << momInert[2][2] << "\n\n";
120			}
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			cm = rcoeffs[dummy1];
134			sm = icoeffs[dummy1];
135			} else {
136			cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
137			sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
138			}
139
140			if (fabs(cm) > biggest) biggest = fabs(cm);
141			if (fabs(sm) > biggest) biggest = fabs(sm);
142			}
143			}
144			for ( l = 0 ; l < bw ; l++ ) {
145			for (m = 0; m < l+1; m++) {
146			dummy1 = seanindex(m, l, bw);
147			dummy2 = seanindex(-m, l, bw);
148
149			if (m == 0) {
150			cm = rcoeffs[dummy1];
151			sm = icoeffs[dummy1];
152			} else {
153			cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
154			sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
155			}
156
157			if (fabs(cm) > 0.01 * biggest) nfuncs++;
158			if (fabs(sm) > 0.01 * biggest) nfuncs++;
159			}
160			}
161
162			switch(index){
163			case 0:{
164			shapes << "\nbegin ContactFunctions\n";
165			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
166			}; break;
167			case 1:{
168			shapes << "\nbegin RangeFunctions\n";
169			shapes << "#l\tm\tsin or cos\tcoeff (Ang)\n";
170			}; break;
171			case 2:{
172			shapes << "\nbegin StrengthFunctions\n";
173			shapes << "#l\tm\tsin or cos\tcoeff (kcal/mol)\n";
174			}; break;
175			}
176
177			for ( l = 0 ; l < bw ; l++ ) {
178			for (m = 0; m < l+1; m++) {
179			dummy1 = seanindex(m, l, bw);
180			dummy2 = seanindex(-m, l, bw);
181
182			if (m == 0) {
183			cm = rcoeffs[dummy1];
184			sm = icoeffs[dummy1];
185			} else {
186			cm = pow(-1.0,(double)m)*rcoeffs[dummy1] + rcoeffs[dummy2];
187			sm = pow(-1.0,(double)m)*icoeffs[dummy1] - icoeffs[dummy2];
188			}
189
190			if (fabs(cm) > 0.01 * biggest)
191			shapes << l << "\t" << m << "\tcos\t" << cm << "\n";
192			if (fabs(sm) > 0.01 * biggest)
193			shapes << l << "\t" << m << "\tsin\t" << sm << "\n";
194			}
195			}
196			switch(index){
197			case 0:{
198			shapes << "\nend ContactFunctions\n";
199			}; break;
200			case 1:{
201			shapes << "\nend RangeFunctions\n";
202			}; break;
203			case 2:{
204			shapes << "\nend StrengthFunctions\n";
205			}; break;
206			}
207			}
208