From 8d444e2d23346d755bc74f5cdff5e62cbc149d36 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Wed, 26 Jul 2017 10:12:22 +0200 Subject: [PATCH] Make loop variables in orires.cpp local Also moved and renamed some variable declarations. Change-Id: I804a1addde81950576ce4a49006ff00ea29db3ab --- src/gromacs/listed-forces/orires.cpp | 250 ++++++++++++++++++----------------- 1 file changed, 127 insertions(+), 123 deletions(-) diff --git a/src/gromacs/listed-forces/orires.cpp b/src/gromacs/listed-forces/orires.cpp index 9996388ee6..4a1f5730e3 100644 --- a/src/gromacs/listed-forces/orires.cpp +++ b/src/gromacs/listed-forces/orires.cpp @@ -38,6 +38,7 @@ #include "orires.h" +#include #include #include "gromacs/gmxlib/network.h" @@ -70,14 +71,6 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, const t_commrec *cr, t_oriresdata *od, t_state *state) { - int i, j, d, ex, nmol, *nr_ex; - double mtot; - rvec com; - gmx_mtop_ilistloop_t iloop; - t_ilist *il; - gmx_mtop_atomloop_all_t aloop; - const gmx_multisim_t *ms; - od->nr = gmx_mtop_ftype_count(mtop, F_ORIRES); if (0 == od->nr) { @@ -102,7 +95,6 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, /* Nothing to do */ return; } - ms = cr->ms; od->fc = ir->orires_fc; od->nex = 0; @@ -111,18 +103,20 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, od->eig = nullptr; od->v = nullptr; - nr_ex = nullptr; - - iloop = gmx_mtop_ilistloop_init(mtop); + int *nr_ex = nullptr; + gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop); + t_ilist *il; + int nmol; while (gmx_mtop_ilistloop_next(iloop, &il, &nmol)) { - for (i = 0; i < il[F_ORIRES].nr; i += 3) + for (int i = 0; i < il[F_ORIRES].nr; i += 3) { - ex = mtop->ffparams.iparams[il[F_ORIRES].iatoms[i]].orires.ex; + int type = il[F_ORIRES].iatoms[i]; + int ex = mtop->ffparams.iparams[type].orires.ex; if (ex >= od->nex) { srenew(nr_ex, ex+1); - for (j = od->nex; j < ex+1; j++) + for (int j = od->nex; j < ex+1; j++) { nr_ex[j] = 0; } @@ -136,6 +130,8 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, * are indentical and the pointers can point to the same memory. */ snew(od->Dinsl, od->nr); + + const gmx_multisim_t *ms = cr->ms; if (ms) { snew(od->Dins, od->nr); @@ -184,17 +180,17 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, } snew(od->tmp, od->nex); snew(od->TMP, od->nex); - for (ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->nex; ex++) { snew(od->TMP[ex], 5); - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { snew(od->TMP[ex][i], 5); } } od->nref = 0; - for (i = 0; i < mtop->natoms; i++) + for (int i = 0; i < mtop->natoms; i++) { if (ggrpnr(&mtop->groups, egcORFIT, i) == 0) { @@ -211,11 +207,12 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, * Copy it to the other nodes after checking multi compatibility, * so we are sure the subsystems match before copying. */ - clear_rvec(com); - mtot = 0.0; - j = 0; - aloop = gmx_mtop_atomloop_all_init(mtop); - const t_atom *atom; + rvec com = { 0, 0, 0 }; + double mtot = 0.0; + int j = 0; + gmx_mtop_atomloop_all_t aloop = gmx_mtop_atomloop_all_init(mtop); + int i = -1; + const t_atom *atom; while (gmx_mtop_atomloop_all_next(aloop, &i, &atom)) { if (mtop->groups.grpnr[egcORFIT] == nullptr || @@ -226,7 +223,7 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, if (ms == nullptr || MASTERSIM(ms)) { copy_rvec(xref[i], od->xref[j]); - for (d = 0; d < DIM; d++) + for (int d = 0; d < DIM; d++) { com[d] += od->mref[j]*xref[i][d]; } @@ -238,14 +235,14 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, svmul(1.0/mtot, com, com); if (ms == nullptr || MASTERSIM(ms)) { - for (j = 0; j < od->nref; j++) + for (int j = 0; j < od->nref; j++) { rvec_dec(od->xref[j], com); } } fprintf(fplog, "Found %d orientation experiments\n", od->nex); - for (i = 0; i < od->nex; i++) + for (int i = 0; i < od->nex; i++) { fprintf(fplog, " experiment %d has %d restraints\n", i+1, nr_ex[i]); } @@ -275,63 +272,63 @@ void init_orires(FILE *fplog, const gmx_mtop_t *mtop, void diagonalize_orires_tensors(t_oriresdata *od) { - int ex, i, j, nrot, ord[DIM], t; - matrix S, TMP; - if (od->M == nullptr) { snew(od->M, DIM); - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { snew(od->M[i], DIM); } snew(od->eig_diag, DIM); snew(od->v, DIM); - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { snew(od->v[i], DIM); } } - for (ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->nex; ex++) { /* Rotate the S tensor back to the reference frame */ + matrix S, TMP; mmul(od->R, od->S[ex], TMP); mtmul(TMP, od->R, S); - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { - for (j = 0; j < DIM; j++) + for (int j = 0; j < DIM; j++) { od->M[i][j] = S[i][j]; } } + int nrot; jacobi(od->M, DIM, od->eig_diag, od->v, &nrot); - for (i = 0; i < DIM; i++) + int ord[DIM]; + for (int i = 0; i < DIM; i++) { ord[i] = i; } - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { - for (j = i+1; j < DIM; j++) + for (int j = i+1; j < DIM; j++) { if (gmx::square(od->eig_diag[ord[j]]) > gmx::square(od->eig_diag[ord[i]])) { - t = ord[i]; + int t = ord[i]; ord[i] = ord[j]; ord[j] = t; } } } - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { od->eig[ex*12 + i] = od->eig_diag[ord[i]]; } - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { - for (j = 0; j < DIM; j++) + for (int j = 0; j < DIM; j++) { od->eig[ex*12 + 3 + 3*i + j] = od->v[j][ord[i]]; } @@ -341,17 +338,16 @@ void diagonalize_orires_tensors(t_oriresdata *od) void print_orires_log(FILE *log, t_oriresdata *od) { - int ex, i; real *eig; diagonalize_orires_tensors(od); - for (ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->nex; ex++) { eig = od->eig + ex*12; fprintf(log, " Orientation experiment %d:\n", ex+1); fprintf(log, " order parameter: %g\n", eig[0]); - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { fprintf(log, " eig: %6.3f %6.3f %6.3f %6.3f\n", (eig[0] != 0) ? eig[i]/eig[0] : eig[i], @@ -368,10 +364,10 @@ real calc_orires_dev(const gmx_multisim_t *ms, const t_mdatoms *md, const rvec x[], const t_pbc *pbc, t_fcdata *fcd, history_t *hist) { - int fa, d, i, j, type, ex, nref; - real edt, edt_1, invn, pfac, r2, invr, corrfac, weight, wsv2, sw, dev; + int nref; + real edt, edt_1, invn, pfac, r2, invr, corrfac, wsv2, sw, dev; tensor *S, R, TMP; - rvec5 *Dinsl, *Dins, *Dtav, *rhs; + rvec5 *rhs; real *mref, ***T; double mtot; rvec *xref, *xtmp, com, r_unrot, r; @@ -391,9 +387,6 @@ real calc_orires_dev(const gmx_multisim_t *ms, edt = od->edt; edt_1 = od->edt_1; S = od->S; - Dinsl = od->Dinsl; - Dins = od->Dins; - Dtav = od->Dtav; T = od->TMP; rhs = od->tmp; nref = od->nref; @@ -425,15 +418,15 @@ real calc_orires_dev(const gmx_multisim_t *ms, } clear_rvec(com); - mtot = 0; - j = 0; - for (i = 0; i < md->nr; i++) + mtot = 0; + int j = 0; + for (int i = 0; i < md->nr; i++) { if (md->cORF[i] == 0) { copy_rvec(x[i], xtmp[j]); mref[j] = md->massT[i]; - for (d = 0; d < DIM; d++) + for (int d = 0; d < DIM; d++) { com[d] += mref[j]*xref[j][d]; } @@ -442,7 +435,7 @@ real calc_orires_dev(const gmx_multisim_t *ms, } } svmul(1.0/mtot, com, com); - for (j = 0; j < nref; j++) + for (int j = 0; j < nref; j++) { rvec_dec(xtmp[j], com); } @@ -450,10 +443,11 @@ real calc_orires_dev(const gmx_multisim_t *ms, calc_fit_R(DIM, nref, mref, xref, xtmp, R); copy_mat(R, od->R); - d = 0; - for (fa = 0; fa < nfa; fa += 3) + /* Index restraint data in order of appearance in forceatoms */ + int restraintIndex = 0; + for (int fa = 0; fa < nfa; fa += 3) { - type = forceatoms[fa]; + int type = forceatoms[fa]; if (pbc) { pbc_dx_aiuc(pbc, x[forceatoms[fa+1]], x[forceatoms[fa+2]], r_unrot); @@ -467,82 +461,88 @@ real calc_orires_dev(const gmx_multisim_t *ms, invr = gmx::invsqrt(r2); /* Calculate the prefactor for the D tensor, this includes the factor 3! */ pfac = ip[type].orires.c*invr*invr*3; - for (i = 0; i < ip[type].orires.power; i++) + for (int i = 0; i < ip[type].orires.power; i++) { pfac *= invr; } - Dinsl[d][0] = pfac*(2*r[0]*r[0] + r[1]*r[1] - r2); - Dinsl[d][1] = pfac*(2*r[0]*r[1]); - Dinsl[d][2] = pfac*(2*r[0]*r[2]); - Dinsl[d][3] = pfac*(2*r[1]*r[1] + r[0]*r[0] - r2); - Dinsl[d][4] = pfac*(2*r[1]*r[2]); + rvec5 &Dinsl = od->Dinsl[restraintIndex]; + Dinsl[0] = pfac*(2*r[0]*r[0] + r[1]*r[1] - r2); + Dinsl[1] = pfac*(2*r[0]*r[1]); + Dinsl[2] = pfac*(2*r[0]*r[2]); + Dinsl[3] = pfac*(2*r[1]*r[1] + r[0]*r[0] - r2); + Dinsl[4] = pfac*(2*r[1]*r[2]); if (ms) { - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { - Dins[d][i] = Dinsl[d][i]*invn; + od->Dins[restraintIndex][i] = Dinsl[i]*invn; } } - d++; + restraintIndex++; } if (ms) { - gmx_sum_sim(5*od->nr, Dins[0], ms); + gmx_sum_sim(5*od->nr, od->Dins[0], ms); } /* Calculate the order tensor S for each experiment via optimization */ - for (ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->nex; ex++) { - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { rhs[ex][i] = 0; - for (j = 0; j <= i; j++) + for (int j = 0; j <= i; j++) { T[ex][i][j] = 0; } } } - d = 0; - for (fa = 0; fa < nfa; fa += 3) + + /* Index restraint data in order of appearance in forceatoms */ + restraintIndex = 0; + for (int fa = 0; fa < nfa; fa += 3) { + rvec5 &Dtav = od->Dtav[restraintIndex]; if (bTAV) { /* Here we update Dtav in t_fcdata using the data in history_t. * Thus the results stay correct when this routine * is called multiple times. */ - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { - Dtav[d][i] = edt*hist->orire_Dtav[d*5+i] + edt_1*Dins[d][i]; + Dtav[i] = + edt*hist->orire_Dtav[restraintIndex*5 + i] + + edt_1*od->Dins[restraintIndex][i]; } } - type = forceatoms[fa]; - ex = ip[type].orires.ex; - weight = ip[type].orires.kfac; + int type = forceatoms[fa]; + int ex = ip[type].orires.ex; + real weight = ip[type].orires.kfac; /* Calculate the vector rhs and half the matrix T for the 5 equations */ - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { - rhs[ex][i] += Dtav[d][i]*ip[type].orires.obs*weight; - for (j = 0; j <= i; j++) + rhs[ex][i] += Dtav[i]*ip[type].orires.obs*weight; + for (int j = 0; j <= i; j++) { - T[ex][i][j] += Dtav[d][i]*Dtav[d][j]*weight; + T[ex][i][j] += Dtav[i]*Dtav[j]*weight; } } - d++; + restraintIndex++; } /* Now we have all the data we can calculate S */ - for (ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->nex; ex++) { /* Correct corrfac and copy one half of T to the other half */ - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { rhs[ex][i] *= corrfac; T[ex][i][i] *= gmx::square(corrfac); - for (j = 0; j < i; j++) + for (int j = 0; j < i; j++) { T[ex][i][j] *= gmx::square(corrfac); T[ex][j][i] = T[ex][i][j]; @@ -555,7 +555,7 @@ real calc_orires_dev(const gmx_multisim_t *ms, S[ex][0][2] = 0; S[ex][1][1] = 0; S[ex][1][2] = 0; - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { S[ex][0][0] += 1.5*T[ex][0][i]*rhs[ex][i]; S[ex][0][1] += 1.5*T[ex][1][i]*rhs[ex][i]; @@ -572,44 +572,49 @@ real calc_orires_dev(const gmx_multisim_t *ms, wsv2 = 0; sw = 0; - d = 0; - for (fa = 0; fa < nfa; fa += 3) + /* Index restraint data in order of appearance in forceatoms */ + restraintIndex = 0; + for (int fa = 0; fa < nfa; fa += 3) { - type = forceatoms[fa]; - ex = ip[type].orires.ex; + int type = forceatoms[fa]; + int ex = ip[type].orires.ex; - od->otav[d] = two_thr* - corrfac*(S[ex][0][0]*Dtav[d][0] + S[ex][0][1]*Dtav[d][1] + - S[ex][0][2]*Dtav[d][2] + S[ex][1][1]*Dtav[d][3] + - S[ex][1][2]*Dtav[d][4]); + const rvec5 &Dtav = od->Dtav[restraintIndex]; + od->otav[restraintIndex] = two_thr* + corrfac*(S[ex][0][0]*Dtav[0] + S[ex][0][1]*Dtav[1] + + S[ex][0][2]*Dtav[2] + S[ex][1][1]*Dtav[3] + + S[ex][1][2]*Dtav[4]); if (bTAV) { - od->oins[d] = two_thr*(S[ex][0][0]*Dins[d][0] + S[ex][0][1]*Dins[d][1] + - S[ex][0][2]*Dins[d][2] + S[ex][1][1]*Dins[d][3] + - S[ex][1][2]*Dins[d][4]); + const rvec5 &Dins = od->Dins[restraintIndex]; + od->oins[restraintIndex] = two_thr* + (S[ex][0][0]*Dins[0] + S[ex][0][1]*Dins[1] + + S[ex][0][2]*Dins[2] + S[ex][1][1]*Dins[3] + + S[ex][1][2]*Dins[4]); } if (ms) { /* When ensemble averaging is used recalculate the local orientation * for output to the energy file. */ - od->oinsl[d] = two_thr* - (S[ex][0][0]*Dinsl[d][0] + S[ex][0][1]*Dinsl[d][1] + - S[ex][0][2]*Dinsl[d][2] + S[ex][1][1]*Dinsl[d][3] + - S[ex][1][2]*Dinsl[d][4]); + const rvec5 &Dinsl = od->Dinsl[restraintIndex]; + od->oinsl[restraintIndex] = two_thr* + (S[ex][0][0]*Dinsl[0] + S[ex][0][1]*Dinsl[1] + + S[ex][0][2]*Dinsl[2] + S[ex][1][1]*Dinsl[3] + + S[ex][1][2]*Dinsl[4]); } - dev = od->otav[d] - ip[type].orires.obs; + dev = od->otav[restraintIndex] - ip[type].orires.obs; wsv2 += ip[type].orires.kfac*gmx::square(dev); sw += ip[type].orires.kfac; - d++; + restraintIndex++; } od->rmsdev = std::sqrt(wsv2/sw); /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */ - for (ex = 0; ex < od->nex; ex++) + for (int ex = 0; ex < od->nex; ex++) { tmmul(R, S[ex], TMP); mmul(TMP, R, S[ex]); @@ -627,8 +632,7 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], const t_mdatoms gmx_unused *md, t_fcdata *fcd, int gmx_unused *global_atom_index) { - int ai, aj; - int fa, d, i, type, ex, power, ki = CENTRAL; + int ex, power, ki = CENTRAL; ivec dt; real r2, invr, invr2, fc, smooth_fc, dev, devins, pfac; rvec r, Sr, fij; @@ -650,12 +654,13 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], smooth_fc *= (1.0 - od->exp_min_t_tau); } - d = 0; - for (fa = 0; fa < nfa; fa += 3) + /* Index restraint data in order of appearance in forceatoms */ + int restraintIndex = 0; + for (int fa = 0; fa < nfa; fa += 3) { - type = forceatoms[fa]; - ai = forceatoms[fa+1]; - aj = forceatoms[fa+2]; + int type = forceatoms[fa]; + int ai = forceatoms[fa + 1]; + int aj = forceatoms[fa + 2]; if (pbc) { ki = pbc_dx_aiuc(pbc, x[ai], x[aj], r); @@ -670,7 +675,7 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], ex = ip[type].orires.ex; power = ip[type].orires.power; fc = smooth_fc*ip[type].orires.kfac; - dev = od->otav[d] - ip[type].orires.obs; + dev = od->otav[restraintIndex] - ip[type].orires.obs; /* NOTE: * there is no real potential when time averaging is applied @@ -680,7 +685,7 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], if (bTAV) { /* Calculate the force as the sqrt of tav times instantaneous */ - devins = od->oins[d] - ip[type].orires.obs; + devins = od->oins[restraintIndex] - ip[type].orires.obs; if (dev*devins <= 0) { dev = 0; @@ -696,12 +701,12 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], } pfac = fc*ip[type].orires.c*invr2; - for (i = 0; i < power; i++) + for (int i = 0; i < power; i++) { pfac *= invr; } mvmul(od->S[ex], r, Sr); - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { fij[i] = -pfac*dev*(4*Sr[i] - 2*(2+power)*invr2*iprod(Sr, r)*r[i]); @@ -713,14 +718,15 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], ki = IVEC2IS(dt); } - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { f[ai][i] += fij[i]; f[aj][i] -= fij[i]; fshift[ki][i] += fij[i]; fshift[CENTRAL][i] -= fij[i]; } - d++; + + restraintIndex++; } } @@ -731,19 +737,17 @@ real orires(int nfa, const t_iatom forceatoms[], const t_iparams ip[], void update_orires_history(t_fcdata *fcd, history_t *hist) { - t_oriresdata *od; - int pair, i; + t_oriresdata *od = &(fcd->orires); - od = &(fcd->orires); if (od->edt != 0) { /* Copy the new time averages that have been calculated * in calc_orires_dev. */ hist->orire_initf = od->exp_min_t_tau; - for (pair = 0; pair < od->nr; pair++) + for (int pair = 0; pair < od->nr; pair++) { - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { hist->orire_Dtav[pair*5+i] = od->Dtav[pair][i]; } -- 2.11.4.GIT