#include <cstring>
#include <algorithm>
+ #include <unordered_map>
#include "gromacs/commandline/pargs.h"
#include "gromacs/commandline/viewit.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/gmxana/gstat.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/listed-forces/disre.h"
+#include "gromacs/listed_forces/disre.h"
#include "gromacs/math/do_fit.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/mdatoms.h"
-#include "gromacs/mdlib/mdrun.h"
#include "gromacs/mdtypes/fcdata.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/topology/index.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/trajectoryanalysis/topologyinformation.h"
#include "gromacs/utility/arraysize.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
ndr, label, i, n);
}
- if (ndr != label)
- {
- gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
- }
do
{
n += nat;
}
typedef struct {
- int index;
+ int label;
gmx_bool bCore;
real up1, r, rT3, rT6, viol, violT3, violT6;
} t_dr_stats;
break;
}
fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
- drs[i].index, yesno_names[drs[i].bCore],
+ drs[i].label, yesno_names[drs[i].bCore],
drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
drs[i].viol, drs[i].violT3, drs[i].violT6);
}
return bIC;
}
- static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
+ static void dump_stats(FILE *log, int nsteps,
+ const t_disresdata &dd,
+ const t_ilist *disres,
t_iparams ip[], t_dr_result *dr,
int isize, int index[], t_atoms *atoms)
{
- int i, j, nra;
t_dr_stats *drs;
fprintf(log, "\n");
fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
- snew(drs, ndr);
- j = 0;
- nra = interaction_function[F_DISRES].nratoms+1;
- for (i = 0; (i < ndr); i++)
+ snew(drs, dd.nres);
+ const int nra = interaction_function[F_DISRES].nratoms + 1;
+ for (int j = 0; j < disres->nr; j += nra)
{
- /* Search for the appropriate restraint */
- for (; (j < disres->nr); j += nra)
- {
- if (ip[disres->iatoms[j]].disres.label == i)
- {
- break;
- }
- }
- drs[i].index = i;
- drs[i].bCore = is_core(i, isize, index);
+ // Note that the restraint i can be used by multiple pairs
+ const int i = disres->iatoms[j] - dd.type_min;
+ GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range");
+
+ drs[i].label = ip[disres->iatoms[j]].disres.label;
+ drs[i].bCore = is_core(drs[i].label, isize, index);
drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
drs[i].r = dr->aver1[i]/nsteps;
drs[i].rT3 = gmx::invcbrt(dr->aver_3[i]/nsteps);
atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
}
}
- dump_viol(log, ndr, drs, FALSE);
+ dump_viol(log, dd.nres, drs, FALSE);
fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
- std::sort(drs, drs+ndr, [](const t_dr_stats &a, const t_dr_stats &b)
+ std::sort(drs, drs + dd.nres, [](const t_dr_stats &a, const t_dr_stats &b)
{return a.viol > b.viol; }); //Reverse sort
- dump_viol(log, ndr, drs, TRUE);
+ dump_viol(log, dd.nres, drs, TRUE);
- dump_dump(log, ndr, drs);
+ dump_dump(log, dd.nres, drs);
sfree(drs);
}
- static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
+ static void dump_clust_stats(FILE *fp,
+ const t_disresdata &dd,
+ const t_ilist *disres,
t_iparams ip[], t_blocka *clust, t_dr_result dr[],
char *clust_name[], int isize, int index[])
{
- int i, j, k, nra, mmm = 0;
+ int k, nra, mmm = 0;
double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
t_dr_stats *drs;
fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
- snew(drs, ndr);
+ snew(drs, dd.nres);
for (k = 0; (k < clust->nr); k++)
{
{
gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
}
- j = 0;
nra = interaction_function[F_DISRES].nratoms+1;
sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
- for (i = 0; (i < ndr); i++)
+
+ // Use a map to process each restraint only once while looping over all pairs
+ std::unordered_map<int, bool> restraintHasBeenProcessed;
+ for (int j = 0; j < dd.nres; j += nra)
{
- /* Search for the appropriate restraint */
- for (; (j < disres->nr); j += nra)
+ // Note that the restraint i can be used by multiple pairs
+ const int i = disres->iatoms[j] - dd.type_min;
+
+ if (restraintHasBeenProcessed[i])
{
- if (ip[disres->iatoms[j]].disres.label == i)
- {
- break;
- }
+ continue;
}
- drs[i].index = i;
- drs[i].bCore = is_core(i, isize, index);
+
+ drs[i].label = ip[disres->iatoms[j]].disres.label;
+ drs[i].bCore = is_core(drs[i].label, isize, index);
drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
drs[i].r = dr[k].aver1[i]/dr[k].nframes;
if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
maxV = std::max(maxV, static_cast<double>(drs[i].viol));
maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
+
+ // We have processed restraint i, mark it as such
+ restraintHasBeenProcessed[i] = true;
}
if (std::strcmp(clust_name[k], "1000") == 0)
{
"because this analysis is done from a trajectory file afterwards",
"it does not make sense to use time averaging. However,",
"the time averaged values per restraint are given in the log file.[PAR]",
- "An index file may be used to select specific restraints for",
+ "An index file may be used to select specific restraints by index group label for",
"printing.[PAR]",
"When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
"amount of average violations.[PAR]",
};
FILE *out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
- t_tpxheader header;
- gmx_mtop_t mtop;
- rvec *xtop;
- gmx_localtop_t *top;
- t_atoms *atoms = nullptr;
+ gmx_localtop_t top;
t_fcdata fcd;
- t_nrnb nrnb;
t_graph *g;
- int ntopatoms, natoms, i, j, kkk;
+ int i, j, kkk;
t_trxstatus *status;
real t;
rvec *x, *xav = nullptr;
init5(ntop);
}
- t_inputrec irInstance;
- t_inputrec *ir = &irInstance;
+ t_inputrec irInstance;
+ t_inputrec *ir = &irInstance;
- read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
- snew(xtop, header.natoms);
- read_tpx(ftp2fn(efTPR, NFILE, fnm), ir, box, &ntopatoms, xtop, nullptr, &mtop);
+ gmx::TopologyInformation topInfo;
+ topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
+ int ntopatoms = topInfo.mtop()->natoms;
+ AtomsDataPtr atoms;
bPDB = opt2bSet("-q", NFILE, fnm);
if (bPDB)
{
ind_fit[kkk] = kkk;
}
- snew(atoms, 1);
- *atoms = gmx_mtop_global_atoms(&mtop);
+ atoms = topInfo.copyAtoms();
if (atoms->pdbinfo == nullptr)
{
atoms->havePdbInfo = TRUE;
}
- top = gmx_mtop_generate_local_top(&mtop, ir->efep != efepNO);
+ gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
g = nullptr;
pbc_null = nullptr;
}
else
{
- g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
+ g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE);
}
}
}
ir->dr_tau = 0.0;
- init_disres(fplog, &mtop, ir, nullptr, nullptr, &fcd, nullptr, FALSE);
+ init_disres(fplog, topInfo.mtop(), ir, nullptr, nullptr, &fcd, nullptr, FALSE);
- natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
+ int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
snew(f, 5*natoms);
init_dr_res(&dr, fcd.disres.nres);
"Largest Violation", "Time (ps)", "nm", oenv);
}
- auto mdAtoms = gmx::makeMDAtoms(fplog, mtop, *ir, false);
- atoms2md(&mtop, ir, -1, nullptr, mtop.natoms, mdAtoms.get());
+ auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
+ atoms2md(topInfo.mtop(), ir, -1, nullptr, ntopatoms, mdAtoms.get());
update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
- init_nrnb(&nrnb);
if (ir->ePBC != epbcNONE)
{
- gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
+ gpbc = gmx_rmpbc_init(&top.idef, ir->ePBC, natoms);
}
j = 0;
}
my_clust = clust->inv_clust[j];
range_check(my_clust, 0, clust->clust->nr);
- check_viol(fplog, &(top->idef.il[F_DISRES]),
- top->idef.iparams,
+ check_viol(fplog, &(top.idef.il[F_DISRES]),
+ top.idef.iparams,
x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
}
else
{
- check_viol(fplog, &(top->idef.il[F_DISRES]),
- top->idef.iparams,
+ check_viol(fplog, &(top.idef.il[F_DISRES]),
+ top.idef.iparams,
x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
}
if (bPDB)
if (clust)
{
- dump_clust_stats(fplog, fcd.disres.nres, &(top.idef.il[F_DISRES]),
- dump_clust_stats(fplog, fcd.disres, &(top->idef.il[F_DISRES]),
- top->idef.iparams, clust->clust, dr_clust,
++ dump_clust_stats(fplog, fcd.disres, &(top.idef.il[F_DISRES]),
+ top.idef.iparams, clust->clust, dr_clust,
clust->grpname, isize, index);
}
else
{
- dump_stats(fplog, j, fcd.disres.nres, &(top.idef.il[F_DISRES]),
- dump_stats(fplog, j, fcd.disres, &(top->idef.il[F_DISRES]),
- top->idef.iparams, &dr, isize, index,
- bPDB ? atoms : nullptr);
++ dump_stats(fplog, j, fcd.disres, &(top.idef.il[F_DISRES]),
+ top.idef.iparams, &dr, isize, index,
+ bPDB ? atoms.get() : nullptr);
if (bPDB)
{
write_sto_conf(opt2fn("-q", NFILE, fnm),
"Coloured by average violation in Angstrom",
- atoms, xav, nullptr, ir->ePBC, box);
+ atoms.get(), xav, nullptr, ir->ePBC, box);
}
dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
- j, &top->idef, &mtop, max_dr, nlevels, bThird);
+ j, &top.idef, topInfo.mtop(), max_dr, nlevels, bThird);
xvgrclose(out);
xvgrclose(aver);
xvgrclose(numv);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff