using gmx::RVec;
+void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
+{
+ FILE *file;
+ file = std::fopen(file_name, "w+");
+ for (int i = 1; i < trajectory.size(); i++) {
+ std::fprintf(file, "MODEL%9d\n", i);
+ for (int j = 0; j < trajectory[i].size(); j++) {
+ std::fprintf(file, "ATOM %5d CA LYS 1 %8.3f%8.3f%8.3f 1.00 0.00\n", j, trajectory[i][j][0] * 10, trajectory[i][j][1] * 10, trajectory[i][j][2] * 10);
+ // 123456 123456789012345678901 8 6 4567890123456
+ //std::fprintf(file, "1234567890123456789012345678901234567890123456789012345678901234567890\n");
+ }
+ std::fprintf(file, "ENDMDL\n");
+ }
+}
+/*
+1 - 6 Record name "ATOM "
+7 - 11 Integer Atom serial number.
+13 - 16 Atom Atom name.
+17 Character Alternate location indicator.
+18 - 20 Residue name Residue name.
+22 Character Chain identifier.
+23 - 26 Integer Residue sequence number.
+27 AChar Code for insertion of residues.
+31 - 38 Real(8.3) Orthogonal coordinates for X in Angstroms.
+39 - 46 Real(8.3) Orthogonal coordinates for Y in Angstroms.
+47 - 54 Real(8.3) Orthogonal coordinates for Z in Angstroms.
+55 - 60 Real(6.2) Occupancy.
+61 - 66 Real(6.2) Temperature factor (Default = 0.0).
+73 - 76 LString(4) Segment identifier, left-justified.
+77 - 78 LString(2) Element symbol, right-justified.
+79 - 80 LString(2) Charge on the atom.
+*/
+
/*! \brief
* Class used to compute free volume in a simulations box.
*
*
* \ingroup module_trajectoryanalysis
*/
+
+
class Domains : public TrajectoryAnalysisModule
{
public:
SelectionList sel_;
std::vector< std::vector< int > > domains;
+ std::vector< std::vector< int > > domains_local;
std::vector< std::vector< RVec > > trajectory;
std::vector< std::vector< RVec > > frankenstein_trajectory;
options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
.store(&fnNdx_).defaultBasename("domains")
.description("Index file from the domains"));
+ // Add option for selection list
options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
.required().dynamicMask().multiValue()
.description("Domains to form rigid skeleton"));
}
trajectory.resize(1);
- trajectory[0].resize(sel_[0].atomCount());
+ trajectory[0].resize(index.size());
- for (int i = 0; i < sel_[0].atomCount(); i++) {
+ for (int i = 0; i < index.size(); i++) {
trajectory[0][i] = fr.x[index[i]];
}
}
}
- frames++;
-}
-
-void
-Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata)
-{
- std::vector< std::vector< int > > domains_local;
- std::vector< bool > flags;
+ for (int i = 0; i < domains.size(); i++) {
+ for (int j = 0; j < domains[i].size(); j++) {
+ int k = 0;
+ while (index[k] != domains[i][j]) {
+ k++;
+ }
+ domains[i][j] = k;
+ }
+ }
- flags.resize(index.size(), true);
- domains_local = domains; // технически оно за О(1) должно делаться, но не точно
+ for (int i = 0; i < domains.size(); i++) {
+ if (domains[i].size() >= 2) {
+ domains_local.push_back(domains[i]);
+ for (int k = 0; k < domains[i].size(); k++) {
+ for (int j = i + 1; j < domains.size(); j++) {
+ for (int x = 0; x < domains[j].size(); x++) {
+ if (domains[j][x] == domains[i][k]) {
+ domains[j].erase(domains[j].begin() + x);
+ }
+ }
+ }
+ }
+ }
+ }
- trajectory.resize(trajectory.size() + 1);
- trajectory.back().resize(index.size());
+ domains.resize(0);
+ domains = domains_local;
- for (int i = 0; i < index.size(); i++) {
- trajectory.back()[i] = fr.x[index[i]];
- }
+ std::vector< bool > flags;
+ flags.resize(index.size(), true);
for (int i = 0; i < domains.size(); i++) {
for (int j = 0; j < domains[i].size(); j++) {
for (int i = 0; i < index.size(); i++) {
if (flags[i]) {
- int a, b, dist = 9001;
+ int a, b, dist = 90000001;
rvec temp;
- for (int j = 0; j < domains.size(); j++) {
- for (int k = 0; k < domains[j].size(); k++) {
- rvec_sub(trajectory.back()[i], trajectory.back()[domains[j][k]], temp);
- if (norm(temp) <= dist) {
- dist = norm(temp);
- a = j;
- b = k;
+ //if (domains.size() != 2) { спорно как ту поступить
+ for (int j = 0; j < domains.size(); j++) {
+ for (int k = 0; k < domains[j].size(); k++) {
+ rvec_sub(trajectory.back()[i], trajectory.back()[domains[j][k]], temp);
+ if (norm(temp) <= dist) {
+ dist = norm(temp);
+ a = j;
+ b = k;
+ }
}
}
- }
+ //}
domains_local[a].push_back(i);
flags[i] = false;
}
}
- snew(w_rls, domains_local.size());
+ frames++;
+}
+
+void
+Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+ TrajectoryAnalysisModuleData *pdata)
+{
+ trajectory.resize(trajectory.size() + 1);
+ trajectory.back().resize(index.size());
+
+ for (int i = 0; i < index.size(); i++) {
+ trajectory.back()[i] = fr.x[index[i]];
+ }
+
+ snew(w_rls, domains_local.size() + 1);
for (int i = 0; i < domains_local.size(); i++) {
snew(w_rls[i], index.size());
for (int j = 0; j < index.size(); j++) {
w_rls[i][domains_local[i][j]] = 1;
}
}
+ snew(w_rls[domains_local.size()], index.size());
+ for (int i = 0; i < index.size(); i++) {
+ w_rls[domains_local.size()][i] = 1;
+ }
frankenstein_trajectory.resize(trajectory.size());
frankenstein_trajectory.back().resize(index.size());
reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[i]);
do_fit(index.size(), w_rls[i], basic, traj);
+ for (int j = 0; j < index.size(); j++) {
+ if (w_rls[i][j] == 0) {
+ copy_rvec(basic[j], traj[j]);
+ }
+ }
+ reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[domains_local.size()]);
+
for (int j = 0; j < domains_local[i].size(); j++) {
frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
}
frames++;
}
-//domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
+//spacetimecorr -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'
void
Domains::finishAnalysis(int nframes)
{
-
+ make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/frank_with_duo_test.pdb");
}
void
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff