Normalization to volume, vdw radii as sigma and opendx output

Signed-off-by: Alexey Shvetsov <alexxy@gentoo.org>

diff --git a/src/sans.cpp b/src/sans.cpp
index 4c932d9c6684d4aa698e77a7deebbf89669cc370..5c13d802c3843c21ee659dbc5f6dbc1451ffb0fc 100644 (file)
--- a/src/sans.cpp
+++ b/src/sans.cpp
@@ -38,6 +38,8 @@
 #include <gromacs/selection/nbsearch.h>
 #include <gromacs/math/vec.h>
 #include <gromacs/math/units.h>
+#include <gromacs/topology/atomprop.h>
+
 
 using namespace gmx;
 
@@ -68,7 +70,6 @@ class SANS : public TrajectoryAnalysisModule
 
         double Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing);
 
-        std::string                      fnNdx_;
         double                           cutoff_;
         double                           grid_;
         double                           boxScale_;
@@ -79,6 +80,7 @@ class SANS : public TrajectoryAnalysisModule
         AnalysisNeighborhood             nb_;
         const TopologyInformation       *top_;
         std::vector < std::vector < std::vector<double>>> gausGrid_;
+        std::vector<double>               vdw_radius_;
 };
 
 SANS::SANS()
@@ -91,6 +93,7 @@ double SANS::Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec Gr
     //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.);
     //return value*std::exp(-distance2(Rpos, Rgrid)/(2*sigma*sigma))/(2*M_PI*std::sqrt(2.0*M_PI)*sigma*sigma*sigma);
     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];
+    //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];
 }
 
 void
@@ -121,6 +124,28 @@ SANS::initAnalysis(const TrajectoryAnalysisSettings  &settings,
 {
     nb_.setCutoff(cutoff_);
     top_ = &top;
+    gmx_atomprop_t aps    = gmx_atomprop_init();
+    t_atoms       *atoms  = &(top.topology()->atoms);
+    real value;
+    for(int i=0; i<top.topology()->atoms.nr;i++) {
+        // Lookup the Van der Waals radius of this atom
+        int resnr = atoms->atom[i].resind;
+        if (gmx_atomprop_query(aps, epropVDW,
+                               *(atoms->resinfo[resnr].name),
+                               *(atoms->atomname[i]),
+                               &value))
+            {
+                vdw_radius_.push_back(value);
+            }
+            else
+            {
+                fprintf(stderr, "Could not determine VDW radius for %s-%s. Set to zero.\n",
+                                *(atoms->resinfo[resnr].name),
+                                *(atoms->atomname[i]));
+                    vdw_radius_.push_back(0.0);
+            }
+
+    }
 
 }
 
@@ -169,16 +194,18 @@ SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                 AnalysisNeighborhoodPair       pair;
                 while (pairSearch.findNextPair(&pair))
                 {
-                    // fprintf(stderr,"point = [ %f %f %f ]\n", point[XX],point[YY],point[ZZ]);
+                    //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]);
-                    // 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_)*AMU/(NANO*NANO*NANO);
-
+                    //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]);
+                    //fprintf(stderr, "Distance^2 = %3.8f\n", pair.distance2());
+                    //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])));
+                    //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]));
+                    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]));
                 }
                 fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i, j, k, gausGrid_[i][j][k]);
             }
@@ -196,7 +223,65 @@ SANS::finishAnalysis(int nframes)
 void
 SANS::writeOutput()
 {
-
+    FILE *fo;
+    // Construct opendx grid
+    fo = fopen("data.dx","w");
+    //object 1 class gridpositions counts  140 140 140
+    //origin 0.000000 0.000000 0.000000
+    //delta  0.500000 0.000000 0.000000
+    //delta  0.000000 0.500000 0.000000
+    //delta  0.000000 0.000000 0.500000
+    //object 2 class gridconnections counts  140 140 140
+    //object 3 class array type "double" rank 0 items 2744000 data follows
+    fprintf(fo, "object 1 class gridpositions counts %d %d %d\n", gridPoints_[XX],gridPoints_[YY],gridPoints_[ZZ]);
+    fprintf(fo, "origin 0.000000 0.000000 0.000000\n");
+    fprintf(fo, "delta %3.6f %3.6f %3.6f\n", gridSpacing_[XX]*NM2A, 0.0, 0.0);
+    fprintf(fo, "delta %3.6f %3.6f %3.6f\n", 0.0, gridSpacing_[YY]*NM2A, 0.0);
+    fprintf(fo, "delta %3.6f %3.6f %3.6f\n", 0.0, 0.0, gridSpacing_[ZZ]*NM2A);
+    fprintf(fo, "object 2 class gridconnections counts  %d %d %d\n", gridPoints_[XX],gridPoints_[YY],gridPoints_[ZZ]);
+    fprintf(fo, "object 3 class array type \"double\" rank 0 items %d data follows\n", gridPoints_[XX]*gridPoints_[YY]*gridPoints_[ZZ]);
+    // now dump data in ordered mode
+    //    u(0,0,0) u(0,0,1) u(0,0,2)
+    //    ...
+    //    u(0,0,nz-3) u(0,0,nz-2) u(0,0,nz-1)
+    //    u(0,1,0) u(0,1,1) u(0,1,2)
+    //    ...
+    //    u(0,1,nz-3) u(0,1,nz-2) u(0,1,nz-1)
+    //    ...
+    //    u(0,ny-1,nz-3) u(0,ny-1,nz-2) u(0,ny-1,nz-1)
+    //    u(1,0,0) u(1,0,1) u(1,0,2)
+    //    ...
+    int NR = 3;
+    int n = 0;
+    for (int i = 0; i < gridPoints_[XX]; i++)
+    {
+        for (int j = 0; j < gridPoints_[YY]; j++)
+        {
+            for (int k = 0; k < gridPoints_[ZZ]; k++)
+            {
+                fprintf(fo, "%3.8f ", gausGrid_[i][j][k]);
+                n++;
+                if (n == NR) {
+                    fprintf(fo, "\n");
+                    n = 0;
+                }
+            }
+        }
+    }
+    if (n != NR) {
+        fprintf(fo, "\n");
+    }
+    //attribute "dep" string "positions"
+    //object "density" class field
+    //component "positions" value 1
+    //component "connections" value 2
+    //component "data" value 3
+    fprintf(fo, "attribute \"dep\" string \"positions\"\n");
+    fprintf(fo, "object \"density\" class field\n");
+    fprintf(fo, "component \"positions\" value 1\n");
+    fprintf(fo, "component \"connections\" value 2\n");
+    fprintf(fo, "component \"data\" value 3\n");
+    fclose(fo);
 }
 
 /*! \brief