From: Berk Hess Date: Tue, 2 Dec 2014 11:16:59 +0000 (+0100) Subject: Add value_ref and value to pull_coord_t X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=4cf66de4b630b29329c9699091c4fa78d96bdf92;p=alexxy%2Fgromacs.git Add value_ref and value to pull_coord_t This is only code refactoring. Change-Id: I64b4cb05f69325a1dd1e17f269894cb538663817 --- diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c index c92eada292..68752443c2 100644 --- a/src/gromacs/gmxpreprocess/readpull.c +++ b/src/gromacs/gmxpreprocess/readpull.c @@ -427,13 +427,9 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l { t_mdatoms *md; t_pull *pull; - t_pull_coord *pcrd; - t_pull_group *pgrp0, *pgrp1; t_pbc pbc; - int c, m; + int c; double t_start; - real init = 0; - double dev, value; init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0); md = init_mdatoms(NULL, mtop, ir->efep); @@ -453,6 +449,11 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n"); for (c = 0; c < pull->ncoord; c++) { + t_pull_coord *pcrd; + t_pull_group *pgrp0, *pgrp1; + double value; + real init = 0; + pcrd = &pull->coord[c]; pgrp0 = &pull->group[pcrd->group[0]]; @@ -468,9 +469,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l pcrd->init = 0; } - get_pull_coord_distance(pull, c, &pbc, t_start, &dev); - /* Calculate the value of the coordinate at t_start */ - value = pcrd->init + t_start*pcrd->rate + dev; + get_pull_coord_value(pull, c, &pbc, &value); fprintf(stderr, " %10.3f nm", value); if (pcrd->bStart) diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h index 01b5451b7b..95cfc402cd 100644 --- a/src/gromacs/legacyheaders/types/inputrec.h +++ b/src/gromacs/legacyheaders/types/inputrec.h @@ -139,6 +139,8 @@ typedef struct { real kB; /* force constant for state B */ /* Variables not present in mdp, but used at run time */ + double value_ref; /* The reference value, usually init+rate*t */ + double value; /* The current value of the coordinate */ dvec dr; /* The distance from the reference group */ double vec_len; /* Length of vec for direction-relative */ dvec ffrad; /* conversion factor from vec to radial force */ diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c index d7cf619295..a252388ef5 100644 --- a/src/gromacs/pulling/pull.c +++ b/src/gromacs/pulling/pull.c @@ -74,46 +74,20 @@ static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp) } static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd, - gmx_bool bPrintRefValue, double t, + gmx_bool bPrintRefValue, gmx_bool bPrintComponents) { - double distance, distance2; - int m; - - if (pcrd->eGeom == epullgDIST) - { - /* Geometry=distance: print the length of the distance vector */ - distance2 = 0; - for (m = 0; m < DIM; m++) - { - if (pcrd->dim[m]) - { - distance2 += pcrd->dr[m]*pcrd->dr[m]; - } - } - distance = sqrt(distance2); - } - else - { - /* Geometry=direction-like: print distance along the pull vector */ - distance = 0; - for (m = 0; m < DIM; m++) - { - if (pcrd->dim[m]) - { - distance += pcrd->dr[m]*pcrd->vec[m]; - } - } - } - fprintf(out, "\t%g", distance); + fprintf(out, "\t%g", pcrd->value); if (bPrintRefValue) { - fprintf(out, "\t%g", pcrd->init + pcrd->rate*t); + fprintf(out, "\t%g", pcrd->value_ref); } if (bPrintComponents) { + int m; + for (m = 0; m < DIM; m++) { if (pcrd->dim[m]) @@ -151,7 +125,7 @@ static void pull_print_x(FILE *out, t_pull *pull, double t) pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]); } pull_print_coord_dr(out, pcrd, - pull->bPrintRefValue, t, pull->bPrintComp); + pull->bPrintRefValue, pull->bPrintComp); } fprintf(out, "\n"); } @@ -469,7 +443,7 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc) */ static void low_get_pull_coord_dr(const t_pull *pull, const t_pull_coord *pcrd, - const t_pbc *pbc, double t, + const t_pbc *pbc, dvec xg, dvec xref, double max_dist2, dvec dr) { @@ -496,7 +470,7 @@ static void low_get_pull_coord_dr(const t_pull *pull, { for (m = 0; m < DIM; m++) { - dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m]; + dref[m] = pcrd->value_ref*pcrd->vec[m]; } /* Add the reference position, so we use the correct periodic image */ dvec_inc(xrefr, dref); @@ -524,9 +498,9 @@ static void low_get_pull_coord_dr(const t_pull *pull, /* This function returns the distance based on the contents of the pull struct. * pull->coord[coord_ind].dr, and potentially vector, are updated. */ -static void get_pull_coord_dr(const t_pull *pull, - int coord_ind, - const t_pbc *pbc, double t) +static void get_pull_coord_dr(t_pull *pull, + int coord_ind, + const t_pbc *pbc) { double md2; t_pull_coord *pcrd; @@ -573,51 +547,95 @@ static void get_pull_coord_dr(const t_pull *pull, } } - low_get_pull_coord_dr(pull, pcrd, pbc, t, + low_get_pull_coord_dr(pull, pcrd, pbc, pull->group[pcrd->group[1]].x, pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x, md2, pcrd->dr); } -void get_pull_coord_distance(const t_pull *pull, - int coord_ind, - const t_pbc *pbc, double t, - double *dev) +static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t) { - static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety, - but is fairly benign */ - const t_pull_coord *pcrd; - int m; - double ref, drs, inpr; + /* With zero rate the reference value is set initially and doesn't change */ + if (pcrd->rate != 0) + { + pcrd->value_ref = pcrd->init + pcrd->rate*t; + } +} + +/* Calculates pull->coord[coord_ind].value. + * This function also updates pull->coord[coord_ind].dr. + */ +static void get_pull_coord_distance(t_pull *pull, + int coord_ind, + const t_pbc *pbc) +{ + t_pull_coord *pcrd; + int m; pcrd = &pull->coord[coord_ind]; - get_pull_coord_dr(pull, coord_ind, pbc, t); + get_pull_coord_dr(pull, coord_ind, pbc); - ref = pcrd->init + pcrd->rate*t; + switch (pcrd->eGeom) + { + case epullgDIST: + /* Pull along the vector between the com's */ + pcrd->value = dnorm(pcrd->dr); + break; + case epullgDIR: + case epullgDIRPBC: + case epullgDIRRELATIVE: + case epullgCYL: + /* Pull along vec */ + pcrd->value = 0; + for (m = 0; m < DIM; m++) + { + pcrd->value += pcrd->vec[m]*pcrd->dr[m]; + } + break; + } +} + +/* Returns the deviation from the reference value. + * Updates pull->coord[coord_ind].dr, .value and .value_ref. + */ +static double get_pull_coord_deviation(t_pull *pull, + int coord_ind, + const t_pbc *pbc, + double t) +{ + static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety, + but is fairly benign */ + t_pull_coord *pcrd; + double dev = 0; + + pcrd = &pull->coord[coord_ind]; + + get_pull_coord_distance(pull, coord_ind, pbc); + + update_pull_coord_reference_value(pcrd, t); switch (pcrd->eGeom) { case epullgDIST: /* Pull along the vector between the com's */ - if (ref < 0 && !bWarned) + if (pcrd->value_ref < 0 && !bWarned) { - fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref); + fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref); bWarned = TRUE; } - drs = dnorm(pcrd->dr); - if (drs == 0) + if (pcrd->value == 0) { /* With no vector we can not determine the direction for the force, * so we set the force to zero. */ - *dev = 0; + dev = 0; } else { /* Determine the deviation */ - *dev = drs - ref; + dev = pcrd->value - pcrd->value_ref; } break; case epullgDIR: @@ -625,14 +643,21 @@ void get_pull_coord_distance(const t_pull *pull, case epullgDIRRELATIVE: case epullgCYL: /* Pull along vec */ - inpr = 0; - for (m = 0; m < DIM; m++) - { - inpr += pcrd->vec[m]*pcrd->dr[m]; - } - *dev = inpr - ref; + dev = pcrd->value - pcrd->value_ref; break; } + + return dev; +} + +void get_pull_coord_value(t_pull *pull, + int coord_ind, + const t_pbc *pbc, + double *value) +{ + get_pull_coord_distance(pull, coord_ind, pbc); + + *value = pull->coord[coord_ind].value; } void clear_pull_forces(t_pull *pull) @@ -661,7 +686,6 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */ dvec *rnew; /* current 'new' positions of the groups */ double *dr_tot; /* the total update of the coords */ - double ref; dvec vec; double inpr; double lambda, rm, invdt = 0; @@ -696,8 +720,12 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, continue; } - get_pull_coord_dr(pull, c, pbc, t); - /* Store the difference vector at time t for printing */ + /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value. + * We don't modify dr and value anymore, so these values are also used + * for printing. + */ + get_pull_coord_distance(pull, c, pbc); + if (debug) { fprintf(debug, "Pull coord %d dr %f %f %f\n", @@ -745,17 +773,17 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, continue; } + update_pull_coord_reference_value(pcrd, t); + pgrp0 = &pull->group[pcrd->group[0]]; pgrp1 = &pull->group[pcrd->group[1]]; /* Get the current difference vector */ - low_get_pull_coord_dr(pull, pcrd, pbc, t, + low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->group[1]], rnew[pcrd->group[0]], -1, unc_ij); - ref = pcrd->init + pcrd->rate*t; - if (debug) { fprintf(debug, "Pull coord %d, iteration %d\n", c, niter); @@ -766,9 +794,9 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, switch (pcrd->eGeom) { case epullgDIST: - if (ref <= 0) + if (pcrd->value_ref <= 0) { - gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref); + gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref); } { @@ -776,7 +804,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, c_a = diprod(r_ij[c], r_ij[c]); c_b = diprod(unc_ij, r_ij[c])*2; - c_c = diprod(unc_ij, unc_ij) - dsqr(ref); + c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref); if (c_b < 0) { @@ -814,7 +842,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, } /* Select only the component along the vector */ dsvmul(a, vec, unc_ij); - lambda = a - ref; + lambda = a - pcrd->value_ref; if (debug) { fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda); @@ -834,15 +862,15 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, g0 = pcrd->group[0]; g1 = pcrd->group[1]; - low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp); - low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3); + low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp); + low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3); fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0], rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp)); fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", - "", "", "", "", "", "", ref); + "", "", "", "", "", "", pcrd->value_ref); fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", dr0[0], dr0[1], dr0[2], @@ -865,9 +893,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, continue; } - ref = pcrd->init + pcrd->rate*t; - - low_get_pull_coord_dr(pull, pcrd, pbc, t, + low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->group[1]], rnew[pcrd->group[0]], -1, unc_ij); @@ -875,7 +901,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, switch (pcrd->eGeom) { case epullgDIST: - bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol; + bConverged = + fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol; break; case epullgDIR: case epullgDIRPBC: @@ -887,7 +914,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, inpr = diprod(unc_ij, vec); dsvmul(inpr, vec, unc_ij); bConverged = - fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol; + fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol; break; } @@ -897,7 +924,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, { fprintf(debug, "NOT CONVERGED YET: Group %d:" "d_ref = %f, current d = %f\n", - g, ref, dnorm(unc_ij)); + g, pcrd->value_ref, dnorm(unc_ij)); } bConverged_all = FALSE; @@ -1017,7 +1044,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, continue; } - get_pull_coord_distance(pull, c, pbc, t, &dev); + dev = get_pull_coord_deviation(pull, c, pbc, t); k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB; dkdl = pcrd->kB - pcrd->k; @@ -1424,6 +1451,13 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], { pull->group[pcrd->group[g]].bCalcCOM = TRUE; } + + /* With non-zero rate the reference value is set at every step */ + if (pcrd->rate == 0) + { + /* Initialize the constant reference value */ + pcrd->value_ref = pcrd->init; + } } pull->ePBC = ir->ePBC; diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index 72a0f8bfe6..f1f49bd7ac 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -59,20 +59,18 @@ extern "C" { struct t_pbc; -/*! \brief Get the deviation for pull coord coord_ind. - * - * This call updates some data in the pull coordinates in \p pull. + +/*! \brief Get the value for pull coord coord_ind. * * \param[in,out] pull The pull struct. * \param[in] coord_ind Number of the pull coordinate. * \param[in] pbc Information structure about periodicity. - * \param[in] t Time. - * \param[out] dev The deviation from the reference distance. + * \param[out] value The value of the pull coordinate. */ -void get_pull_coord_distance(const t_pull *pull, - int coord_ind, - const struct t_pbc *pbc, double t, - double *dev); +void get_pull_coord_value(t_pull *pull, + int coord_ind, + const struct t_pbc *pbc, + double *value); /*! \brief Set the all the pull forces to zero. diff --git a/src/gromacs/pulling/pullutil.c b/src/gromacs/pulling/pullutil.c index 135edf649c..abfcda8ae6 100644 --- a/src/gromacs/pulling/pullutil.c +++ b/src/gromacs/pulling/pullutil.c @@ -158,9 +158,14 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, * we reduced the other group COM over the ranks. This resolves * any PBC issues and we don't need to use a PBC-atom here. */ + if (pcrd->rate != 0) + { + /* With rate=0, value_ref is set initially */ + pcrd->value_ref = pcrd->init + pcrd->rate*t; + } for (m = 0; m < DIM; m++) { - g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t); + g_x[m] = pgrp->x[m] - pcrd->vec[m]*pcrd->value_ref; } /* loop over all atoms in the main ref group */ @@ -280,7 +285,7 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, pcrd->cyl_dev = 0; for (m = 0; m < DIM; m++) { - g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t); + g_x[m] = pgrp->x[m] - pcrd->vec[m]*pcrd->value_ref; dist = -pcrd->vec[m]*pull->dbuf_cyl[c*stride+2]*pdyna->mwscale; pdyna->x[m] = g_x[m] - dist; pcrd->cyl_dev += dist;