From b0fa8b53f898c551aa0e111e767032216930d9b9 Mon Sep 17 00:00:00 2001 From: Alexey Shvetsov Date: Tue, 18 Sep 2018 10:26:59 +0300 Subject: [PATCH] New update Signed-off-by: Alexey Shvetsov --- include/gromacs/compat/make_unique.h | 118 +++++++++++++ src/sans.cpp | 243 ++++++++++++++++++++++++++- 2 files changed, 353 insertions(+), 8 deletions(-) create mode 100644 include/gromacs/compat/make_unique.h diff --git a/include/gromacs/compat/make_unique.h b/include/gromacs/compat/make_unique.h new file mode 100644 index 0000000..663b15d --- /dev/null +++ b/include/gromacs/compat/make_unique.h @@ -0,0 +1,118 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2017,2018, by the GROMACS development team, led by + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, + * and including many others, as listed in the AUTHORS file in the + * top-level source directory and at http://www.gromacs.org. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +#ifndef GMX_COMPAT_MAKE_UNIQUE_H +#define GMX_COMPAT_MAKE_UNIQUE_H +/*! \libinternal + * \defgroup module_compat C++ standard library compatibility helpers. + * \brief Provide uniform interface to selected C++ standard library features. + * + * For some features not available on all platforms supported by + * \Gromacs, provide back-ports or mappings to available standard + * library implementations as appropriate. + */ + +/*! + * \file + * \brief Provides template gmx::compat::make_unique + * + * The implementation of gmx::compat::make_unique is copied with little + * modification from C++ standardization doc N3656 + * at http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3656.htm + * though additional wrapping has been added for use in \Gromacs. + * + * \author M. Eric Irrgang + * \ingroup module_compat + * \inlibraryapi + */ +/*! \addtogroup module_compat + * ### gmx::compat::make_unique + * + * Provide a trivially adapted implementation of the C++ standard library `make_unique` function template. + * When All supported \Gromacs build platforms provide `std::make_unique`, this should be removed. + * + */ +#include + +#include +#include +#include + +namespace gmx +{ +namespace compat +{ + +///\cond + +// All gmx::compat code should use std::unique_ptr +using ::std::unique_ptr; + +template struct Unique_if { + typedef unique_ptr Single_object; +}; + +template struct Unique_if { + typedef unique_ptr Unknown_bound; +}; + +template struct Unique_if { + typedef void Known_bound; +}; + +template +typename Unique_if::Single_object +make_unique(Args && ... args) +{ + return unique_ptr(new T(::std::forward(args) ...)); +} + +template +typename Unique_if::Unknown_bound +make_unique(size_t n) +{ + typedef typename ::std::remove_extent::type U; + return unique_ptr(new U[n]()); +} + +template +typename Unique_if::Known_bound +make_unique(Args && ...) = delete; + +///\endcond + +} // namespace compat +} // namespace gmx + +#endif // header guard diff --git a/src/sans.cpp b/src/sans.cpp index 3a398e0..7b59339 100644 --- a/src/sans.cpp +++ b/src/sans.cpp @@ -33,15 +33,158 @@ * the research papers on the package. Check out http://www.gromacs.org. */ #include +#include #include #include #include #include #include +#include +#include +#include +#include +#include +#include +#include + using namespace gmx; +gmx_bool get_a_line(FILE *fp, char line[], int n) +{ + char *line0; + char *dum; + + snew(line0, n+1); + + do + { + if (!fgets(line0, n+1, fp)) + { + sfree(line0); + return FALSE; + } + dum = std::strchr(line0, '\n'); + if (dum) + { + dum[0] = '\0'; + } + else if (static_cast(std::strlen(line0)) == n) + { + fprintf(stderr, "Warning: line length exceeds buffer length (%d), data might be corrupted\n", n); + line0[n-1] = '\0'; + } + else + { + fprintf(stderr, "Warning: file does not end with a newline, last line:\n%s\n", + line0); + } + dum = std::strchr(line0, ';'); + if (dum) + { + dum[0] = '\0'; + } + std::strncpy(line, line0, n); + dum = line0; + ltrim(dum); + } + while (dum[0] == '\0'); + + sfree(line0); + return TRUE; +} + + + + +/*! \internal \brief + * Neutron scattering factor parameters for an atom type. + */ +struct NeutronStructureFactor { + + //! atomic isotope + std::string isotope; + //! proton number + int proton; + //! neutron number + int neutron; + //! neutron scattering length + double scatterLength; +}; + + +std::vector readNeutronScatteringFactors() +{ + std::vector neutronScatteringFactors; + FILE *fp; + char line[1000]; + // all non-header lines + fp = libopen("nsfactor.dat"); + while (get_a_line(fp, line, 1000)) + { + int proton; + int neutron; + char currentAtomType[8]; + double scatterLength; + + if (sscanf(line, "%s %d %d %lf", currentAtomType, &proton, &neutron, &scatterLength) == 4) + { + NeutronStructureFactor nsf = {currentAtomType, proton, neutron, scatterLength}; + neutronScatteringFactors.push_back(nsf); + } + } + return neutronScatteringFactors; +}; + +/*! \internal \brief + * Derrived class for computing neutron scattering + */ +class NeutronInfo +{ + public: + //! constructor + NeutronInfo(std::vector isotopes); + + //! retrieves scattering length based on atom index + double getScatterLength(int i); + + private: + //! refId of each atom in selection, should work for dynamic selections + std::vector isotopes_; + + //! scattering length of each atom in selection; same order as atomIndex_ + std::unordered_map scatterFactors_; + +}; + +std::vector getIsotopes(const t_atoms *atoms) +{ + std::vector isotopes; + for (int i = 0; i < atoms->nr; ++i) + { + std::string isotope = atoms->atom[i].elem; + isotopes.push_back(isotope); + } + return isotopes; +}; + +NeutronInfo::NeutronInfo(std::vector isotopes) +{ + std::vector neutronScatterFactors = readNeutronScatteringFactors(); + for (auto scatter : neutronScatterFactors) + { + scatterFactors_[scatter.isotope] = scatter.scatterLength; + } + isotopes_ = std::move(isotopes); +}; + +double NeutronInfo::getScatterLength(int i) +{ + std::string isotope = isotopes_[i]; + double scattering = scatterFactors_[isotope]; + return scattering; +}; /*! \brief * Template class to serve as a basis for user analysis tools. @@ -71,21 +214,30 @@ class SANS : public TrajectoryAnalysisModule double Gaussian(double value, double sigma, RVec Rpos, RVec Rgrid, RVec GridSpacing); double cutoff_; + double cutoffsol_; double grid_; double boxScale_; + double d2oRatio_; IVec gridPoints_; RVec gridSpacing_; AnalysisNeighborhood nb_; + AnalysisNeighborhood nbSol_; const TopologyInformation *top_; - std::vector< std::vector < std::vector>> gausGrid_; + std::vector>> gausGrid_; std::vector>> isSolvent_; std::vector vdw_radius_; + friend class NeutronInfo; + std::unique_ptr < NeutronInfo> neutronInfo_; }; SANS::SANS() - : cutoff_(1.0), grid_(0.05), boxScale_(1.0) + : cutoff_(1.0), + cutoffsol_(0.5), + grid_(0.05), + boxScale_(1.0), + d2oRatio_(1.0) { } @@ -111,12 +263,18 @@ SANS::initOptions(IOptionsContainer *options, options->addOption(DoubleOption("cutoff") .store(&cutoff_) .description("cutoff for neighbour search")); + options->addOption(DoubleOption("cutoffsol") + .store(&cutoffsol_) + .description("cutoff for pure solvent 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")); + options->addOption(DoubleOption("d2oratio") + .store(&d2oRatio_) + .description("D2O/H2O ration")); } void @@ -127,10 +285,14 @@ SANS::initAnalysis(const TrajectoryAnalysisSettings &settings, gmx_atomprop_t aps = gmx_atomprop_init(); t_atoms *atoms = &(top.topology()->atoms); real value; + std::vector isotopes = getIsotopes(atoms); nb_.setCutoff(cutoff_); + nbSol_.setCutoff(cutoffsol_); top_ = ⊤ + neutronInfo_ = compat::make_unique(isotopes); + for (int i = 0; i < top.topology()->atoms.nr; i++) { // Lookup the Van der Waals radius of this atom @@ -144,10 +306,10 @@ SANS::initAnalysis(const TrajectoryAnalysisSettings &settings, } else { - fprintf(stderr, "Could not determine VDW radius for %s-%s. Set to zero.\n", + fprintf(stderr, "Could not determine VDW radius for %s-%s. Set to 0.15.\n", *(atoms->resinfo[resnr].name), *(atoms->atomname[i])); - vdw_radius_.push_back(0.0); + vdw_radius_.push_back(0.15); } } @@ -167,16 +329,30 @@ SANS::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, // prepare gausGrid gausGrid_.resize(gridPoints_[XX], std::vector < std::vector < double>>(gridPoints_[YY], std::vector(gridPoints_[ZZ]))); isSolvent_.resize(gridPoints_[XX], std::vector < std::vector < bool>>(gridPoints_[YY], std::vector(gridPoints_[ZZ]))); + // set isSolvent default to true... + // 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++) + { + isSolvent_[i][j][k] = true; + } + } + } } void SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata) { - RVec GridSpacing(0.1, 0.1, 0.1); AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, AnalysisNeighborhoodPositions(fr.x, fr.natoms)); + AnalysisNeighborhoodSearch nbsearchSol = nbSol_.initSearch(pbc, AnalysisNeighborhoodPositions(fr.x, fr.natoms)); t_topology *top = top_->topology(); + int nsol = 0; + double solDens = 0.0 ; fprintf(stderr, "Box dimensions:\n"); for (int i = 0; i < DIM; i++) @@ -196,6 +372,7 @@ SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, 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)) @@ -211,16 +388,66 @@ SANS::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, //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])); - // Mark if this is a pure solvent point e.g. no non-solvent atoms within cut-off + //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,"NeutronInfo %6.6f\n",neutronInfo_->getScatterLength(pair.refIndex())); + if (!strcmp(*(top_->topology()->atoms.resinfo[top_->topology()->atoms.atom[pair.refIndex()].resind].name), "SOL") + & !(strcmp(top_->topology()->atoms.atom[pair.refIndex()].elem,"H"))) + { + gausGrid_[i][j][k] += Gaussian((1.-d2oRatio_)*(-3.7406) + d2oRatio_*6.6710, vdw_radius_[pair.refIndex()], refPos, point, gridSpacing_)/(gridSpacing_[XX]*gridSpacing_[YY]*gridSpacing_[ZZ]); + } else { + gausGrid_[i][j][k] += Gaussian(neutronInfo_->getScatterLength(pair.refIndex()), vdw_radius_[pair.refIndex()], refPos, point, gridSpacing_)/(gridSpacing_[XX]*gridSpacing_[YY]*gridSpacing_[ZZ]); + } + } + AnalysisNeighborhoodPairSearch pairSearchSol = nbsearchSol.startPairSearch(point.as_vec()); + AnalysisNeighborhoodPair pairSol; + while (pairSearchSol.findNextPair(&pairSol)) { + // Mark if this is a pure solvent point e.g. no non-solvent atoms within cut-off + //fprintf(stderr,"atom %d resname %s\n", pair.refIndex(), *(top_->topology()->atoms.resinfo[top_->topology()->atoms.atom[pair.refIndex()].resind].name)); + //fprintf(stderr,"Is sol %d\n", strcmp(*(top_->topology()->atoms.resinfo[top_->topology()->atoms.atom[pair.refIndex()].resind].name), "SOL")); + isSolvent_[i][j][k] = isSolvent_[i][j][k]&(!strcmp(*(top_->topology()->atoms.resinfo[top_->topology()->atoms.atom[pairSol.refIndex()].resind].name), "SOL")); + //fprintf(stderr, "bool = %s\n", isSolvent_[i][j][k]?"true":"false"); + //fprintf(stderr, "strncmp = %s", !strcmp(*(top_->topology()->atoms.resinfo[top_->topology()->atoms.atom[pair.refIndex()].resind].name), "SOL")?"true": "false" ); + } + //fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i, j, k, gausGrid_[i][j][k]); + //fprintf(stderr, "isSolvent[%d, %d, %d] = %s\n", i, j, k, isSolvent_[i][j][k]?"true":"false"); + if (isSolvent_[i][j][k]) { + nsol++; } - fprintf(stderr, "GausGrid[%d, %d, %d] = %f\n", i, j, k, gausGrid_[i][j][k]); } } } // only for orthogonal boxes... + for (int i = 0; i < gridPoints_[XX]; i++) + { + for (int j = 0; j < gridPoints_[YY]; j++) + { + for (int k = 0; k < gridPoints_[ZZ]; k++) + { + if (isSolvent_[i][j][k]) { + solDens += gausGrid_[i][j][k]/float(nsol); + } + } + } + } + + // + // + for (int i = 0; i < gridPoints_[XX]; i++) + { + for (int j = 0; j < gridPoints_[YY]; j++) + { + for (int k = 0; k < gridPoints_[ZZ]; k++) + { + if (isSolvent_[i][j][k]) { + gausGrid_[i][j][k] -= solDens; + } + } + } + } fprintf(stderr, "NX = %d, NY = %d, NZ = %d\n", gridPoints_[XX], gridPoints_[YY], gridPoints_[ZZ]); + fprintf(stderr, "nsol = %d\n", nsol); + fprintf(stderr, "solDens = %6.6f\n", solDens); } void -- 2.22.0