--- /dev/null
+/*
+ * 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 <ericirrgang@gmail.com>
+ * \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 <cstddef>
+
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+namespace gmx
+{
+namespace compat
+{
+
+///\cond
+
+// All gmx::compat code should use std::unique_ptr
+using ::std::unique_ptr;
+
+template<class T> struct Unique_if {
+ typedef unique_ptr<T> Single_object;
+};
+
+template<class T> struct Unique_if<T[]> {
+ typedef unique_ptr<T[]> Unknown_bound;
+};
+
+template<class T, size_t N> struct Unique_if<T[N]> {
+ typedef void Known_bound;
+};
+
+template<class T, class ... Args>
+typename Unique_if<T>::Single_object
+make_unique(Args && ... args)
+{
+ return unique_ptr<T>(new T(::std::forward<Args>(args) ...));
+}
+
+template<class T>
+typename Unique_if<T>::Unknown_bound
+make_unique(size_t n)
+{
+ typedef typename ::std::remove_extent<T>::type U;
+ return unique_ptr<T>(new U[n]());
+}
+
+template<class T, class ... Args>
+typename Unique_if<T>::Known_bound
+make_unique(Args && ...) = delete;
+
+///\endcond
+
+} // namespace compat
+} // namespace gmx
+
+#endif // header guard
* the research papers on the package. Check out http://www.gromacs.org.
*/
#include <iostream>
+#include <cstring>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/selection/nbsearch.h>
#include <gromacs/math/vec.h>
#include <gromacs/math/units.h>
#include <gromacs/topology/atomprop.h>
+#include <gromacs/utility/futil.h>
+#include <gromacs/utility/smalloc.h>
+#include <gromacs/utility/cstringutil.h>
+#include <cstdio>
+#include <unordered_map>
+#include <algorithm>
+#include <gromacs/compat/make_unique.h>
+
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<int>(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<NeutronStructureFactor> readNeutronScatteringFactors()
+{
+ std::vector<NeutronStructureFactor> 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<std::string> 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<std::string> isotopes_;
+
+ //! scattering length of each atom in selection; same order as atomIndex_
+ std::unordered_map<std::string, double> scatterFactors_;
+
+};
+
+std::vector<std::string> getIsotopes(const t_atoms *atoms)
+{
+ std::vector<std::string> 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<std::string> isotopes)
+{
+ std::vector<NeutronStructureFactor> 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.
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<double>>> gausGrid_;
+ std::vector<std::vector<std::vector<double>>> gausGrid_;
std::vector<std::vector<std::vector<bool>>> isSolvent_;
std::vector<double> 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)
{
}
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
gmx_atomprop_t aps = gmx_atomprop_init();
t_atoms *atoms = &(top.topology()->atoms);
real value;
+ std::vector<std::string> isotopes = getIsotopes(atoms);
nb_.setCutoff(cutoff_);
+ nbSol_.setCutoff(cutoffsol_);
top_ = ⊤
+ neutronInfo_ = compat::make_unique<NeutronInfo>(isotopes);
+
for (int i = 0; i < top.topology()->atoms.nr; i++)
{
// Lookup the Van der Waals radius of this atom
}
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);
}
}
// prepare gausGrid
gausGrid_.resize(gridPoints_[XX], std::vector < std::vector < double>>(gridPoints_[YY], std::vector<double>(gridPoints_[ZZ])));
isSolvent_.resize(gridPoints_[XX], std::vector < std::vector < bool>>(gridPoints_[YY], std::vector<bool>(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++)
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, "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
biod.pnpi.spb.ru / alexxy / gromacs-sans.git / commitdiff