From bec9c8757e59cae58fc61ed841c0bb73c84079db Mon Sep 17 00:00:00 2001 From: Viveca Lindahl Date: Mon, 21 Sep 2015 18:18:38 +0200 Subject: [PATCH] Added pull coordinate geometry angle-axis The new geometry is described in the docs. Some checks in readpull.cpp where reorganized since adding new geometries made some old logic a bit convoluted. Change-Id: I310891740d65e3c47948a06726e7a9e7c974bc44 --- docs/manual/special.tex | 1 + docs/user-guide/mdp-options.rst | 5 +++ src/gromacs/gmxana/gmx_wham.cpp | 1 + src/gromacs/gmxpreprocess/readpull.cpp | 63 ++++++++++++++++++---------------- src/gromacs/math/vec.h | 16 +++++++++ src/gromacs/mdtypes/md_enums.cpp | 2 +- src/gromacs/mdtypes/md_enums.h | 2 +- src/gromacs/pulling/pull.cpp | 46 ++++++++++++++++++++++--- src/gromacs/pulling/pull_internal.h | 2 +- src/gromacs/pulling/pullutil.cpp | 2 +- 10 files changed, 101 insertions(+), 39 deletions(-) diff --git a/docs/manual/special.tex b/docs/manual/special.tex index 14ce7da762..aa6b959f16 100644 --- a/docs/manual/special.tex +++ b/docs/manual/special.tex @@ -332,6 +332,7 @@ $i$ and $i+1$ define a vector connecting the COMs of groups $i$ and $i+1$. The a E.g., the {\tt .mdp} option {\tt pull-coord?-groups = 1 2 2 4} defines the angle between the vector from the COM of group 1 to the COM of group 2 and the vector from the COM of group 2 to the COM of group 4. The angle takes values in the closed interval [0, 180] deg. +For {\tt pull-coord?-geometry = angle-axis} the angle is defined with respect to a reference axis given by {\tt pull-coord?-vec} and only two groups need to be given. The dihedral geometry requires six pull groups. These pair up in the same way as described above and so define three vectors. The dihedral angle is defined as the angle between the two planes spanned by the two first and the two last vectors. Equivalently, the dihedral angle can be seen as the angle between the first and the third vector when these vectors are projected onto a plane normal to the second vector (the axis vector). diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index a4cbd0aa35..ddd63f95e1 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -1700,6 +1700,11 @@ applicable pulling coordinate. the vector connecting the COM of the third group to the COM of the fourth group. + .. mdp-value:: angle-axis + + As :mdp-value:`angle` but the second vector is given by :mdp:`pull-coord1-vec`. + Thus, only the two groups that define the first vector need to be given. + .. mdp-value:: dihedral Pull along a dihedral angle defined by six groups. These pairwise diff --git a/src/gromacs/gmxana/gmx_wham.cpp b/src/gromacs/gmxana/gmx_wham.cpp index ef82cdf55a..15a8f4c034 100644 --- a/src/gromacs/gmxana/gmx_wham.cpp +++ b/src/gromacs/gmxana/gmx_wham.cpp @@ -2157,6 +2157,7 @@ void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions case epullgDIRPBC: case epullgANGLE: case epullgDIHEDRAL: + case epullgANGLEAXIS: header->umbInitDist[i] = header->init_dist[i]; break; default: diff --git a/src/gromacs/gmxpreprocess/readpull.cpp b/src/gromacs/gmxpreprocess/readpull.cpp index 5bc6ce7c65..15da010e02 100644 --- a/src/gromacs/gmxpreprocess/readpull.cpp +++ b/src/gromacs/gmxpreprocess/readpull.cpp @@ -123,9 +123,10 @@ static void process_pull_dim(char *dim_buf, ivec dim, const t_pull_coord *pcrd) { gmx_fatal(FARGS, "Pull geometry dihedral is only useful with pull-dim = Y Y Y"); } - if ((pcrd->eGeom == epullgANGLE) && (ndim < 2)) + if ((pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgANGLEAXIS ) && (ndim < 2)) { - gmx_fatal(FARGS, "Pull geometry angle is only useful with pull-dim = Y for at least 2 dimensions"); + gmx_fatal(FARGS, "Pull geometry %s is only useful with pull-dim = Y for at least 2 dimensions", + EPULLGEOM(pcrd->eGeom)); } } @@ -141,6 +142,7 @@ static void init_pull_coord(t_pull_coord *pcrd, if (pcrd->eType == epullCONSTRAINT && (pcrd->eGeom == epullgCYL || pcrd->eGeom == epullgDIRRELATIVE || pcrd->eGeom == epullgANGLE || + pcrd->eGeom == epullgANGLEAXIS || pcrd->eGeom == epullgDIHEDRAL)) { gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.", @@ -157,7 +159,7 @@ static void init_pull_coord(t_pull_coord *pcrd, gmx_fatal(FARGS, "The pull origin can only be set with an absolute reference"); } - clear_dvec(vec); + /* Check the given initial reference value and warn for dangerous values */ if (pcrd->eGeom == epullgDIST) { if (pcrd->bStart && pcrd->init < 0) @@ -169,7 +171,7 @@ static void init_pull_coord(t_pull_coord *pcrd, warning(wi, buf); } } - else if (pcrd->eGeom == epullgANGLE) + else if (pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgANGLEAXIS) { if (pcrd->bStart && (pcrd->init < 0 || pcrd->init > 180)) { @@ -190,19 +192,38 @@ static void init_pull_coord(t_pull_coord *pcrd, warning(wi, buf); } } - else if (pcrd->eGeom != epullgDIRRELATIVE) + + /* Check and set the pull vector */ + clear_dvec(vec); + string2dvec(vec_buf, vec); + + if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgCYL || pcrd->eGeom == epullgDIRPBC || pcrd->eGeom == epullgANGLEAXIS) { - /* Check and set the pull vector */ - string2dvec(vec_buf, vec); if (dnorm2(vec) == 0) { gmx_fatal(FARGS, "With pull geometry %s the pull vector can not be 0,0,0", epullg_names[pcrd->eGeom]); } - if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgCYL) + for (int d = 0; d < DIM; d++) + { + if (vec[d] != 0 && pcrd->dim[d] == 0) + { + gmx_fatal(FARGS, "pull-coord-vec has non-zero %c-component while pull_dim for the %c-dimension is set to N", 'x'+d, 'x'+d); + } + } + + /* Normalize the direction vector */ + dsvmul(1/dnorm(vec), vec, vec); + } + else /* This case is for are all the geometries where the pull vector is not used */ + { + if (dnorm2(vec) > 0) { - /* Normalize the direction vector */ - dsvmul(1/dnorm(vec), vec, vec); + sprintf(buf, "A pull vector is given (%g %g %g) but will not be used with geometry %s. If you really want to use this " + "vector, consider using geometry %s instead.", + vec[0], vec[1], vec[2], EPULLGEOM(pcrd->eGeom), + pcrd->eGeom == epullgANGLE ? EPULLGEOM(epullgANGLEAXIS) : EPULLGEOM(epullgDIR)); + warning(wi, buf); } } for (m = 0; m < DIM; m++) @@ -410,7 +431,7 @@ void make_pull_groups(pull_params_t *pull, void make_pull_coords(pull_params_t *pull) { - int c, d; + int c; t_pull_coord *pcrd; for (c = 0; c < pull->ncoord; c++) @@ -435,24 +456,6 @@ void make_pull_coords(pull_params_t *pull) gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling"); } } - - if (pcrd->eGeom != epullgDIST) - { - for (d = 0; d < DIM; d++) - { - if (pcrd->vec[d] != 0 && pcrd->dim[d] == 0) - { - gmx_fatal(FARGS, "ERROR: pull-group%d-vec has non-zero %c-component while pull_dim for the %c-dimension is N\n", c+1, 'x'+d, 'x'+d); - } - } - } - - if ((pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgCYL) && - norm2(pcrd->vec) == 0) - { - gmx_fatal(FARGS, "pull-group%d-vec can not be zero with geometry %s", - c+1, EPULLGEOM(pcrd->eGeom)); - } } } @@ -530,7 +533,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l * generalization of the pull code makes pull dim available here. */ } - else if (pcrd->eGeom == epullgANGLE) + else if (pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgANGLEAXIS) { if (pcrd->init < 0 || pcrd->init > 180) { diff --git a/src/gromacs/math/vec.h b/src/gromacs/math/vec.h index a67dc5df61..c842597338 100644 --- a/src/gromacs/math/vec.h +++ b/src/gromacs/math/vec.h @@ -52,6 +52,8 @@ void rvec_dec(rvec a,rvec b) a -= b void copy_rvec(const rvec a,rvec b) b = a (reals) void copy_dvec(const dvec a,dvec b) b = a (reals) + void copy_rvec_to_dvec(const rvec a,dvec b) b = a (reals) + void copy_dvec_to_rvec(const dvec a,rvec b) b = a (reals) void copy_ivec(const ivec a,ivec b) b = a (integers) void ivec_sub(const ivec a,const ivec b,ivec c) c = a - b void svmul(real a,rvec v1,rvec v2) v2 = a * v1 @@ -214,6 +216,20 @@ static inline void copy_rvec(const rvec a, rvec b) b[ZZ] = a[ZZ]; } +static inline void copy_rvec_to_dvec(const rvec a, dvec b) +{ + b[XX] = a[XX]; + b[YY] = a[YY]; + b[ZZ] = a[ZZ]; +} + +static inline void copy_dvec_to_rvec(const dvec a, rvec b) +{ + b[XX] = a[XX]; + b[YY] = a[YY]; + b[ZZ] = a[ZZ]; +} + static inline void copy_rvecn(const rvec *a, rvec *b, int startn, int endn) { int i; diff --git a/src/gromacs/mdtypes/md_enums.cpp b/src/gromacs/mdtypes/md_enums.cpp index 9b3e0accf4..6c0be143b8 100644 --- a/src/gromacs/mdtypes/md_enums.cpp +++ b/src/gromacs/mdtypes/md_enums.cpp @@ -210,7 +210,7 @@ const char *epull_names[epullNR+1] = { }; const char *epullg_names[epullgNR+1] = { - "distance", "direction", "cylinder", "direction-periodic", "direction-relative", "angle", "dihedral", nullptr + "distance", "direction", "cylinder", "direction-periodic", "direction-relative", "angle", "dihedral", "angle-axis", nullptr }; const char *erotg_names[erotgNR+1] = { diff --git a/src/gromacs/mdtypes/md_enums.h b/src/gromacs/mdtypes/md_enums.h index fafbf16b38..462910511c 100644 --- a/src/gromacs/mdtypes/md_enums.h +++ b/src/gromacs/mdtypes/md_enums.h @@ -516,7 +516,7 @@ extern const char *epull_names[epullNR+1]; //! Control of pull groups enum { - epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgDIRRELATIVE, epullgANGLE, epullgDIHEDRAL, epullgNR + epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgDIRRELATIVE, epullgANGLE, epullgDIHEDRAL, epullgANGLEAXIS, epullgNR }; //! String for pull groups extern const char *epullg_names[epullgNR+1]; diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index fccd73db12..e6f9d8a5c1 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -78,7 +78,8 @@ static bool pull_coordinate_is_angletype(const t_pull_coord *pcrd) { return (pcrd->eGeom == epullgANGLE || - pcrd->eGeom == epullgDIHEDRAL); + pcrd->eGeom == epullgDIHEDRAL || + pcrd->eGeom == epullgANGLEAXIS); } const char *pull_coordinate_units(const t_pull_coord *pcrd) @@ -724,7 +725,7 @@ static void low_set_pull_coord_reference_value(pull_coord_work_t *pcrd, gmx_fatal(FARGS, "Pull reference distance for coordinate %d (%f) needs to be non-negative", coord_ind + 1, value_ref); } } - else if (pcrd->params.eGeom == epullgANGLE) + else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS) { if (value_ref < 0 || value_ref > M_PI) { @@ -808,6 +809,9 @@ static void get_pull_coord_distance(struct pull_t *pull, case epullgDIHEDRAL: pcrd->value = get_dihedral_angle_coord(pcrd); break; + case epullgANGLEAXIS: + pcrd->value = gmx_angle_between_dvecs(pcrd->dr01, pcrd->vec); + break; default: gmx_incons("Unsupported pull type in get_pull_coord_distance"); } @@ -1359,6 +1363,36 @@ static void calc_pull_coord_force(pull_coord_work_t *pcrd, clear_pull_forces_coord(pcrd); } } + else if (pcrd->params.eGeom == epullgANGLEAXIS) + { + double cos_theta, cos_theta2; + + /* The angle-axis force is exactly the same as the angle force (above) except that in + this case the second vector (dr23) is replaced by the pull vector. */ + cos_theta = cos(pcrd->value); + cos_theta2 = gmx::square(cos_theta); + + if (cos_theta2 < 1) + { + double a, b; + double invdr01; + dvec normalized_dr01; + + invdr01 = 1./dnorm(pcrd->dr01); + dsvmul(invdr01, pcrd->dr01, normalized_dr01); + a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */ + b = a*cos_theta; + + for (m = 0; m < DIM; m++) + { + pcrd->f01[m] = pcrd->f_scal*invdr01*(a*pcrd->vec[m] - b*normalized_dr01[m]); + } + } + else + { + clear_pull_forces_coord(pcrd); + } + } else if (pcrd->params.eGeom == epullgDIHEDRAL) { double m2, n2, tol, sqrdist_32; @@ -1966,7 +2000,8 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, case epullgDIR: case epullgDIRPBC: case epullgCYL: - copy_rvec(pull_params->coord[c].vec, pcrd->vec); + case epullgANGLEAXIS: + copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->vec); break; default: /* We allow reading of newer tpx files with new pull geometries, @@ -1988,7 +2023,8 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE || pcrd->params.eGeom == epullgANGLE || - pcrd->params.eGeom == epullgDIHEDRAL) + pcrd->params.eGeom == epullgDIHEDRAL || + pcrd->params.eGeom == epullgANGLEAXIS) { gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.", epull_names[pcrd->params.eType], @@ -2007,7 +2043,7 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, { pull->bCylinder = TRUE; } - else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL) + else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL || pcrd->params.eGeom == epullgANGLEAXIS) { pull->bAngle = TRUE; } diff --git a/src/gromacs/pulling/pull_internal.h b/src/gromacs/pulling/pull_internal.h index 7b59301151..65c2707969 100644 --- a/src/gromacs/pulling/pull_internal.h +++ b/src/gromacs/pulling/pull_internal.h @@ -90,7 +90,7 @@ typedef struct dvec dr01; /* The direction vector of group 1 relative to group 0 */ dvec dr23; /* The direction vector of group 3 relative to group 2 */ dvec dr45; /* The direction vector of group 5 relative to group 4 */ - rvec vec; /* The pull direction */ + dvec vec; /* The pull direction */ double vec_len; /* Length of vec for direction-relative */ dvec ffrad; /* conversion factor from vec to radial force */ double cyl_dev; /* The deviation from the reference position */ diff --git a/src/gromacs/pulling/pullutil.cpp b/src/gromacs/pulling/pullutil.cpp index ba3434cdf7..503b40f55a 100644 --- a/src/gromacs/pulling/pullutil.cpp +++ b/src/gromacs/pulling/pullutil.cpp @@ -243,7 +243,7 @@ static void make_cyl_refgrps(t_commrec *cr, struct pull_t *pull, t_mdatoms *md, pref = &pull->group[pcrd->params.group[0]]; pgrp = &pull->group[pcrd->params.group[1]]; pdyna = &pull->dyna[c]; - copy_rvec(pcrd->vec, dir); + copy_dvec_to_rvec(pcrd->vec, dir); pdyna->nat_loc = 0; /* We calculate distances with respect to the reference location -- 2.11.4.GIT