From dbc91f6f0c817ca7d72bfea39fb472c7bf5d1dd3 Mon Sep 17 00:00:00 2001 From: Anatoly Titov Date: Tue, 12 Dec 2017 12:29:07 +0300 Subject: [PATCH] added evaluations for each frame --- src/spirals.cpp | 194 ++++++++++++++++++++++++------------------------ 1 file changed, 99 insertions(+), 95 deletions(-) diff --git a/src/spirals.cpp b/src/spirals.cpp index 62f1adb..15debf0 100644 --- a/src/spirals.cpp +++ b/src/spirals.cpp @@ -38,6 +38,8 @@ #include #include #include +#include "gromacs/selection/selection.h" +#include "gromacs/selection/selectionoption.h" using namespace gmx; @@ -55,12 +57,8 @@ crcl return_crcl(long double x1, long double y1, long double z1, long double a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12; c = (z3 - z1 * x3 / x1 - (z2 - z1 * x2 / x1) / (y2 - y1 * x2 / x1)) / (1 - x3 / x1 - (1 - x2 / x1) / (y2 - y1 * x2 / x1)); - //std::cout << c << " "; - //std::cout << y1 * (x2 / x1) << " "; b = (y2 - y1 * x2 / x1) / (1 - x2 / x1 - 1 / c * (z2 - z1 * x2 / x1)); - //std::cout << b << " "; a = (x1) / (1 - y1 / b - z1 / c); - //std::cout << a << " "; a1 = 1 / a; a2 = 1 / b; @@ -77,19 +75,11 @@ crcl return_crcl(long double x1, long double y1, long double z1, a11 = 2 * (z3 - z2); a12 = (x3 * x3 - x2 * x2) + (y3 * y3 - y2 * y2) + (z3 * z3 - z2 * z2); - /*std::cout << a1 << " " << a2 << " " << a3 << " " << a4 << " " - << a5 << " " << a6 << " " << a7 << " " << a8 << " " - << a9 << " " << a10 << " " << a11 << " " << a12 << " ";*/ - z = (a12 - a9 * a4 / a1 - (a10 - a9 * a2 / a1) * (a8 - a5 * a4 / a1) / (a6 - a5 * a2 / a1)) / (a11 - a9 * a3 / a1 - (a10 - a9 * a2 / a1) * (a7 - a5 * a3 / a1) / (a6 - a5 * a2 / a1)); - //std::cout << z << " "; y = 1 / (a6 - a5 * a2 / a1) * (a8 - a5 * a4 / a1 - z * (a7 - a5 * a3 / a1)); - //std::cout << y << " "; x = 1 / a1 * (a4 - a2 * y - a3 * z); - //std::cout << x << " "; r = std::sqrt((x - x3) * (x - x3) + (y - y3) * (y - y3) + (z - z3) * (z - z3)); - //std::cout << r << " "; crcl rtrn; @@ -101,6 +91,15 @@ crcl return_crcl(long double x1, long double y1, long double z1, return rtrn; } +long double lr_turn(std::vector< long double > a, std::vector< long double > b, std::vector< long double > c) { + return a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + b[0] * c[1] * a[2] - + (a[2] * b[1] * c[0] + a[0] * b[2] * c[1] + a[1] * b[0] * c[2]); +} + +long double kernel_dist(crcl a, crcl b) { + return std::sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z)); +} + /*! \brief * Template class to serve as a basis for user analysis tools. */ @@ -132,13 +131,11 @@ class Spirals : public TrajectoryAnalysisModule AnalysisData data_; AnalysisDataAverageModulePointer avem_; - - Selection selec; SelectionList sel_; - std::vector< std::vector < RVec > > trajectory; - std::vector< std::vector< int > > index_all; - std::vector< int > index; int frames = 0; + std::vector< std::vector< crcl > > kernel; + std::vector< std::vector< std::vector< int > > > circles; + std::vector< std::vector< long double > > spiral_dist; }; Spirals::Spirals() @@ -162,9 +159,9 @@ Spirals::initOptions(IOptionsContainer *options, .store(&fnNdx_).defaultBasename("rcore") .description("Index file from the rcore")); // Add option for selection list - options->addOption(SelectionOption("Select groups").storeVector(&sel_) + options->addOption(SelectionOption("select").storeVector(&sel_) .required().dynamicMask().multiValue() - .description("Groups for masses")); + .description("Position to calculate distances for")); } // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test.ndx' -on '/home/toluk/Develop/samples/reca_rd/core.ndx' @@ -173,98 +170,105 @@ void Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation & /*top*/) { - std::vector< int > a; - index_all.resize(0); - index.resize(0); - for (int i = 0; i < sel_.size(); i++) { - a.resize(0); - ConstArrayRef< int > atomind = sel_[i].atomIndices(); - for (ConstArrayRef< int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) { - a.push_back(*ai); - index.push_back(*ai); - } - index_all.push_back(a); - std::cout << "one selection parsed\n"; - } - std::cout << "First part finished\n"; + } void Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata) { - trajectory.resize(frames + 1); - trajectory[frames].resize(index.size()); - for (int i = 0; i < index.size(); i++) - { - trajectory[frames][i] = fr.x[index[i]]; + kernel.resize(frames + 1); + circles.resize(frames + 1); + spiral_dist.resize(frames + 1); + + const SelectionList &sel = pdata->parallelSelections(sel_); + std::vector< RVec > temp; + temp.resize(sel_.size()); + for (int i = 0; i < sel.size(); i++) { + copy_rvec(sel[i].position(0).x(), temp[i]); } - frames++; -} -void -Spirals::finishAnalysis(int /*nframes*/) -{ - std::cout << "trajectory finished\n"; - long double x1, y1, z1, x2, y2, z2, x3, y3, z3; - crcl test; - - for (long double r = 1; r < 50; r++) { - for (long double x0 = -5; x0 < 6; x0+=0.99) { - for (long double y0 = -5; y0 < 6; y0+=0.99) { - x1 = x0 + r - r/100; - x2 = x0 + r/125; - x3 = x0 - r/7; - z1 = 1.01; - z2 = 1; - z3 = 0.99; - y1 = std::sqrt(r * r - (x1 - x0) * (x1 - x0)) + y0; - y2 = std::sqrt(r * r - (x2 - x0) * (x2 - x0)) + y0; - y3 = std::sqrt(r * r - (x3 - x0) * (x3 - x0)) + y0; - - /*std::cout << x1 << " " << y1 << " " << z1 << " " - << x2 << " " << y2 << " " << z2 << " " - << x3 << " " << y3 << " " << z3 << "\n";*/ - - test = return_crcl(x1, y1, z1, x2, y2, z2, x3, y3, z3); - - //std::cout.precision(1); - //std::cout.width(3); - //std::cout << std::abs(x0 - test.x) << " " << std::abs(y0 - test.y) << " " << std::abs(1 - test.z) << " " << std::abs(r - test.r) << "\n"; + // центры окружностей и их радиусы + for (int i = 0; i < temp.size() - 2; i++) { + kernel[frames].push_back(return_crcl( temp[i][0], temp[i][1], temp[i][2], + temp[i + 1][0], temp[i + 1][1], temp[i + 1][2], + temp[i + 2][0], temp[i + 2][1], temp[i + 2][2])); + } + + // распределение точек по виткам + std::vector< long double > a, b, c; + a.resize(3); + b.resize(3); + c.resize(3); + + a[0] = temp[0][0] - kernel[frames].begin()->x; + a[1] = temp[0][1] - kernel[frames].begin()->y; + a[2] = temp[0][2] - kernel[frames].begin()->z; + + c[0] = kernel[frames].end()->x - kernel[frames].begin()->x; + c[1] = kernel[frames].end()->y - kernel[frames].begin()->y; + c[2] = kernel[frames].end()->z - kernel[frames].begin()->z; + + std::vector< int > empty; + empty.resize(0); + + + circles[frames].push_back(empty); + circles[frames][0].push_back(0); + + long double sign = 0, prev_sign = 0; + for (int i = 1; i < temp.size(); i++) { + if (i < kernel[frames].size()) { + b[0] = temp[i][0] - kernel[frames][i].x; + b[1] = temp[i][1] - kernel[frames][i].y; + b[2] = temp[i][2] - kernel[frames][i].z; + } else { + b[0] = temp[i][0] - kernel[frames].end()->x; + b[1] = temp[i][1] - kernel[frames].end()->y; + b[2] = temp[i][2] - kernel[frames].end()->z; + } + long double local_sign = lr_turn(a, b ,c); + if (sign == 0) { + sign = local_sign; + prev_sign = local_sign; + circles[frames].back().push_back(i); + } else { + if (std::signbit(local_sign) == std::signbit(prev_sign)) { + circles[frames].back().push_back(i); + } else { + if (std::signbit(local_sign) != std::signbit(sign)) { + circles[frames].back().push_back(i); + } else { + circles[frames].push_back(empty); + circles[frames].back().push_back(i); + } } + prev_sign = local_sign; } } - std::vector< std::vector< long double > > masses; - masses.resize(0); - int count = 0; - for (int i = 0; i < index_all.size(); i++) { - std::vector< long double > center; - center.resize(3, 0); - for (int j = 0; j < index_all[i].size(); j++) { - center[0] += trajectory[0][count][0]; - center[1] += trajectory[0][count][1]; - center[2] += trajectory[0][count][2]; - count++; + // шаг спирали + spiral_dist[frames].resize(circles[frames].size(), 0); + for (int i = 0; i < circles[frames].size(); i++) { + if (i == 0) { + spiral_dist[frames][i] = kernel_dist(kernel[frames][circles[frames][i].front()], kernel[frames][circles[frames][i].back()]) + + kernel_dist(kernel[frames][circles[frames][i].back()], kernel[frames][circles[frames][i + 1].front()]) / 2; + } else if (i != circles[frames].size() - 1) { + spiral_dist[frames][i] = kernel_dist(kernel[frames][circles[frames][i - 1].back()], kernel[frames][circles[frames][i].front()]) / 2 + + kernel_dist(kernel[frames][circles[frames][i].front()], kernel[frames][circles[frames][i].back()]) + + kernel_dist(kernel[frames][circles[frames][i].back()], kernel[frames][circles[frames][i + 1].front()]) / 2; + } else if (i == circles[frames].size() - 1) { + spiral_dist[frames][i] = kernel_dist(kernel[frames][circles[frames][i - 1].back()], kernel[frames][circles[frames][i].front()]) / 2 + + kernel_dist(kernel[frames][circles[frames][i].front()], kernel[frames][circles[frames][i].back()]); } - center[0] /= index_all[i].size(); - center[1] /= index_all[i].size(); - center[2] /= index_all[i].size(); - masses.push_back(center); } + frames++; +} - for (int i = 0; i < 1; i++) { - for (int j = 0; j < masses.size() - 2; j++) { - /*test = return_crcl( trajectory[i][j][0], trajectory[i][j][1], trajectory[i][j][2], - trajectory[i][j + 1][0], trajectory[i][j + 1][1], trajectory[i][j + 1][2], - trajectory[i][j + 2][0], trajectory[i][j + 2][1], trajectory[i][j + 2][2]);*/ - test = return_crcl( masses[j][0], masses[j][1], masses[j][2], - masses[j + 1][0], masses[j + 1][1], masses[j + 1][2], - masses[j + 2][0], masses[j + 2][1], masses[j + 2][2]); - std::cout << test.x << " " << test.y << " " << test.z << " " << test.r << "\n"; - } - } +void +Spirals::finishAnalysis(int /*nframes*/) +{ } -- 2.22.0