1 |
#include "SphereHarm.hpp" |
2 |
#include "cospmls.h" |
3 |
#include "makeweights.h" |
4 |
#include "FST_semi_memo.h" |
5 |
|
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+(4*bw*bw), |
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 |
|