From 377c6ce7ab18c5c470623f237240a24010cc5371 Mon Sep 17 00:00:00 2001
From: Alexey Shvetsov <alexxy@gentoo.org>
Date: Wed, 11 Jul 2018 14:42:54 +0300
Subject: [PATCH] Fix gaus

---
 src/sans.cpp | 84 +++++++++++++++++++++++++++++-----------------------
 1 file changed, 47 insertions(+), 37 deletions(-)

diff --git a/src/sans.cpp b/src/sans.cpp
index 354d9d4..7ff2645 100644
--- a/src/sans.cpp
+++ b/src/sans.cpp
@@ -53,7 +53,7 @@ class SANS : public TrajectoryAnalysisModule
         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,
@@ -65,7 +65,7 @@ class SANS : public TrajectoryAnalysisModule
     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_;
@@ -76,8 +76,8 @@ class SANS : public TrajectoryAnalysisModule
         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()
@@ -87,12 +87,12 @@ 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."
@@ -102,19 +102,19 @@ SANS::initOptions(IOptionsContainer          *options,
     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_ = &top;
@@ -123,50 +123,60 @@ SANS::initAnalysis(const TrajectoryAnalysisSettings &settings,
 
 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]);
             }
         }
     }
-- 
2.22.0