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>
44 * Template class to serve as a basis for user analysis tools.
46 class SANS : public TrajectoryAnalysisModule
51 virtual void initOptions(IOptionsContainer *options,
52 TrajectoryAnalysisSettings *settings);
53 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
54 const TopologyInformation &top);
55 virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
56 const t_trxframe &fr);
59 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
60 TrajectoryAnalysisModuleData *pdata);
62 virtual void finishAnalysis(int nframes);
63 virtual void writeOutput();
68 double Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing);
78 AnalysisNeighborhood nb_;
79 const TopologyInformation *top_;
80 std::vector<std::vector<std::vector<double>>> gausGrid_;
84 : cutoff_(1.0), grid_(0.05), boxScale_(1.0)
88 double SANS::Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing)
90 return 2.*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*216.);
94 SANS::initOptions(IOptionsContainer *options,
95 TrajectoryAnalysisSettings *settings)
97 static const char *const desc[] = {
98 "Analysis tool for finding molecular core."
101 // Add the descriptive text (program help text) to the options
102 settings->setHelpText(desc);
103 // Add option for cutoff constant
104 options->addOption(DoubleOption("cutoff")
106 .description("cutoff for neighbour search"));
107 options->addOption(DoubleOption("grid")
109 .description("grid spacing for gaus grid"));
110 options->addOption(DoubleOption("boxscale")
112 .description("box scaling for gaus grid"));
116 SANS::initAnalysis(const TrajectoryAnalysisSettings &settings,
117 const TopologyInformation & top)
119 nb_.setCutoff(cutoff_);
125 SANS::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
126 const t_trxframe &fr)
128 // set gridpoints and grid spacing for orthogonal boxes only....
129 for (int i = 0; i<DIM;i++) {
130 gridSpacing_[i] = grid_;
131 gridPoints_[i] = static_cast<int>(std::ceil(fr.box[i][i]*boxScale_/grid_));
134 gausGrid_.resize(gridPoints_[XX], std::vector<std::vector<double>>(gridPoints_[YY], std::vector<double>(gridPoints_[ZZ])));
138 SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
139 TrajectoryAnalysisModuleData *pdata)
141 RVec GridSpacing(0.1, 0.1, 0.1);
142 AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, AnalysisNeighborhoodPositions(fr.x, fr.natoms));
144 t_topology *top = top_->topology();
146 fprintf(stderr,"Box dimensions:\n");
147 for(int i = 0; i<DIM; i++) {
148 fprintf(stderr,"[ ");
149 for(int j = 0 ; j<DIM; j++) {
150 fprintf(stderr," %f ", fr.box[i][j]);
152 fprintf(stderr, " ]\n");
154 // main loop over grid points
155 for (int i=0; i<gridPoints_[XX];i++){
156 for(int j=0; j<gridPoints_[YY];j++) {
157 for (int k=0; k<gridPoints_[ZZ];k++) {
158 RVec point(i*gridSpacing_[XX], j*gridSpacing_[YY], k*gridSpacing_[ZZ]);
159 AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(point.as_vec());
160 AnalysisNeighborhoodPair pair;
161 while (pairSearch.findNextPair(&pair)) {
162 // fprintf(stderr,"Index %d\n", pair.refIndex());
163 // fprintf(stderr,"dx = (%f, %f, %f)\n", pair.dx()[XX], pair.dx()[YY], pair.dx()[ZZ] );
164 // fprintf(stderr,"ref = (%f, %f, %f)\n", fr.x[pair.refIndex()][XX], fr.x[pair.refIndex()][YY], fr.x[pair.refIndex()][ZZ]);
165 // fprintf(stderr,"ref_dx = (%f, %f, %f)\n", GridSpacing[XX]+pair.dx()[XX], GridSpacing[YY]+pair.dx()[YY], GridSpacing[ZZ]+pair.dx()[ZZ] );
166 RVec refPos(point[XX]+pair.dx()[XX], point[YY]+pair.dx()[YY], point[ZZ]+pair.dx()[ZZ]);
167 gausGrid_[i][j][k] += Gaussian(top->atoms.atom[pair.refIndex()].m, 2.0, refPos, point, gridSpacing_);
169 fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i,j,k, gausGrid_[i][j][k]);
173 // only for orthogonal boxes...
174 fprintf(stderr, "NX = %d, NY = %d, NZ = %d\n", gridPoints_[XX], gridPoints_[YY], gridPoints_[ZZ]);
178 SANS::finishAnalysis(int nframes)
189 * The main function for the analysis template.
192 main(int argc, char *argv[])
194 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SANS>(argc, argv);