Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include <iostream>
23 : #include <complex>
24 : #include "gridtools/ActionWithGrid.h"
25 : #include "core/ActionRegister.h"
26 : #ifdef __PLUMED_HAS_FFTW
27 : #include <fftw3.h> // FFTW interface
28 : #endif
29 :
30 : namespace PLMD {
31 : namespace fourier {
32 :
33 : //+PLUMEDOC GRIDANALYSIS FOURIER_TRANSFORM
34 : /*
35 : Compute the Discrete Fourier Transform (DFT) by means of FFTW of data stored on a 2D grid.
36 :
37 : This action takes a function on a two-dimensional grid as input and computes a fourier transform upon the input function using [FFTW](https://www.fftw.org).
38 : Currently, this actions performs a complex fourier transition even if the input data is real. The functionality here was developed and used in the paper cited
39 : below. The following input was used in that paper:
40 :
41 : ```plumed
42 : UNITS NATURAL
43 :
44 : # These two commands calculate one symmetry function for each atom. These
45 : # symmetry functions tell us whether the environment around each atom resembles
46 : # the environment in the solid or the environment in the liquid.
47 : fcc: FCCUBIC SPECIES=1-20736 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
48 : smapfcc: MORE_THAN ARG=fcc SWITCH={SMAP R_0=0.5 A=8 B=8}
49 : smapfcc_grp: GROUP ATOMS=1-20736
50 : # This calculates the position of the center of the solid region of the simulation box. What we are computing here a weighted average position
51 : # the weights are the order parameters computed using the two commands above.
52 : center: CENTER PHASES ATOMS=fcc WEIGHTS=smapfcc
53 : # This calculates the phase field that tells us whether the structure in each part of the simulation box is solid-like or liquid like.
54 : dens: MULTICOLVARDENS DATA=smapfcc ORIGIN=center DIR=xyz NBINS=50,80,80 BANDWIDTH=1.0,1.0,1.0 GRID_MIN=0.0,auto,auto GRID_MAX=20.0,auto,auto STRIDE=1 CLEAR=1
55 : # This finds the instantaneous location of the interface between the solid and liquid phases
56 : contour: FIND_CONTOUR_SURFACE ARG=dens CONTOUR=0.5 SEARCHDIR=dens_dist.x
57 : DUMPGRID ARG=contour FILE=contour.dat
58 : # This does the fourier transform of the location of the interface. We can extract the interfacial stiffness from the average of this fourier transform
59 : ft: FOURIER_TRANSFORM ARG=contour FT_TYPE=norm FOURIER_PARAMETERS=-1,1
60 : DUMPGRID ARG=ft FILE=fourier.dat STRIDE=10
61 : ```
62 :
63 : We do not think this action has been used in any other paper. If you are interested in using the functionality here we would recommend you carefully
64 : read the documentation for the FFTW library. You may find it necessary to modify the code in this action for your particular purpose.
65 :
66 : To what FFTW computes we would recommend reading [this page](http://www.fftw.org/doc/What-FFTW-Really-Computes.html#What-FFTW-Really-Computes).
67 :
68 : Notice that the keyword "FOURIER_PARAMETERS" specifies how the Fourier transform will be normalized. The keyword takes two numerical parameters $a$ and $b$ that are both set equal to one by default.
69 : The normalization is then defined as:
70 :
71 : $$
72 : \frac{1}{n^{(1-a)/2}} \sum_{j=0}^{n-1} X_j e^{2\pi b\, j k \sqrt{-1}/n}
73 : $$
74 :
75 : */
76 : //+ENDPLUMEDOC
77 :
78 :
79 : class FourierTransform : public gridtools::ActionWithGrid {
80 : private:
81 : bool firsttime;
82 : std::string output_type;
83 : bool real_output, store_norm;
84 : std::vector<int> fourier_params;
85 : gridtools::GridCoordinatesObject gridcoords;
86 : public:
87 : static void registerKeywords( Keywords& keys );
88 : explicit FourierTransform(const ActionOptions&ao);
89 : unsigned getNumberOfDerivatives() override ;
90 : const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
91 : std::vector<std::string> getGridCoordinateNames() const override ;
92 : void calculate() override ;
93 : };
94 :
95 : PLUMED_REGISTER_ACTION(FourierTransform,"FOURIER_TRANSFORM")
96 :
97 5 : void FourierTransform::registerKeywords( Keywords& keys ) {
98 5 : ActionWithGrid::registerKeywords( keys );
99 10 : keys.addInputKeyword("compulsory","ARG","grid","the label of the grid that you want to fourer transform");
100 5 : keys.add("optional","FT_TYPE","choose what kind of data you want as output on the grid. Possible values are: ABS = compute the complex modulus of Fourier coefficients (DEFAULT); NORM = compute the norm (i.e. ABS^2) of Fourier coefficients; COMPLEX = store the FFTW complex output on the grid (as a vector).");
101 5 : keys.add("compulsory","FOURIER_PARAMETERS","default","what kind of normalization is applied to the output and if the Fourier transform in FORWARD or BACKWARD. This keyword takes the form FOURIER_PARAMETERS=A,B, where A and B can be 0, 1 or -1. The default values are A=1 (no normalization at all) and B=1 (forward FFT). Other possible choices for A are: "
102 : "A=-1: normalize by the number of data, "
103 : "A=0: normalize by the square root of the number of data (one forward and followed by backward FFT recover the original data). ");
104 10 : keys.addOutputComponent("real","FT_TYPE","grid","the real part of the function");
105 10 : keys.addOutputComponent("imag","FT_TYPE","grid","the imaginary part of the function");
106 10 : keys.setValueDescription("grid","the fourier transform of the input grid");
107 5 : keys.addDOI("10.1088/1361-648X/aa893d");
108 5 : }
109 :
110 1 : FourierTransform::FourierTransform(const ActionOptions&ao):
111 : Action(ao),
112 : ActionWithGrid(ao),
113 1 : firsttime(true),
114 1 : real_output(true),
115 1 : store_norm(false),
116 1 : fourier_params(2) {
117 1 : if( getPntrToArgument(0)->getRank()!=2 ) {
118 0 : error("fourier transform currently only works with two dimensional grids");
119 : }
120 :
121 : // Get the type of FT
122 2 : parse("FT_TYPE",output_type);
123 1 : if (output_type.length()==0) {
124 0 : log<<" keyword FT_TYPE unset. By default output grid will contain REAL Fourier coefficients\n";
125 1 : } else if ( output_type=="ABS" || output_type=="abs") {
126 0 : log << " keyword FT_TYPE is '"<< output_type << "' : will compute the MODULUS of Fourier coefficients\n";
127 1 : } else if ( output_type=="NORM" || output_type=="norm") {
128 0 : log << " keyword FT_TYPE is '"<< output_type << "' : will compute the NORM of Fourier coefficients\n";
129 0 : store_norm=true;
130 1 : } else if ( output_type=="COMPLEX" || output_type=="complex" ) {
131 1 : log<<" keyword FT_TYPE is '"<< output_type <<"' : output grid will contain the COMPLEX Fourier coefficients\n";
132 1 : real_output=false;
133 : } else {
134 0 : error("keyword FT_TYPE unrecognized!");
135 : }
136 :
137 : // Normalize output?
138 : std::string params_str;
139 1 : parse("FOURIER_PARAMETERS",params_str);
140 1 : if (params_str=="default") {
141 0 : fourier_params.assign( fourier_params.size(), 1 );
142 0 : log.printf(" default values of Fourier parameters A=%i, B=%i : the output will NOT be normalized and BACKWARD Fourier transform is computed \n", fourier_params[0],fourier_params[1]);
143 : } else {
144 1 : std::vector<std::string> fourier_str = Tools::getWords(params_str, "\t\n ,");
145 1 : if (fourier_str.size()>2) {
146 0 : error("FOURIER_PARAMETERS can take just two values");
147 : }
148 3 : for (unsigned i=0; i<fourier_str.size(); ++i) {
149 2 : Tools::convert(fourier_str[i],fourier_params[i]);
150 2 : if (fourier_params[i]>1 || fourier_params[i]<-1) {
151 0 : error("values accepted for FOURIER_PARAMETERS are only -1, 1 or 0");
152 : }
153 : }
154 1 : log.printf(" Fourier parameters are A=%i, B=%i \n", fourier_params[0],fourier_params[1]);
155 1 : }
156 :
157 1 : std::vector<std::size_t> shape( getPntrToArgument(0)->getRank() );
158 1 : if (real_output) {
159 0 : addValueWithDerivatives( shape );
160 : } else {
161 1 : addComponentWithDerivatives( "real", shape );
162 2 : addComponentWithDerivatives( "imag", shape );
163 : }
164 :
165 1 : unsigned dimension = getPntrToArgument(0)->getRank();
166 1 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
167 1 : if( !ag ) {
168 0 : error("input action should be a grid");
169 : }
170 1 : const gridtools::GridCoordinatesObject & gcoords( ag->getGridCoordinatesObject() );
171 1 : if( gcoords.getGridType()=="fibonacci" ) {
172 0 : error("cannot fourier transform fibonacci grids");
173 : }
174 1 : std::vector<bool> ipbc( dimension );
175 3 : for(unsigned i=0; i<dimension; ++i) {
176 2 : ipbc[i] = gcoords.isPeriodic(i);
177 : }
178 1 : gridcoords.setup( "flat", ipbc, 0, 0.0 );
179 1 : checkRead();
180 : #ifndef __PLUMED_HAS_FFTW
181 : error("this feature is only available if you compile PLUMED with FFTW");
182 : #endif
183 2 : }
184 :
185 4 : unsigned FourierTransform::getNumberOfDerivatives() {
186 4 : return 2;
187 : }
188 :
189 7 : const gridtools::GridCoordinatesObject& FourierTransform::getGridCoordinatesObject() const {
190 7 : return gridcoords;
191 : }
192 :
193 2 : std::vector<std::string> FourierTransform::getGridCoordinateNames() const {
194 2 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
195 : // already checked in the constructor
196 2 : plumed_assert( ag );
197 2 : return ag->getGridCoordinateNames();
198 : }
199 :
200 1 : void FourierTransform::calculate() {
201 1 : if( firsttime ) {
202 1 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
203 : // already checked in the constructor
204 1 : plumed_assert( ag );
205 1 : const gridtools::GridCoordinatesObject & gcoords( ag->getGridCoordinatesObject() );
206 1 : std::vector<double> fspacing;
207 1 : std::vector<std::size_t> snbins( getGridCoordinatesObject().getDimension() );
208 1 : std::vector<std::string> smin( gcoords.getDimension() ), smax( gcoords.getDimension() );
209 3 : for(unsigned i=0; i<getGridCoordinatesObject().getDimension(); ++i) {
210 4 : smin[i]=gcoords.getMin()[i];
211 4 : smax[i]=gcoords.getMax()[i];
212 : // Compute k-grid extents
213 : double dmin, dmax;
214 2 : snbins[i]=gcoords.getNbin(false)[i];
215 2 : Tools::convert(smin[i],dmin);
216 2 : Tools::convert(smax[i],dmax);
217 2 : dmax=2.0*pi*snbins[i]/( dmax - dmin );
218 2 : dmin=0.0;
219 2 : Tools::convert(dmin,smin[i]);
220 2 : Tools::convert(dmax,smax[i]);
221 : }
222 1 : gridcoords.setBounds( smin, smax, snbins, fspacing );
223 1 : firsttime=false;
224 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
225 2 : getPntrToComponent(i)->setShape( gcoords.getNbin(true) );
226 : }
227 1 : }
228 :
229 : #ifdef __PLUMED_HAS_FFTW
230 : // *** CHECK CORRECT k-GRID BOUNDARIES ***
231 : //log<<"Real grid boundaries: \n"
232 : // <<" min_x: "<<mygrid->getMin()[0]<<" min_y: "<<mygrid->getMin()[1]<<"\n"
233 : // <<" max_x: "<<mygrid->getMax()[0]<<" max_y: "<<mygrid->getMax()[1]<<"\n"
234 : // <<"K-grid boundaries:"<<"\n"
235 : // <<" min_x: "<<ft_min[0]<<" min_y: "<<ft_min[1]<<"\n"
236 : // <<" max_x: "<<ft_max[0]<<" max_y: "<<ft_max[1]<<"\n";
237 :
238 : // Get the size of the input data arrays (to allocate FFT data)
239 1 : std::vector<std::size_t> N_input_data( gridcoords.getNbin(true) );
240 : size_t fft_dimension=1;
241 3 : for(unsigned i=0; i<N_input_data.size(); ++i) {
242 2 : fft_dimension*=static_cast<size_t>( N_input_data[i] );
243 : }
244 : // FFT arrays
245 1 : std::vector<std::complex<double> > input_data(fft_dimension), fft_data(fft_dimension);
246 :
247 : // Fill real input with the data on the grid
248 : Value* arg=getPntrToArgument(0);
249 1 : std::vector<unsigned> ind( arg->getRank() );
250 10202 : for (unsigned i=0; i<arg->getNumberOfValues(); ++i) {
251 : // Get point indices
252 10201 : gridcoords.getIndices(i, ind);
253 : // Fill input data in row-major order
254 10201 : input_data[ind[0]*N_input_data[0]+ind[1]].real( arg->get( i ) );
255 10201 : input_data[ind[0]*N_input_data[0]+ind[1]].imag( 0.0 );
256 : }
257 :
258 : // *** HERE is the only clear limitation: I'm computing explicitly a 2D FT. It should not happen to deal with other than two-dimensional grid ...
259 1 : fftw_plan plan_complex = fftw_plan_dft_2d(N_input_data[0], N_input_data[1], reinterpret_cast<fftw_complex*>(&input_data[0]), reinterpret_cast<fftw_complex*>(&fft_data[0]), fourier_params[1], FFTW_ESTIMATE);
260 :
261 : // Compute FT
262 1 : fftw_execute( plan_complex );
263 :
264 : // Compute the normalization constant
265 : double norm=1.0;
266 3 : for (unsigned i=0; i<N_input_data.size(); ++i) {
267 2 : norm *= pow( N_input_data[i], (1-fourier_params[0])/2 );
268 : }
269 :
270 : // Save FT data to output grid
271 1 : std::vector<std::size_t> N_out_data ( getGridCoordinatesObject().getNbin(true) );
272 1 : std::vector<unsigned> out_ind ( getPntrToArgument(0)->getRank() );
273 10202 : for(unsigned i=0; i<getPntrToArgument(0)->getNumberOfValues(); ++i) {
274 10201 : gridcoords.getIndices( i, out_ind );
275 10201 : if (real_output) {
276 : double ft_value;
277 : // Compute abs/norm and fix normalization
278 0 : if (!store_norm) {
279 0 : ft_value=std::abs( fft_data[out_ind[0]*N_out_data[0]+out_ind[1]] / norm );
280 : } else {
281 0 : ft_value=std::norm( fft_data[out_ind[0]*N_out_data[0]+out_ind[1]] / norm );
282 : }
283 : // Set the value
284 0 : getPntrToComponent(0)->set( i, ft_value);
285 : } else {
286 : double ft_value_real, ft_value_imag;
287 10201 : ft_value_real=fft_data[out_ind[0]*N_out_data[0]+out_ind[1]].real() / norm;
288 10201 : ft_value_imag=fft_data[out_ind[0]*N_out_data[0]+out_ind[1]].imag() / norm;
289 : // Set values
290 10201 : getPntrToComponent(0)->set( i, ft_value_real );
291 10201 : getPntrToComponent(1)->set( i, ft_value_imag );
292 : }
293 : }
294 :
295 : // Free FFTW stuff
296 1 : fftw_destroy_plan(plan_complex);
297 : #endif
298 2 : }
299 :
300 : } // end namespace 'gridtools'
301 : } // end namespace 'PLMD'
|
