using namespace gmx;
+struct crcl {
+ long double x;
+ long double y;
+ long double z;
+ long double r;
+};
+
+crcl return_crcl(long double x1, long double y1, long double z1,
+ long double x2, long double y2, long double z2,
+ long double x3, long double y3, long double z3) {
+ long double c, b, a, z, y, x, r;
+ 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;
+ a3 = 1 / c;
+ a4 = 1;
+
+ a5 = 2 * (x3 - x1);
+ a6 = 2 * (y3 - y1);
+ a7 = 2 * (z3 - z1);
+ a8 = (x3 * x3 - x1 * x1) + (y3 * y3 - y1 * y1) + (z3 * z3 - z1 * z1);
+
+ a9 = 2 * (x3 - x2);
+ a10 = 2 * (y3 - y2);
+ 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;
+
+ rtrn.x = x;
+ rtrn.y = y;
+ rtrn.z = z;
+ rtrn.r = r;
+
+ return rtrn;
+}
+
/*! \brief
* Template class to serve as a basis for user analysis tools.
*/
-class AnalysisTemplate : public TrajectoryAnalysisModule
+class Spirals : public TrajectoryAnalysisModule
{
public:
- AnalysisTemplate();
+ Spirals();
virtual void initOptions(IOptionsContainer *options,
TrajectoryAnalysisSettings *settings);
Selection selec;
std::vector< std::vector < RVec > > trajectory;
- std::vector< double > noise;
std::vector< int > index;
- double epsi = 0.1; // selectable
- int etalon_frame = 0; // selectable
int frames = 0;
};
-AnalysisTemplate::AnalysisTemplate()
+Spirals::Spirals()
: cutoff_(0.0)
{
registerAnalysisDataset(&data_, "avedist");
}
void
-AnalysisTemplate::initOptions(IOptionsContainer *options,
+Spirals::initOptions(IOptionsContainer *options,
TrajectoryAnalysisSettings *settings)
{
static const char *const desc[] = {
options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
.store(&fnNdx_).defaultBasename("rcore")
.description("Index file from the rcore"));
- // Add option for epsi constant
- options->addOption(DoubleOption("epsilon")
- .store(&epsi)
- .description("thermal vibrations' constant"));
- // Add option for etalon_frame constant
- options->addOption(gmx::IntegerOption("etalon_frame")
- .store(&etalon_frame)
- .description("basic frame to base evaluation on"));
}
+// -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'
+
void
-AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings,
+Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation & /*top*/)
{
index.resize(0);
}
void
-AnalysisTemplate::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{
trajectory.resize(frames + 1);
}
void
-AnalysisTemplate::finishAnalysis(int /*nframes*/)
+Spirals::finishAnalysis(int /*nframes*/)
{
- real *w_rls;
- rvec *x1, *x2, temp;
- std::vector< bool > flag;
- double noise_mid = 9999;
- int count = 0;
- FILE *fpNdx_;
-
- snew(x1, index.size());
- snew(x2, index.size());
- snew(w_rls, index.size());
- flag.resize(index.size(), true);
- noise.resize(index.size(), 0);
-
- while (noise_mid > epsi)
- {
- noise_mid = 0;
- count = 0;
- noise.assign(index.size(), 0);
-
- for (int i = 0; i < index.size(); i++)
- {
- if (flag[i])
- {
- w_rls[i] = 1;
- }
- else
- {
- w_rls[i] = 0;
- }
- }
-
- for (int i = 0; i < frames; i++)
- {
- for (int j = 0; j < trajectory[i].size(); j++)
- {
- copy_rvec(trajectory[i][j].as_vec(), x2[j]);
- copy_rvec(trajectory[etalon_frame][j].as_vec(), x1[j]);
- }
- reset_x(index.size(), NULL, index.size(), NULL, x1, w_rls);
- reset_x(index.size(), NULL, index.size(), NULL, x2, w_rls);
- do_fit(index.size(), w_rls, x1, x2);
- if (i % 100 == 0)
- {
- std::cout << i << " / " << frames << "\n";
- }
- for (int j = 0; j < index.size(); j++)
- {
- rvec_sub(x1[j], x2[j], temp);
- noise[j] += norm(temp);
- }
- }
-
- for (int i = 0; i < index.size(); i++)
- {
- noise[i] /= frames;
- }
-
- for (int i = 0; i < index.size(); i++)
- {
- if (flag[i])
- {
- count++;
- noise_mid += noise[i];
- }
- }
- noise_mid /= count;
-
- for (int i = 0; i < index.size(); i++)
- {
- if (noise[i] > noise_mid)
- {
- flag[i] = false;
+ 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";
}
}
}
- int write_count = 0;
- fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
- std::fprintf(fpNdx_, "[ rcore ]\n");
- for (int i = 0; i < noise.size(); i++)
- {
- if (noise[i] <= epsi)
- {
- write_count++;
- if (write_count > 20) {
- write_count -= 20;
- std::fprintf(fpNdx_, "\n");
- }
- std::fprintf(fpNdx_, "%5d ", index[i] + 1);
+ for (int i = 0; i < 1; i++) {
+ for (int j = 0; j < index.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]);
+ std::cout << test.x << " " << test.y << " " << test.z << " " << test.r << "\n";
}
}
- std::fprintf(fpNdx_,"\n");
- std::fclose(fpNdx_);
- sfree(x1);
- sfree(x2);
- sfree(w_rls);
+
}
void
-AnalysisTemplate::writeOutput()
+Spirals::writeOutput()
{
}
int
main(int argc, char *argv[])
{
- return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<AnalysisTemplate>(argc, argv);
+ return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Spirals>(argc, argv);
}
biod.pnpi.spb.ru / alexxy / gromacs-spirals.git / commitdiff