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+(4*bw*bw), |
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 shapes(name); |
115 |
shapes << "begin ShapeInfo\n"; |
116 |
shapes << "#name\t\tmass\tI_xx\tI_yy\tI_zz\n"; |
117 |
shapes << particle << "\t" << mass << "\t" << momInert[0][0] << "\t" |
118 |
<< momInert[1][1] << "\t" << momInert[2][2] << "\n"; |
119 |
shapes << "end ShapeInfo\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 = normFactor(l,m)*rcoeffs[dummy1]; |
134 |
sm = normFactor(l,m)*icoeffs[dummy1]; |
135 |
} else { |
136 |
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 |
} |
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 |
cm = normFactor(l,m)*rcoeffs[dummy1]; |
153 |
sm = normFactor(l,m)*icoeffs[dummy1]; |
154 |
} else { |
155 |
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 |
} |
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 |
cm = normFactor(l,m)*rcoeffs[dummy1]; |
188 |
sm = normFactor(l,m)*icoeffs[dummy1]; |
189 |
} else { |
190 |
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 |
} |
195 |
|
196 |
if (fabs(cm) > 0.01 * biggest) |
197 |
shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n"; |
198 |
if (fabs(sm) > 0.01 * biggest) |
199 |
shapes << l << "\t" << m << "\tsin\t\t" << sm << "\n"; |
200 |
} |
201 |
} |
202 |
switch(index){ |
203 |
case 0:{ |
204 |
shapes << "end ContactFunctions\n"; |
205 |
}; break; |
206 |
case 1:{ |
207 |
shapes << "end RangeFunctions\n"; |
208 |
}; break; |
209 |
case 2:{ |
210 |
shapes << "end StrengthFunctions\n"; |
211 |
}; break; |
212 |
} |
213 |
} |
214 |
|
215 |
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( (2*L+1)/(4.0*M_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 |
} |