#include <gromacs/trajectoryanalysis.h>
#include <gromacs/math/do_fit.h>
#include <gromacs/utility/smalloc.h>
+#include "gromacs/selection/selection.h"
+#include "gromacs/selection/selectionoption.h"
using namespace gmx;
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;
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;
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.
*/
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()
.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'
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*/)
+{
}
biod.pnpi.spb.ru / alexxy / gromacs-spirals.git / blobdiff