1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024.
6 //
7 // This file is part of Force Field X.
8 //
9 // Force Field X is free software; you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License version 3 as published by
11 // the Free Software Foundation.
12 //
13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16 // details.
17 //
18 // You should have received a copy of the GNU General Public License along with
19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20 // Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // Linking this library statically or dynamically with other modules is making a
23 // combined work based on this library. Thus, the terms and conditions of the
24 // GNU General Public License cover the whole combination.
25 //
26 // As a special exception, the copyright holders of this library give you
27 // permission to link this library with independent modules to produce an
28 // executable, regardless of the license terms of these independent modules, and
29 // to copy and distribute the resulting executable under terms of your choice,
30 // provided that you also meet, for each linked independent module, the terms
31 // and conditions of the license of that module. An independent module is a
32 // module which is not derived from or based on this library. If you modify this
33 // library, you may extend this exception to your version of the library, but
34 // you are not obligated to do so. If you do not wish to do so, delete this
35 // exception statement from your version.
36 //
37 // ******************************************************************************
38 package ffx.numerics.fft;
39
40 /**
41 * Compute the 3D FFT of complex, double precision input of arbitrary dimensions via 1D Mixed Radix
42 * FFTs.
43 *
44 * <p>The location of the input point [i, j, k] within the input array must be: <br>
45 * double real = input[x*nextX + y*nextY + z*nextZ] <br> double imag = input[x*nextX + y*nextY +
46 * z*nextZ + 1] <br> where <br> int nextX = 2 <br> int nextY = 2*nX <br> int nextZ = 2*nX*nY<br>
47 *
48 * <p>
49 *
50 * @author Michal J. Schnieders
51 * @see Complex
52 * @since 1.0
53 */
54 public class Complex3D {
55
56 private final int nX, nY, nZ;
57 private final int nZ2;
58 private final int nextX, nextY, nextZ;
59 private final double[] work;
60 private final Complex fftX, fftY, fftZ;
61 private final double[] recip;
62
63 /**
64 * Initialize the 3D FFT for complex 3D matrix.
65 *
66 * @param nX X-dimension.
67 * @param nY Y-dimension.
68 * @param nZ Z-dimension.
69 */
70 public Complex3D(int nX, int nY, int nZ) {
71 this.nX = nX;
72 this.nY = nY;
73 this.nZ = nZ;
74 nZ2 = 2 * nZ;
75 nextX = 2;
76 nextY = 2 * nX;
77 nextZ = nextY * nY;
78
79 work = new double[nZ2];
80 recip = new double[nX * nY * nZ];
81
82 fftX = new Complex(nX);
83 fftY = new Complex(nY);
84 fftZ = new Complex(nZ);
85 }
86
87 /**
88 * Determine the index of the complex number in the 1D array from the X, Y and Z indices.
89 *
90 * @param i the index along the X-axis.
91 * @param j the index along the Y-axis.
92 * @param k the index along the Z-axis.
93 * @param nX the number of points along the X-axis.
94 * @param nY the number of points along the Y-axis.
95 * @return the index of the complex number in the 1D array.
96 */
97 public static int iComplex3D(int i, int j, int k, int nX, int nY) {
98 return 2 * (i + nX * (j + nY * k));
99 }
100
101 /**
102 * Perform a convolution.
103 *
104 * @param input The input array.
105 */
106 public void convolution(final double[] input) {
107 for (int z = 0; z < nZ; z++) {
108 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
109 fftX.fft(input, offset, nextX);
110 }
111 for (int offset = z * nextZ, x = 0; x < nX; x++, offset += nextX) {
112 fftY.fft(input, offset, nextY);
113 }
114 }
115 for (int offset = 0, index = 0, y = 0; y < nY; y++) {
116 for (int x = 0; x < nX; x++, offset += nextX) {
117 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
118 work[i] = input[z];
119 work[i + 1] = input[z + 1];
120 }
121 fftZ.fft(work, 0, 2);
122 for (int i = 0; i < nZ2; i += 2) {
123 double r = recip[index++];
124 work[i] *= r;
125 work[i + 1] *= r;
126 }
127 fftZ.ifft(work, 0, 2);
128 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
129 input[z] = work[i];
130 input[z + 1] = work[i + 1];
131 }
132 }
133 }
134 for (int z = 0; z < nZ; z++) {
135 for (int offset = z * nextZ, x = 0; x < nX; x++, offset += nextX) {
136 fftY.ifft(input, offset, nextY);
137 }
138 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
139 fftX.ifft(input, offset, nextX);
140 }
141 }
142 }
143
144 /**
145 * Compute the 3D FFT.
146 *
147 * @param input The input array must be of size 2 * nX * nY * nZ.
148 */
149 public void fft(final double[] input) {
150 for (int z = 0; z < nZ; z++) {
151 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
152 fftX.fft(input, offset, nextX);
153 }
154 for (int offset = z * nextZ, x = 0; x < nX; x++, offset += nextX) {
155 fftY.fft(input, offset, nextY);
156 }
157 }
158 for (int x = 0, offset = 0; x < nX; x++) {
159 for (int y = 0; y < nY; y++, offset += nextX) {
160 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
161 work[i] = input[z];
162 work[i + 1] = input[z + 1];
163 }
164 fftZ.fft(work, 0, 2);
165 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
166 input[z] = work[i];
167 input[z + 1] = work[i + 1];
168 }
169 }
170 }
171 }
172
173 /**
174 * Compute the inverse 3D FFT.
175 *
176 * @param input The input array must be of size 2 * nX * nY * nZ.
177 */
178 public void ifft(final double[] input) {
179 for (int offset = 0, x = 0; x < nX; x++) {
180 for (int y = 0; y < nY; y++, offset += nextX) {
181 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
182 work[i] = input[z];
183 work[i + 1] = input[z + 1];
184 }
185 fftZ.ifft(work, 0, 2);
186 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
187 input[z] = work[i];
188 input[z + 1] = work[i + 1];
189 }
190 }
191 }
192 for (int z = 0; z < nZ; z++) {
193 for (int offset = z * nextZ, x = 0; x < nX; x++, offset += nextX) {
194 fftY.ifft(input, offset, nextY);
195 }
196 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
197 fftX.ifft(input, offset, nextX);
198 }
199 }
200 }
201
202 /**
203 * Setter for the field <code>recip</code>.
204 *
205 * @param recip an array of double.
206 */
207 public void setRecip(double[] recip) {
208 // Reorder the reciprocal space data into the order it is needed by the convolution routine.
209 for (int offset = 0, index = 0, y = 0; y < nY; y++) {
210 for (int x = 0; x < nX; x++, offset += 1) {
211 for (int i = 0, z = offset; i < nZ2; i += 2, z += nX * nY) {
212 this.recip[index++] = recip[z];
213 }
214 }
215 }
216 }
217 }