crl[i][j].resize(traj.front().size(), 0);
}
}
- rvec temp_zero;
+ RVec temp_zero;
temp_zero[0] = 0;
temp_zero[1] = 0;
temp_zero[2] = 0;
std::vector< float > d;
d.resize(traj.front().size(), 0);
- #pragma omp parallel shared(crl, temp_zero, d)
- {
- #pragma omp for schedule(dynamic)
+ //#pragma omp parallel shared(crl, temp_zero, d)
+ //{
+ //#pragma omp for schedule(dynamic)
for (int i = tauS; i <= tauE; i += 1) {
- rvec *medx, *medy, temp1, temp2;
- real tmp;
- snew(medx, traj.front().size());
- snew(medy, traj.front().size());
- for (int i = 0; i < traj.front().size(); i++) {
- copy_rvec(temp_zero, medx[i]);
- copy_rvec(temp_zero, medy[i]);
- }
+ std::vector< RVec > medx, medy;
+ RVec temp1, temp2;
+ float tmp;
+ medx.resize(traj.front().size(), temp_zero);
+ medy.resize(traj.front().size(), temp_zero);
for (int j = 0; j < traj.size() - i - 1; j++) {
for (int f = 0; f < traj.front().size(); f++) {
- rvec_sub(traj[b_frame][f], traj[j][f], temp1);
- rvec_inc(medx[f], temp1);
-
- rvec_sub(traj[b_frame][f], traj[j + i][f], temp1);
- rvec_inc(medy[f], temp1);
+ medx[f] += (traj[b_frame][f] - traj[j][f]);
+ medy[f] += (traj[b_frame][f] - traj[j + i][f]);
}
}
for (int j = 0; j < traj.front().size(); j++) {
tmp -= i;
tmp = 1 / tmp;
//tmp = 1 / (traj.size() - 1 - i);
- copy_rvec(medx[j], temp1);
- svmul(tmp, temp1, medx[j]);
- copy_rvec(medy[j], temp1);
- svmul(tmp, temp1, medy[j]);
+ medx[j] *= tmp;
+ medy[j] *= tmp;
}
std::vector< std::vector< float > > a, b, c;
a.resize(traj.front().size(), d);
for (int j = 0; j < traj.size() - i - 1; j++) {
for (int f1 = 0; f1 < traj.front().size(); f1++) {
for (int f2 = 0; f2 < traj.front().size(); f2++) {
- rvec_sub(traj[b_frame][f1], traj[j][f1], temp1);
- rvec_dec(temp1, medx[f1]);
-
- rvec_sub(traj[b_frame][f2], traj[j + i][f2], temp2);
- rvec_dec(temp2, medy[f2]);
+ temp1 = traj[b_frame][f1] - traj[j][f1] - medx[f1];
+ temp2 = traj[b_frame][f2] - traj[j + i][f2] - medy[f2];
a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
}
}
- sfree(medx);
- sfree(medy);
+ medx.resize(0);
+ medy.resize(0);
std::cout << i << " corr done\n";
}
- }
- #pragma omp barrier
+ //}
+ //#pragma omp barrier
}
void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
for (int i = 0; i < traj.front().size(); i++) {
graph[i].resize(traj.front().size(), std::make_pair(0, -1));
}
- rvec temp;
+ RVec temp;
for (int i = 1; i <= tauE; i++) {
for (int j = 0; j < traj.front().size(); j++) {
for (int f = j; f < traj.front().size(); f++) {
- copy_rvec(traj[b_frame][j], temp);
- rvec_dec(temp, traj[b_frame][f]);
+ temp = traj[b_frame][j] - traj[b_frame][f];
if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b && norm(temp) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
graph[j][f].first = crl[i][j][f];
* \ingroup module_trajectoryanalysis
*/
-class Domains : public TrajectoryAnalysisModule
+class SpaceTimeCorr : public TrajectoryAnalysisModule
{
public:
- Domains();
- virtual ~Domains();
+ SpaceTimeCorr();
+ virtual ~SpaceTimeCorr();
//! Set the options and setting
virtual void initOptions(IOptionsContainer *options,
private:
- std::string fnNdx_;
SelectionList sel_;
-
- std::vector< std::vector< int > > domains_local;
std::vector< std::vector< RVec > > trajectory;
- std::vector< std::vector< RVec > > frankenstein_trajectory;
std::vector< int > index;
- std::vector< int > numbers;
int frames = 0;
int basic_frame = 0;
int tau = 0;
float crl_border = 0;
float eff_rad = 1.5;
- real **w_rls;
- std::vector< std::vector< std::pair< int, int > > > w_rls2;
// Copy and assign disallowed by base.
};
-Domains::Domains(): TrajectoryAnalysisModule()
+SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
{
}
-Domains::~Domains()
+SpaceTimeCorr::~SpaceTimeCorr()
{
}
void
-Domains::initOptions(IOptionsContainer *options,
+SpaceTimeCorr::initOptions(IOptionsContainer *options,
TrajectoryAnalysisSettings *settings)
{
static const char *const desc[] = {
// Add the descriptive text (program help text) to the options
settings->setHelpText(desc);
// Add option for output file name
- options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
- .store(&fnNdx_).defaultBasename("domains")
- .description("Index file from the domains"));
+ //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
+ // .store(&fnNdx_).defaultBasename("domains")
+ // .description("Index file from the domains"));
// Add option for tau constant
options->addOption(gmx::IntegerOption("tau")
.store(&tau)
}
void
-Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
+SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
-{
-}
-
-void
-Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr)
{
ArrayRef< const int > atomind = sel_[0].atomIndices();
index.resize(0);
for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
index.push_back(*ai);
}
- trajectory.resize(1);
- trajectory.back().resize(index.size());
+ trajectory.resize(0);
+}
- for (int i = 0; i < index.size(); i++) {
- trajectory.back()[i] = fr.x[index[i]];
- }
- frames++;
+void
+SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
+ const t_trxframe &fr)
+{
}
// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5' -tau 5 -crl 0.10
// -s '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest1.pdb' -f '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest.pdb' -n '/home/toluk/Develop/samples/CubeTermal/testCa.ndx' -sf '/home/toluk/Develop/samples/CubeTermal/SelListCa' -tau 900 -crl 0.20 -ef_rad 1
void
-Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{
trajectory.resize(trajectory.size() + 1);
}
void
-Domains::finishAnalysis(int nframes)
+SpaceTimeCorr::finishAnalysis(int nframes)
{
std::vector< std::vector< std::vector< float > > > crltns;
std::vector< std::vector< std::pair< float, int > > > graph;
std::cout << "\nCorrelation's evaluation - start\n";
- correlation_evaluation(frankenstein_trajectory, basic_frame, crltns, m, k);
+ correlation_evaluation(trajectory, basic_frame, crltns, m, k);
std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
- graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, k);
+ graph_calculation(graph, sub_graph, sub_graph_rbr, trajectory, basic_frame, crltns, crl_border, eff_rad, k);
std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
std::cout << "routs evaluation: end\n";
- make_correlation_file(crltns, "CubeT.txt", 0);
- make_rout_file2(crl_border, index, rout_new, "Routs_Cube.txt");
- make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_Cube.txt");
+ //make_correlation_file(crltns, "CubeT.txt", 0);
+ //make_rout_file2(crl_border, index, rout_new, "Routs_Cube.txt");
+ //make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_Cube.txt");
std::cout << "Finish Analysis - end\n\n";
}
void
-Domains::writeOutput()
+SpaceTimeCorr::writeOutput()
{
}
int
main(int argc, char *argv[])
{
- return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
+ return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);
}
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff