2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS 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 GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include <gromacs/trajectoryanalysis.h>
38 #include <gromacs/selection/nbsearch.h>
39 #include <gromacs/math/vec.h>
40 #include <gromacs/math/units.h>
41 #include <gromacs/topology/atomprop.h>
47 * Template class to serve as a basis for user analysis tools.
49 class SANS : public TrajectoryAnalysisModule
54 virtual void initOptions(IOptionsContainer *options,
55 TrajectoryAnalysisSettings *settings);
56 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
57 const TopologyInformation &top);
58 virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
59 const t_trxframe &fr);
62 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
63 TrajectoryAnalysisModuleData *pdata);
65 virtual void finishAnalysis(int nframes);
66 virtual void writeOutput();
71 double Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing);
80 AnalysisNeighborhood nb_;
81 const TopologyInformation *top_;
82 std::vector< std::vector < std::vector<double>>> gausGrid_;
83 std::vector<std::vector<std::vector<bool>>> isSolvent_;
84 std::vector<double> vdw_radius_;
88 : cutoff_(1.0), grid_(0.05), boxScale_(1.0)
92 double SANS::Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing)
94 //return value*M_PI*M_PI*GridSpacing[XX]*GridSpacing[YY]*GridSpacing[ZZ]*exp(-distance2(Rpos, Rgrid)/(2*sigma*sigma))/(std::sqrt(2.0)*sigma*sigma*sigma*108.);
95 //return value*std::exp(-distance2(Rpos, Rgrid)/(2*sigma*sigma))/(2*M_PI*std::sqrt(2.0*M_PI)*sigma*sigma*sigma);
96 return value*std::exp(-distance2(Rpos, Rgrid)/(2*sigma*sigma))/(2*M_PI*std::sqrt(2.0*M_PI)*sigma*sigma*sigma)*GridSpacing[XX]*GridSpacing[YY]*GridSpacing[ZZ];
97 //return value*std::exp(-distance2(Rpos, Rgrid)/(sigma*sigma))/(2*M_PI*std::sqrt(2.0*M_PI)*sigma*sigma*sigma)*GridSpacing[XX]*GridSpacing[YY]*GridSpacing[ZZ];
101 SANS::initOptions(IOptionsContainer *options,
102 TrajectoryAnalysisSettings *settings)
104 static const char *const desc[] = {
105 "Analysis tool for finding molecular core."
108 // Add the descriptive text (program help text) to the options
109 settings->setHelpText(desc);
110 // Add option for cutoff constant
111 options->addOption(DoubleOption("cutoff")
113 .description("cutoff for neighbour search"));
114 options->addOption(DoubleOption("grid")
116 .description("grid spacing for gaus grid"));
117 options->addOption(DoubleOption("boxscale")
119 .description("box scaling for gaus grid"));
123 SANS::initAnalysis(const TrajectoryAnalysisSettings &settings,
124 const TopologyInformation &top)
127 gmx_atomprop_t aps = gmx_atomprop_init();
128 t_atoms *atoms = &(top.topology()->atoms);
131 nb_.setCutoff(cutoff_);
134 for (int i = 0; i < top.topology()->atoms.nr; i++)
136 // Lookup the Van der Waals radius of this atom
137 int resnr = atoms->atom[i].resind;
138 if (gmx_atomprop_query(aps, epropVDW,
139 *(atoms->resinfo[resnr].name),
140 *(atoms->atomname[i]),
143 vdw_radius_.push_back(value);
147 fprintf(stderr, "Could not determine VDW radius for %s-%s. Set to zero.\n",
148 *(atoms->resinfo[resnr].name),
149 *(atoms->atomname[i]));
150 vdw_radius_.push_back(0.0);
158 SANS::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
159 const t_trxframe &fr)
161 // set gridpoints and grid spacing for orthogonal boxes only....
162 for (int i = 0; i < DIM; i++)
164 gridSpacing_[i] = grid_;
165 gridPoints_[i] = static_cast<int>(std::ceil(fr.box[i][i]*boxScale_/grid_));
168 gausGrid_.resize(gridPoints_[XX], std::vector < std::vector < double>>(gridPoints_[YY], std::vector<double>(gridPoints_[ZZ])));
169 isSolvent_.resize(gridPoints_[XX], std::vector < std::vector < bool>>(gridPoints_[YY], std::vector<bool>(gridPoints_[ZZ])));
173 SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
174 TrajectoryAnalysisModuleData *pdata)
176 RVec GridSpacing(0.1, 0.1, 0.1);
177 AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, AnalysisNeighborhoodPositions(fr.x, fr.natoms));
179 t_topology *top = top_->topology();
181 fprintf(stderr, "Box dimensions:\n");
182 for (int i = 0; i < DIM; i++)
184 fprintf(stderr, "[ ");
185 for (int j = 0; j < DIM; j++)
187 fprintf(stderr, " %f ", fr.box[i][j]);
189 fprintf(stderr, " ]\n");
191 // main loop over grid points
192 for (int i = 0; i < gridPoints_[XX]; i++)
194 for (int j = 0; j < gridPoints_[YY]; j++)
196 for (int k = 0; k < gridPoints_[ZZ]; k++)
198 RVec point(i*gridSpacing_[XX], j*gridSpacing_[YY], k*gridSpacing_[ZZ]);
199 AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(point.as_vec());
200 AnalysisNeighborhoodPair pair;
201 while (pairSearch.findNextPair(&pair))
203 //fprintf(stderr,"point = [ %f %f %f ]\n", point[XX],point[YY],point[ZZ]);
204 // fprintf(stderr,"Index %d\n", pair.refIndex());
205 // fprintf(stderr,"dx = (%f, %f, %f)\n", pair.dx()[XX], pair.dx()[YY], pair.dx()[ZZ] );
206 // fprintf(stderr,"ref = (%f, %f, %f)\n", fr.x[pair.refIndex()][XX], fr.x[pair.refIndex()][YY], fr.x[pair.refIndex()][ZZ]);
207 // fprintf(stderr,"ref_dx = (%f, %f, %f)\n", GridSpacing[XX]+pair.dx()[XX], GridSpacing[YY]+pair.dx()[YY], GridSpacing[ZZ]+pair.dx()[ZZ] );
208 RVec refPos(point[XX]+pair.dx()[XX], point[YY]+pair.dx()[YY], point[ZZ]+pair.dx()[ZZ]);
209 //fprintf(stderr,"refPos = [ %f %f %f ]\n", refPos[XX], refPos[YY], refPos[ZZ]);
210 //fprintf(stderr, "Mass = %f, Radius = %f\n", top->atoms.atom[pair.refIndex()].m, top->atomtypes.radius[top->atoms.atom[pair.refIndex()].type]);
211 //fprintf(stderr, "Distance^2 = %3.8f\n", pair.distance2());
212 //fprintf(stderr, "Gausian = %3.8f\n",Gaussian(top->atoms.atom[pair.refIndex()].m, 0.2, refPos, point, gridSpacing_)*AMU/((NANO*NANO*NANO)*(gridSpacing_[XX]*gridSpacing_[YY]*gridSpacing_[ZZ])));
213 //gausGrid_[i][j][k] += Gaussian(top->atoms.atom[pair.refIndex()].m, 0.1, refPos, point, gridSpacing_)*AMU/((NANO*NANO*NANO)*(gridSpacing_[XX]*gridSpacing_[YY]*gridSpacing_[ZZ]));
214 gausGrid_[i][j][k] += Gaussian(top->atoms.atom[pair.refIndex()].m, vdw_radius_[pair.refIndex()], refPos, point, gridSpacing_)*AMU/((NANO*NANO*NANO)*(gridSpacing_[XX]*gridSpacing_[YY]*gridSpacing_[ZZ]));
215 // Mark if this is a pure solvent point e.g. no non-solvent atoms within cut-off
218 fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i, j, k, gausGrid_[i][j][k]);
222 // only for orthogonal boxes...
223 fprintf(stderr, "NX = %d, NY = %d, NZ = %d\n", gridPoints_[XX], gridPoints_[YY], gridPoints_[ZZ]);
227 SANS::finishAnalysis(int nframes)
235 // Construct opendx grid
236 fo = fopen("data.dx", "w");
237 //object 1 class gridpositions counts 140 140 140
238 //origin 0.000000 0.000000 0.000000
239 //delta 0.500000 0.000000 0.000000
240 //delta 0.000000 0.500000 0.000000
241 //delta 0.000000 0.000000 0.500000
242 //object 2 class gridconnections counts 140 140 140
243 //object 3 class array type "double" rank 0 items 2744000 data follows
244 fprintf(fo, "object 1 class gridpositions counts %d %d %d\n", gridPoints_[XX], gridPoints_[YY], gridPoints_[ZZ]);
245 fprintf(fo, "origin 0.000000 0.000000 0.000000\n");
246 fprintf(fo, "delta %3.6f %3.6f %3.6f\n", gridSpacing_[XX]*NM2A, 0.0, 0.0);
247 fprintf(fo, "delta %3.6f %3.6f %3.6f\n", 0.0, gridSpacing_[YY]*NM2A, 0.0);
248 fprintf(fo, "delta %3.6f %3.6f %3.6f\n", 0.0, 0.0, gridSpacing_[ZZ]*NM2A);
249 fprintf(fo, "object 2 class gridconnections counts %d %d %d\n", gridPoints_[XX], gridPoints_[YY], gridPoints_[ZZ]);
250 fprintf(fo, "object 3 class array type \"double\" rank 0 items %d data follows\n", gridPoints_[XX]*gridPoints_[YY]*gridPoints_[ZZ]);
251 // now dump data in ordered mode
252 // u(0,0,0) u(0,0,1) u(0,0,2)
254 // u(0,0,nz-3) u(0,0,nz-2) u(0,0,nz-1)
255 // u(0,1,0) u(0,1,1) u(0,1,2)
257 // u(0,1,nz-3) u(0,1,nz-2) u(0,1,nz-1)
259 // u(0,ny-1,nz-3) u(0,ny-1,nz-2) u(0,ny-1,nz-1)
260 // u(1,0,0) u(1,0,1) u(1,0,2)
264 for (int i = 0; i < gridPoints_[XX]; i++)
266 for (int j = 0; j < gridPoints_[YY]; j++)
268 for (int k = 0; k < gridPoints_[ZZ]; k++)
270 fprintf(fo, "%3.8f ", gausGrid_[i][j][k]);
284 //attribute "dep" string "positions"
285 //object "density" class field
286 //component "positions" value 1
287 //component "connections" value 2
288 //component "data" value 3
289 fprintf(fo, "attribute \"dep\" string \"positions\"\n");
290 fprintf(fo, "object \"density\" class field\n");
291 fprintf(fo, "component \"positions\" value 1\n");
292 fprintf(fo, "component \"connections\" value 2\n");
293 fprintf(fo, "component \"data\" value 3\n");
298 * The main function for the analysis template.
301 main(int argc, char *argv[])
303 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SANS>(argc, argv);