From 4bda0baf70b921fcb0655f672488ea5bf03ac191 Mon Sep 17 00:00:00 2001 From: Alexey Shvetsov Date: Mon, 10 Sep 2018 09:14:24 +0300 Subject: [PATCH] Normalization to volume, vdw radii as sigma and opendx output Signed-off-by: Alexey Shvetsov --- src/sans.cpp | 99 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 92 insertions(+), 7 deletions(-) diff --git a/src/sans.cpp b/src/sans.cpp index 4c932d9..5c13d80 100644 --- a/src/sans.cpp +++ b/src/sans.cpp @@ -38,6 +38,8 @@ #include #include #include +#include + 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>> gausGrid_; + std::vector 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_ = ⊤ + gmx_atomprop_t aps = gmx_atomprop_init(); + t_atoms *atoms = &(top.topology()->atoms); + real value; + for(int i=0; iatoms.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 -- 2.22.0