Working gausian grid...

diff --git a/src/sans.cpp b/src/sans.cpp
index bfd4d0efd95f8477239eb348d3ca6aaa49534dfb..354d9d46f7c2901a923923fedd5f2b91e1901043 100644 (file)
--- a/src/sans.cpp
+++ b/src/sans.cpp
@@ -52,6 +52,9 @@ class SANS : public TrajectoryAnalysisModule
                                  TrajectoryAnalysisSettings *settings);
         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
                                   const TopologyInformation        &top);
+        virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
+                                                 const t_trxframe                 &fr);
+
 
         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                                   TrajectoryAnalysisModuleData *pdata);
@@ -67,14 +70,18 @@ class SANS : public TrajectoryAnalysisModule
         std::string                      fnNdx_;
         double                           cutoff_;
         double                           grid_;
+        double                           boxScale_;
 
-        AnalysisNeighborhood             nb_;
-        const TopologyInformation          *top_;
+        IVec gridPoints_;
+        RVec gridSpacing_;
 
+        AnalysisNeighborhood             nb_;
+        const TopologyInformation      *top_;
+        std::vector<std::vector<std::vector<double>>> gausGrid_;
 };
 
 SANS::SANS()
-    : cutoff_(1.0), grid_(0.05)
+    : cutoff_(1.0), grid_(0.05), boxScale_(1.0)
 {
 }
 
@@ -93,6 +100,16 @@ SANS::initOptions(IOptionsContainer          *options,
 
     // Add the descriptive text (program help text) to the options
     settings->setHelpText(desc);
+    // Add option for cutoff constant
+    options->addOption(DoubleOption("cutoff")
+                                .store(&cutoff_)
+                                .description("cutoff for neighbour search"));
+    options->addOption(DoubleOption("grid")
+                       .store(&grid_)
+                       .description("grid spacing for gaus grid"));
+    options->addOption(DoubleOption("boxscale")
+                       .store(&boxScale_)
+                       .description("box scaling for gaus grid"));
 }
 
 void
@@ -101,6 +118,20 @@ SANS::initAnalysis(const TrajectoryAnalysisSettings &settings,
 {
     nb_.setCutoff(cutoff_);
     top_ = &top;
+
+}
+
+void
+SANS::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
+                         const t_trxframe                 &fr)
+{
+    // set gridpoints and grid spacing for orthogonal boxes only....
+    for (int i = 0; i<DIM;i++) {
+        gridSpacing_[i] = 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])));
 }
 
 void
@@ -109,17 +140,38 @@ SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 {
     RVec GridSpacing(0.1, 0.1, 0.1);
     AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, AnalysisNeighborhoodPositions(fr.x, fr.natoms));
-    AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(GridSpacing.as_vec());
-    AnalysisNeighborhoodPair pair;
+
     t_topology *top = top_->topology();
-    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] );
-
-      }
+
+    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]);
+                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] );
+                    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, "GausGrid[%d, %d, %d] = %f\n", i,j,k, gausGrid_[i][j][k]);
+            }
+        }
+    }
+    // only for orthogonal boxes...
+    fprintf(stderr, "NX = %d, NY = %d, NZ = %d\n", gridPoints_[XX], gridPoints_[YY], gridPoints_[ZZ]);
 }
 
 void