View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
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 real, double precision input of arbitrary dimensions.
42   *
43   * @author Michal J. Schnieders
44   * @see Real
45   * @since 1.0
46   */
47  public class Real3D {
48  
49    private final int nextX, nextY, nextZ;
50    private final int nX, nY, nZ;
51    private final int nX1, nZ2;
52    private final double[] work;
53    private final double[] recip;
54    private final Real fftX;
55    private final Complex fftY, fftZ;
56  
57    /**
58     * Initialize the 3D FFT for complex 3D matrix.
59     *
60     * @param nX X-dimension.
61     * @param nY Y-dimension.
62     * @param nZ Z-dimension.
63     */
64    public Real3D(int nX, int nY, int nZ) {
65      this.nX = nX / 2;
66      this.nY = nY;
67      this.nZ = nZ;
68      nX1 = this.nX + 1;
69      nZ2 = 2 * nZ;
70      nextX = 2;
71      nextY = nX + 2;
72      nextZ = nextY * nY;
73      work = new double[nZ2];
74      recip = new double[nX1 * nY * nZ];
75      fftX = new Real(nX);
76      fftY = new Complex(nY);
77      fftZ = new Complex(nZ);
78    }
79  
80    /**
81     * Determine the index of the real number in the 1D array from the X, Y and Z indices.
82     *
83     * @param i the index along the X-axis.
84     * @param j the index along the Y-axis.
85     * @param k the index along the Z-axis.
86     * @param nX the number of points along the X-axis.
87     * @param nY the number of points along the Y-axis.
88     * @return the index of the real number in the 1D array.
89     */
90    public static int iReal3D(int i, int j, int k, int nX, int nY) {
91      int xSide = nX + 2;
92      int xySlice = xSide * nY;
93      return i + j * xSide + k * xySlice;
94    }
95  
96    /**
97     * convolution
98     *
99     * @param input an array of double.
100    */
101   public void convolution(final double[] input) {
102     for (int z = 0; z < nZ; z++) {
103       for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
104         fftX.fft(input, offset);
105       }
106       for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
107         fftY.fft(input, offset, nextY);
108       }
109     }
110     for (int index = 0, offset = 0, y = 0; y < nY; y++) {
111       for (int x = 0; x < nX1; x++, offset += nextX) {
112         for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
113           work[i] = input[z];
114           work[i + 1] = input[z + 1];
115         }
116         fftZ.fft(work, 0, 2);
117         for (int i = 0; i < nZ2; i += 2) {
118           double r = recip[index++];
119           work[i] *= r;
120           work[i + 1] *= r;
121         }
122         fftZ.ifft(work, 0, 2);
123         for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
124           input[z] = work[i];
125           input[z + 1] = work[i + 1];
126         }
127       }
128     }
129     for (int z = 0; z < nZ; z++) {
130       for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
131         fftY.ifft(input, offset, nextY);
132       }
133       for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
134         fftX.ifft(input, offset);
135       }
136     }
137   }
138 
139   /**
140    * Compute the 3D FFT.
141    *
142    * @param input The input array must be of size (nX + 2) * nY * nZ.
143    */
144   public void fft(final double[] input) {
145     for (int z = 0; z < nZ; z++) {
146       for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
147         fftX.fft(input, offset);
148       }
149       for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
150         fftY.fft(input, offset, nextY);
151       }
152     }
153     for (int x = 0; x < nX1; x++) {
154       for (int offset = x * 2, y = 0; y < nY; y++, offset += nextY) {
155         for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
156           work[i] = input[z];
157           work[i + 1] = input[z + 1];
158         }
159         fftZ.fft(work, 0, 2);
160         for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
161           input[z] = work[i];
162           input[z + 1] = work[i + 1];
163         }
164       }
165     }
166   }
167 
168   /**
169    * Compute the inverse 3D FFT.
170    *
171    * @param input The input array must be of size (nX + 2) * nY * nZ.
172    */
173   public void ifft(final double[] input) {
174     for (int x = 0; x < nX1; x++) {
175       for (int offset = x * 2, y = 0; y < nY; y++, offset += nextY) {
176         for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
177           work[i] = input[z];
178           work[i + 1] = input[z + 1];
179         }
180         fftZ.ifft(work, 0, 2);
181         for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
182           input[z] = work[i];
183           input[z + 1] = work[i + 1];
184         }
185       }
186     }
187     for (int z = 0; z < nZ; z++) {
188       for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
189         fftY.ifft(input, offset, nextY);
190       }
191       for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
192         fftX.ifft(input, offset);
193       }
194     }
195   }
196 
197   /**
198    * Setter for the field <code>recip</code>.
199    *
200    * @param recip an array of double.
201    */
202   public void setRecip(double[] recip) {
203     // Reorder the reciprocal space data into the order it is needed by the convolution routine.
204     for (int index = 0, offset = 0, y = 0; y < nY; y++) {
205       for (int x = 0; x < nX1; x++, offset += 1) {
206         for (int i = 0, z = offset; i < nZ; i++, z += nX1 * nY) {
207           this.recip[index++] = recip[z];
208         }
209       }
210     }
211   }
212 }