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+(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 |
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 |
chrisfen |
1313 |
void SphereHarm::printToShapesFile(char name[200], int index, double tolVal){ |
123 |
chrisfen |
1287 |
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 |
chrisfen |
1313 |
if (fabs(cm) > tolVal * biggest) nfuncs++; |
162 |
|
|
if (fabs(sm) > tolVal * biggest) nfuncs++; |
163 |
chrisfen |
1287 |
} |
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 |
chrisfen |
1313 |
if (fabs(cm) > tolVal * biggest) |
197 |
chrisfen |
1292 |
shapes << l << "\t" << m << "\tcos\t\t" << cm << "\n"; |
198 |
chrisfen |
1313 |
if (fabs(sm) > tolVal * 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( (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 |
|
|
} |