diff --git a/README.md b/README.md index dae8a2c..216bd28 100644 --- a/README.md +++ b/README.md @@ -3,40 +3,16 @@ [![Build Status](https://travis-ci.org/votca/csgapps.svg?branch=master)](https://travis-ci.org/votca/csgapps) [![pipeline status](https://gitlab.com/votca/csgapps/badges/master/pipeline.svg)](https://gitlab.com/votca/csgapps/commits/master) -This repository contains some small analysis programs which are written based on the votca_csg framework +This repository contains some small analysis programs which are written based on the votca csg framework # Programs -* __radii__ - - calculate gyration radius + hydrodynamic radius of a molecule or a set of molecules - -* __sphericalorder__ - - calculates a spherical order parameter i.e. the distribution of e_r\*u where e_r is the unit vector from - a reference molecule to the solvent molecules and u is the principle axis of inertia of the solvent molecule - -* __fluctuations__ - - calculate the number density fluctuations in subvolumes. Subvolumes are either cubic slabs along one of the - sim box axis or speherical subvolumes around a solute - -* __orientcorr__ - - calculates the distance dependent orientational correlation function of a polymer melt - -* __part_dist__ - - outputs the time-averaged number of particles, listed by particle types (was a part of csg before) - -* __partial_rdf__ - - calculates the rdf in a spherical subvolume - -* __traj_force__ - - add/subtracts reference forces from a given trajectory and stores in a new trajectory - - +* _radii_: calculate gyration radius + hydrodynamic radius of a molecule or a set of molecules +* _sphericalorder_: calculates a spherical order parameter i.e. the distribution of e\_r\*u where e\_r is the unit vector from a reference molecule to the solvent molecules and u is the principle axis of inertia of the solvent molecule +* _fluctuations_: calculate the number density fluctuations in subvolumes. Subvolumes are either cubic slabs along one of the sim box axis or speherical subvolumes around a solute +* _orientcorr_: calculates the distance dependent orientational correlation function of a polymer melt +* _part_dist_: outputs the time-averaged number of particles, listed by particle types (was a part of csg before) +* _partial\_rdf_: calculates the rdf in a spherical subvolume +* _traj\_force_: add/subtracts reference forces from a given trajectory and stores in a new trajectory To add your own program just create a new folder and put your *.cc files in there diff --git a/fluctuations/fluctuations.cc b/fluctuations/fluctuations.cc index 3273513..88ffbaa 100644 --- a/fluctuations/fluctuations.cc +++ b/fluctuations/fluctuations.cc @@ -16,16 +16,11 @@ */ #include -#include -#include -#include -#include #include #include -#include +#include #include -// using namespace votca::tools; using namespace std; using namespace votca::csg; @@ -77,7 +72,7 @@ class CsgFluctuations : public CsgApplication { bool DoTrajectory() override { return true; } bool DoMapping() override { return true; } - void BeginEvaluate(Topology *top, Topology *top_atom) override { + void BeginEvaluate(Topology *top, Topology *) override { _filter = OptionsMap()["filter"].as(); _refmol = OptionsMap()["refmol"].as(); _rmin = OptionsMap()["rmin"].as(); @@ -102,17 +97,11 @@ class CsgFluctuations : public CsgApplication { cout << "Doing slabs along z-axis" << endl; _dim = 2; } else { - cout << "Unrecognized geometry option. (sphere|x|y|z)" << endl; - exit(0); + throw std::runtime_error("Unrecognized geometry option. (sphere|x|y|z)"); } - _N_avg = new double[_nbins]; - _N_sq_avg = new double[_nbins]; - N = new int[_nbins]; - for (int i = 0; i < _nbins; i++) { - _N_avg[i] = 0; - _N_sq_avg[i] = 0; - } + _N_avg = Eigen::VectorXd::Zero(_nbins); + _N_sq_avg = Eigen::VectorXd::Zero(_nbins); if (_do_spherical) { cout << "Calculating fluctions for " << _rmin << "getBox(); - Eigen::Vector3d a = box.col(0); - Eigen::Vector3d b = box.col(1); - Eigen::Vector3d c = box.col(2); - _ref = (a + b + c) / 2; + _ref = box.rowwise().sum() / 2; - cout << "Refernce is center of box " << _ref << endl; + cout << "Reference is center of box " << _ref << endl; } - _outfile.open(_outfilename.c_str()); + _outfile.open(_outfilename); if (!_outfile) { - throw runtime_error("cannot open outfile for output"); + throw runtime_error("cannot open" + _outfilename + " for output"); } } @@ -147,10 +133,9 @@ class CsgFluctuations : public CsgApplication { protected: // number of particles in dV int _nbins; - double *_N_avg; + Eigen::VectorXd _N_avg; // sqare - double *_N_sq_avg; - int *N; + Eigen::VectorXd _N_sq_avg; string _filter; string _refmol; double _rmax; @@ -170,11 +155,7 @@ int main(int argc, char **argv) { return app.Exec(argc, argv); } -void CsgFluctuations::EvalConfiguration(Topology *conf, - Topology *conf_atom = nullptr) { - Eigen::Vector3d eR; - double r = 0; - int rbin; +void CsgFluctuations::EvalConfiguration(Topology *conf, Topology *) { if (_refmol != "") { for (Bead *bead : conf->Beads()) { @@ -185,21 +166,20 @@ void CsgFluctuations::EvalConfiguration(Topology *conf, } } - for (int i = 0; i < _nbins; i++) { - N[i] = 0; - } + votca::tools::HistogramNew hist; + hist.Initialize(_rmin, _rmax, _nbins); /* check how many molecules are in each bin*/ for (Bead *bead : conf->Beads()) { if (!votca::tools::wildcmp(_filter.c_str(), bead->getName().c_str())) { continue; } - + double r = 0; if (_do_spherical) { - eR = bead->getPos() - _ref; + Eigen::Vector3d eR = bead->getPos() - _ref; r = eR.norm(); } else { - eR = bead->getPos(); + Eigen::Vector3d eR = bead->getPos(); if (_dim == 0) { r = eR.x(); } else if (_dim == 1) { @@ -208,17 +188,12 @@ void CsgFluctuations::EvalConfiguration(Topology *conf, r = eR.z(); } } - if (r > _rmin && r < _rmax) { - rbin = (int)_nbins * (double)((r - _rmin) / (_rmax - _rmin)); - N[rbin]++; - } + hist.Process(r); } /* update averages*/ - for (int i = 0; i < _nbins; i++) { - _N_avg[i] += N[i]; - _N_sq_avg[i] += N[i] * N[i]; - } + _N_avg += hist.data().y(); + _N_sq_avg += hist.data().y().cwiseAbs2(); _nframes++; } diff --git a/orientcorr/orientcorr.cc b/orientcorr/orientcorr.cc index f38ec8f..df9796a 100644 --- a/orientcorr/orientcorr.cc +++ b/orientcorr/orientcorr.cc @@ -80,7 +80,6 @@ string OrientCorrApp::_nbmethod; class MyWorker : public CsgApplication::Worker { public: ~MyWorker() override = default; - ; // evaluate the current frame void EvalConfiguration(Topology *top, Topology *top_ref) override; @@ -135,7 +134,7 @@ NBList *OrientCorrApp::CreateNBSearch() { } // initialize the histograms -void OrientCorrApp::BeginEvaluate(Topology *top, Topology *top_ref) { +void OrientCorrApp::BeginEvaluate(Topology *, Topology *) { _cor.Initialize(0, _cut_off, _nbins); _count.Initialize(0, _cut_off, _nbins); _cor_excl.Initialize(0, _cut_off, _nbins); @@ -155,7 +154,7 @@ CsgApplication::Worker *OrientCorrApp::ForkWorker() { } // evaluates a frame -void MyWorker::EvalConfiguration(Topology *top, Topology *top_ref) { +void MyWorker::EvalConfiguration(Topology *top, Topology *) { // first genearate a mapped topology // the beads are sitting on the bonds and have an orientation which @@ -193,8 +192,8 @@ void MyWorker::EvalConfiguration(Topology *top, Topology *top_ref) { cout << "done\n"; // the neighbor search only finds pairs, add self-self correlation parts here - _cor.Process(0.0f, mapped.BeadCount()); - _count.Process(0.0f, mapped.BeadCount()); + _cor.Process(0.0f, (double)mapped.BeadCount()); + _count.Process(0.0f, (double)mapped.BeadCount()); // search for all beads BeadList b; @@ -215,7 +214,7 @@ void MyWorker::EvalConfiguration(Topology *top, Topology *top_ref) { // process a pair, since return value is falsed, pairs are not cached which // saves a lot of memory for the big systems -bool MyWorker::FoundPair(Bead *b1, Bead *b2, const Eigen::Vector3d &r, +bool MyWorker::FoundPair(Bead *b1, Bead *b2, const Eigen::Vector3d &, const double dist) { double tmp = b1->getV().dot(b2->getV()); double P2 = 3. / 2. * tmp * tmp - 0.5; diff --git a/partial_rdf/rdf_calculator.cc b/partial_rdf/rdf_calculator.cc index 430d11d..4d7e74e 100644 --- a/partial_rdf/rdf_calculator.cc +++ b/partial_rdf/rdf_calculator.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include namespace votca { @@ -31,6 +32,7 @@ namespace csg { RDFCalculator::RDFCalculator() : _write_every(0), _do_blocks(false), + _nblock(0), _subvol_rad(0), _do_vol_corr(false), _processed_some_frames(false) {} @@ -56,14 +58,11 @@ void RDFCalculator::Initialize() { interaction_t *i = AddInteraction(prop); i->_is_bonded = false; } -}; +} -void RDFCalculator::BeginEvaluate(Topology *top, Topology *top_atom) { +void RDFCalculator::BeginEvaluate(Topology *top, Topology *) { Eigen::Matrix3d box = top->getBox(); - Eigen::Vector3d a = box.col(0); - Eigen::Vector3d b = box.col(1); - Eigen::Vector3d c = box.col(2); - _boxc = a / 2 + b / 2 + c / 2; + _boxc = box.rowwise().sum() / 2.0; std::cout << "Using center of box: " << _boxc << std::endl; // we didn't process any frames so far @@ -99,21 +98,12 @@ void RDFCalculator::BeginEvaluate(Topology *top, Topology *top_atom) { } // calculate normalization factor for rdf - /*if ((*iter)->get("type1").value() == (*iter)->get("type2").value()) - i._norm = 1. / (4. * M_PI * i._step * beads1.size()*(beads2.size() - 1.) - / 2.); else i._norm = 1. / (4. * M_PI * i._step * beads1.size() * - beads2.size());*/ - /*if ((*iter)->get("type1").value() == (*iter)->get("type2").value()) - i._norm = 1. / ( (beads2.size() - 1.) / 2.); - else - i._norm = 1. / ( beads2.size());*/ - if (_do_vol_corr) { std::cout << "Volume correction on" << std::endl; - i._norm = 1. / (4.0 * M_PI * i._step); + i._norm = 1. / (4.0 * votca::tools::conv::Pi * i._step); } else { std::cout << "Volume correction off" << std::endl; - i._norm = 1. / (4. * M_PI * i._step); + i._norm = 1. / (4. * votca::tools::conv::Pi * i._step); } } } @@ -167,9 +157,8 @@ void RDFCalculator::LoadOptions(const std::string &file) { } // evaluate current conformation -void RDFCalculator::Worker::EvalConfiguration(Topology *top, - Topology *top_atom) { - _cur_vol = 4.0 / 3.0 * M_PI * _rdfcalculator->_subvol_rad * +void RDFCalculator::Worker::EvalConfiguration(Topology *top, Topology *) { + _cur_vol = 4.0 / 3.0 * votca::tools::conv::Pi * _rdfcalculator->_subvol_rad * _rdfcalculator->_subvol_rad * _rdfcalculator->_subvol_rad; // process non-bonded interactions DoNonbonded(top); @@ -178,18 +167,14 @@ void RDFCalculator::Worker::EvalConfiguration(Topology *top, } void RDFCalculator::ClearAverages() { - std::map::iterator ic_iter; - std::map::iterator group_iter; _nframes = 0; - for (ic_iter = _interactions.begin(); ic_iter != _interactions.end(); - ++ic_iter) { - ic_iter->second->_average.Clear(); + for (auto &_interaction : _interactions) { + _interaction.second->_average.Clear(); } - for (group_iter = _groups.begin(); group_iter != _groups.end(); - ++group_iter) { - group_iter->second->_corr.setZero(); + for (auto &_group : _groups) { + _group.second->_corr.setZero(); } } @@ -207,8 +192,7 @@ class IMCNBSearchHandler { Eigen::Vector3d _boxc; // center of box bool _do_vol_corr; - bool FoundPair(Bead *b1, Bead *b2, const Eigen::Vector3d &r, - const double dist) { + bool FoundPair(Bead *b1, Bead *, const Eigen::Vector3d &, const double dist) { if (_do_vol_corr) { double dr = (b1->Pos() - _boxc).norm(); @@ -250,8 +234,8 @@ void RDFCalculator::Worker::DoNonbonded(Topology *top) { _rdfcalculator->_boxc, _rdfcalculator->_subvol_rad); - _cur_beadlist_1_count = beads1.size(); - _cur_beadlist_2_count = beads2.size(); + _cur_beadlist_1_count = (double)beads1.size(); + _cur_beadlist_2_count = (double)beads2.size(); // same types, so put factor 1/2 because of already counted interactions if (prop->get("type1").value() == prop->get("type2").value()) { @@ -259,7 +243,7 @@ void RDFCalculator::Worker::DoNonbonded(Topology *top) { } // generate the neighbour list - NBList *nb; + std::unique_ptr nb; bool gridsearch = true; @@ -275,9 +259,9 @@ void RDFCalculator::Worker::DoNonbonded(Topology *top) { } } if (gridsearch) { - nb = new NBListGrid(); + nb = std::make_unique(NBListGrid()); } else { - nb = new NBList(); + nb = std::make_unique(NBList()); } nb->setCutoff(i._max + i._step); @@ -300,14 +284,6 @@ void RDFCalculator::Worker::DoNonbonded(Topology *top) { // store particle number in subvolume for each interaction i._avg_beadlist_1_count.Process(_cur_beadlist_1_count); i._avg_beadlist_2_count.Process(_cur_beadlist_2_count); - - // process all pairs - /*NBList::iterator pair_iter; - for(pair_iter = nb->begin(); pair_iter!=nb->end();++pair_iter) { - _current_hists[i._index].Process((*pair_iter)->dist()); - }*/ - - delete nb; } } @@ -324,9 +300,7 @@ void RDFCalculator::Worker::DoBonded(Topology *top) { // now fill with new data std::list list = top->InteractionsInGroup(name); - std::list::iterator ic_iter; - for (ic_iter = list.begin(); ic_iter != list.end(); ++ic_iter) { - Interaction *ic = *ic_iter; + for (auto ic : list) { double v = ic->EvaluateVar(*top); _current_hists[i._index].Process(v); } @@ -345,27 +319,29 @@ RDFCalculator::group_t *RDFCalculator::getGroup(const std::string &name) { // write the distribution function void RDFCalculator::WriteDist(const std::string &suffix) { - std::map::iterator iter; // for all interactions - for (iter = _interactions.begin(); iter != _interactions.end(); ++iter) { + for (auto &_interaction : _interactions) { // calculate the rdf - Table &t = iter->second->_average.data(); + Table &t = _interaction.second->_average.data(); Table dist(t); - iter->second->_norm /= (iter->second->_avg_beadlist_1_count.getAvg() * - iter->second->_avg_beadlist_2_count.getAvg()); - dist.y() = _avg_vol.getAvg() * iter->second->_norm * + _interaction.second->_norm /= + (_interaction.second->_avg_beadlist_1_count.getAvg() * + _interaction.second->_avg_beadlist_2_count.getAvg()); + dist.y() = _avg_vol.getAvg() * _interaction.second->_norm * dist.y().cwiseQuotient(dist.x().cwiseAbs2()); - dist.Save((iter->first) + suffix + ".dist.new"); - std::cout << "written " << (iter->first) + suffix + ".dist.new\n"; + dist.Save((_interaction.first) + suffix + ".dist.new"); + std::cout << "written " << (_interaction.first) + suffix + ".dist.new\n"; - std::cout << "Avg. number of particles in subvol for " << (iter->first) - << std::endl; - std::cout << "beadlist 1: " << iter->second->_avg_beadlist_1_count.getAvg() + std::cout << "Avg. number of particles in subvol for " + << (_interaction.first) << std::endl; + std::cout << "beadlist 1: " + << _interaction.second->_avg_beadlist_1_count.getAvg() << std::endl; - std::cout << "beadlist 2: " << iter->second->_avg_beadlist_2_count.getAvg() + std::cout << "beadlist 2: " + << _interaction.second->_avg_beadlist_2_count.getAvg() << std::endl; } @@ -376,14 +352,12 @@ void RDFCalculator::WriteDist(const std::string &suffix) { CsgApplication::Worker *RDFCalculator::ForkWorker() { RDFCalculator::Worker *worker; worker = new RDFCalculator::Worker; - std::map::iterator ic_iter; worker->_current_hists.resize(_interactions.size()); worker->_rdfcalculator = this; - for (ic_iter = _interactions.begin(); ic_iter != _interactions.end(); - ++ic_iter) { - interaction_t *i = ic_iter->second; + for (auto &_interaction : _interactions) { + interaction_t *i = _interaction.second; worker->_current_hists[i->_index].Initialize( i->_average.getMin(), i->_average.getMax(), i->_average.getNBins()); } @@ -395,16 +369,13 @@ void RDFCalculator::MergeWorker(CsgApplication::Worker *worker_) { RDFCalculator::Worker *worker = dynamic_cast(worker_); // update the average - std::map::iterator ic_iter; - // map::iterator group_iter; ++_nframes; _avg_vol.Process(worker->_cur_vol); - for (ic_iter = _interactions.begin(); ic_iter != _interactions.end(); - ++ic_iter) { - interaction_t *i = ic_iter->second; + for (auto &_interaction : _interactions) { + interaction_t *i = _interaction.second; i->_average.data().y() = (((double)_nframes - 1.0) * i->_average.data().y() + worker->_current_hists[i->_index].data().y()) / diff --git a/partial_rdf/rdf_calculator.h b/partial_rdf/rdf_calculator.h index 19e7c3c..46e5e3e 100644 --- a/partial_rdf/rdf_calculator.h +++ b/partial_rdf/rdf_calculator.h @@ -79,7 +79,7 @@ class RDFCalculator { /// struct to store collected information for interactions struct interaction_t { - int _index; + long int _index; Property *_p; HistogramNew _average; double _min, _max, _step; diff --git a/radii/radii.cc b/radii/radii.cc index 62a2524..5ffc795 100644 --- a/radii/radii.cc +++ b/radii/radii.cc @@ -56,7 +56,7 @@ int main(int argc, char **argv) { return app.Exec(argc, argv); } -void CsgTestApp::EvalConfiguration(Topology *top, Topology *top_ref) { +void CsgTestApp::EvalConfiguration(Topology *top, Topology *) { // loop over all molecules for (Molecule *mol : top->Molecules()) { // does the id match if given? @@ -73,7 +73,7 @@ void CsgTestApp::EvalConfiguration(Topology *top, Topology *top_ref) { } // Number of beads in the molecule - int N = mol->BeadCount(); + long int N = mol->BeadCount(); // sqared tensor of gyration for current snapshot double r_gyr_sq = 0; diff --git a/sphericalorder/sphericalorder.cc b/sphericalorder/sphericalorder.cc index 5196a94..9fdda18 100644 --- a/sphericalorder/sphericalorder.cc +++ b/sphericalorder/sphericalorder.cc @@ -71,7 +71,7 @@ class CGOrderParam : public CsgApplication { bool DoTrajectory() override { return true; } bool DoMapping() override { return true; } - void BeginEvaluate(Topology *top, Topology *top_atom) override { + void BeginEvaluate(Topology *top, Topology *) override { string filter; @@ -107,19 +107,17 @@ class CGOrderParam : public CsgApplication { Eigen::Matrix3d box = top->getBox(); Eigen::Vector3d a = box.col(0); - Eigen::Vector3d b = box.col(1); - Eigen::Vector3d c = box.col(2); if (_refmol == "") { - _ref = (a + b + c) / 2; + _ref = box.rowwise().sum() / 2; cout << "Refernce is center of box " << _ref << endl; } boxl = a.norm() / 2; if (_rbinw > 0) { - _rbins = boxl / _rbinw + 1; + _rbins = (int)(boxl / _rbinw) + 1; cout << "radial bins " << _rbins << endl; } else { _rbins = 1; @@ -141,9 +139,6 @@ class CGOrderParam : public CsgApplication { for (int i = 0; i < _rbins; i++) { _nmol[i] = 0; } - - // cout << "Test" << endl; - // _hist_u[1][10] =0.0; } void EndEvaluate() override { @@ -178,8 +173,7 @@ class CGOrderParam : public CsgApplication { _file_w.close(); }; - void EvalConfiguration(Topology *conf, - Topology *conf_atom = nullptr) override { + void EvalConfiguration(Topology *conf, Topology * = nullptr) override { Eigen::Vector3d eR; int nu, nv, nw; diff --git a/traj_force/traj_force.cc b/traj_force/traj_force.cc index c150576..3c6ae44 100644 --- a/traj_force/traj_force.cc +++ b/traj_force/traj_force.cc @@ -52,7 +52,7 @@ bool TrajForce::EvaluateOptions() { return true; } -void TrajForce::BeginEvaluate(Topology *top, Topology *top_atom) { +void TrajForce::BeginEvaluate(Topology *top, Topology *) { _top_force.CopyTopologyData(top); _trjreader_force = TrjReaderFactory().Create(_op_vm["trj-force"].as()); @@ -85,14 +85,13 @@ void TrajForce::EndEvaluate() { void TrajForce::WriteOutFiles() {} -void TrajForce::EvalConfiguration(Topology *conf, Topology *conf_atom) { +void TrajForce::EvalConfiguration(Topology *conf, Topology *) { if (conf->BeadCount() != _top_force.BeadCount()) { throw std::runtime_error( "number of beads in topology and reference force topology does not " "match"); } for (int i = 0; i < conf->BeadCount(); ++i) { - // conf->getBead(i)->F() += _scale*_top_force.getBead(i)->getF(); // \todo check why "conf" HasForce() is false // Since "conf" topology Force is set to false @@ -108,8 +107,7 @@ void TrajForce::EvalConfiguration(Topology *conf, Topology *conf_atom) { "trajectory differ by more than 1e-6"); } } - // cout << conf->HasForce() << endl; - // cout << _top_force.HasForce() << endl; + _trjwriter->Write(&_top_force); _trjreader_force->NextFrame(_top_force); }