virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top);
virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr);
+ const t_trxframe &fr);
virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
private:
class ModuleData;
- double Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing);
+ double Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing);
std::string fnNdx_;
double cutoff_;
RVec gridSpacing_;
AnalysisNeighborhood nb_;
- const TopologyInformation *top_;
- std::vector<std::vector<std::vector<double>>> gausGrid_;
+ const TopologyInformation *top_;
+ std::vector < std::vector < std::vector<double>>> gausGrid_;
};
SANS::SANS()
double SANS::Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing)
{
- 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.);
+ 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.);
}
void
SANS::initOptions(IOptionsContainer *options,
- TrajectoryAnalysisSettings *settings)
+ TrajectoryAnalysisSettings *settings)
{
static const char *const desc[] = {
"Analysis tool for finding molecular core."
settings->setHelpText(desc);
// Add option for cutoff constant
options->addOption(DoubleOption("cutoff")
- .store(&cutoff_)
- .description("cutoff for neighbour search"));
+ .store(&cutoff_)
+ .description("cutoff for neighbour search"));
options->addOption(DoubleOption("grid")
- .store(&grid_)
- .description("grid spacing for gaus grid"));
+ .store(&grid_)
+ .description("grid spacing for gaus grid"));
options->addOption(DoubleOption("boxscale")
- .store(&boxScale_)
- .description("box scaling for gaus grid"));
+ .store(&boxScale_)
+ .description("box scaling for gaus grid"));
}
void
-SANS::initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation & top)
+SANS::initAnalysis(const TrajectoryAnalysisSettings &settings,
+ const TopologyInformation &top)
{
nb_.setCutoff(cutoff_);
top_ = ⊤
void
SANS::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr)
+ const t_trxframe &fr)
{
// set gridpoints and grid spacing for orthogonal boxes only....
- for (int i = 0; i<DIM;i++) {
+ for (int i = 0; i < DIM; i++)
+ {
gridSpacing_[i] = grid_;
- gridPoints_[i] = static_cast<int>(std::ceil(fr.box[i][i]*boxScale_/grid_));
+ gridPoints_[i] = static_cast<int>(std::ceil(fr.box[i][i]*boxScale_/grid_));
}
// prepare gausGrid
- gausGrid_.resize(gridPoints_[XX], std::vector<std::vector<double>>(gridPoints_[YY], std::vector<double>(gridPoints_[ZZ])));
+ gausGrid_.resize(gridPoints_[XX], std::vector < std::vector < double>>(gridPoints_[YY], std::vector<double>(gridPoints_[ZZ])));
}
void
SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata)
+ TrajectoryAnalysisModuleData *pdata)
{
RVec GridSpacing(0.1, 0.1, 0.1);
AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, AnalysisNeighborhoodPositions(fr.x, fr.natoms));
- t_topology *top = top_->topology();
+ t_topology *top = top_->topology();
- fprintf(stderr,"Box dimensions:\n");
- for(int i = 0; i<DIM; i++) {
- fprintf(stderr,"[ ");
- for(int j = 0 ; j<DIM; j++) {
- fprintf(stderr," %f ", fr.box[i][j]);
+ fprintf(stderr, "Box dimensions:\n");
+ for (int i = 0; i < DIM; i++)
+ {
+ fprintf(stderr, "[ ");
+ for (int j = 0; j < DIM; j++)
+ {
+ fprintf(stderr, " %f ", fr.box[i][j]);
}
fprintf(stderr, " ]\n");
}
// main loop over grid points
- for (int i=0; i<gridPoints_[XX];i++){
- for(int j=0; j<gridPoints_[YY];j++) {
- for (int k=0; k<gridPoints_[ZZ];k++) {
- RVec point(i*gridSpacing_[XX], j*gridSpacing_[YY], k*gridSpacing_[ZZ]);
+ for (int i = 0; i < gridPoints_[XX]; i++)
+ {
+ for (int j = 0; j < gridPoints_[YY]; j++)
+ {
+ for (int k = 0; k < gridPoints_[ZZ]; k++)
+ {
+ RVec point(i*gridSpacing_[XX], j*gridSpacing_[YY], k*gridSpacing_[ZZ]);
AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(point.as_vec());
- AnalysisNeighborhoodPair pair;
- while (pairSearch.findNextPair(&pair)) {
- // fprintf(stderr,"Index %d\n", pair.refIndex());
- // fprintf(stderr,"dx = (%f, %f, %f)\n", pair.dx()[XX], pair.dx()[YY], pair.dx()[ZZ] );
- // fprintf(stderr,"ref = (%f, %f, %f)\n", fr.x[pair.refIndex()][XX], fr.x[pair.refIndex()][YY], fr.x[pair.refIndex()][ZZ]);
- // fprintf(stderr,"ref_dx = (%f, %f, %f)\n", GridSpacing[XX]+pair.dx()[XX], GridSpacing[YY]+pair.dx()[YY], GridSpacing[ZZ]+pair.dx()[ZZ] );
+ AnalysisNeighborhoodPair pair;
+ while (pairSearch.findNextPair(&pair))
+ {
+ // fprintf(stderr,"point = [ %f %f %f ]\n", point[XX],point[YY],point[ZZ]);
+ // fprintf(stderr,"Index %d\n", pair.refIndex());
+ // fprintf(stderr,"dx = (%f, %f, %f)\n", pair.dx()[XX], pair.dx()[YY], pair.dx()[ZZ] );
+ // fprintf(stderr,"ref = (%f, %f, %f)\n", fr.x[pair.refIndex()][XX], fr.x[pair.refIndex()][YY], fr.x[pair.refIndex()][ZZ]);
+ // fprintf(stderr,"ref_dx = (%f, %f, %f)\n", GridSpacing[XX]+pair.dx()[XX], GridSpacing[YY]+pair.dx()[YY], GridSpacing[ZZ]+pair.dx()[ZZ] );
RVec refPos(point[XX]+pair.dx()[XX], point[YY]+pair.dx()[YY], point[ZZ]+pair.dx()[ZZ]);
- gausGrid_[i][j][k] += Gaussian(top->atoms.atom[pair.refIndex()].m, 2.0, refPos, point, gridSpacing_);
+ // fprintf(stderr,"refPos = [ %f %f %f ]\n", refPos[XX], refPos[YY], refPos[ZZ]);
+ // fprintf(stderr, "Mass = %f, Radius = %f\n", top->atoms.atom[pair.refIndex()].m, top->atomtypes.radius[top->atoms.atom[pair.refIndex()].type]);
+ gausGrid_[i][j][k] += Gaussian(top->atoms.atom[pair.refIndex()].m, 2., refPos, point, gridSpacing_);
}
- fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i,j,k, gausGrid_[i][j][k]);
+ fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i, j, k, gausGrid_[i][j][k]);
}
}
}