src/gromacs/pulling/output.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "output.h"
  40
  41 #include <cstdio>
  42
  43 #include "gromacs/commandline/filenm.h"
  44 #include "gromacs/fileio/gmxfio.h"
  45 #include "gromacs/fileio/xvgr.h"
  46 #include "gromacs/pulling/pull.h"
  47 #include "gromacs/utility/cstringutil.h"
  48 #include "gromacs/utility/exceptions.h"
  49 #include "gromacs/utility/fatalerror.h"
  50
  51 #include "pull_internal.h"
  52
  53 static std::string append_before_extension(const std::string &pathname,
  54                                            const std::string &to_append)
  55 {
  56     /* Appends to_append before last '.' in pathname */
  57     size_t extPos = pathname.find_last_of('.');
  58     if (extPos == std::string::npos)
  59     {
  60         return pathname + to_append;
  61     }
  62     else
  63     {
  64         return pathname.substr(0, extPos) + to_append +
  65                pathname.substr(extPos, std::string::npos);
  66     }
  67 }
  68
  69 static void pull_print_group_x(FILE *out, const ivec dim,
  70                                const pull_group_work_t *pgrp)
  71 {
  72     int m;
  73
  74     for (m = 0; m < DIM; m++)
  75     {
  76         if (dim[m])
  77         {
  78             fprintf(out, "\t%g", pgrp->x[m]);
  79         }
  80     }
  81 }
  82
  83 static void pull_print_coord_dr_components(FILE *out, const ivec dim, const dvec dr)
  84 {
  85     for (int m = 0; m < DIM; m++)
  86     {
  87         if (dim[m])
  88         {
  89             fprintf(out, "\t%g", dr[m]);
  90         }
  91     }
  92 }
  93
  94 static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd,
  95                                 gmx_bool bPrintRefValue,
  96                                 gmx_bool bPrintComponents)
  97 {
  98     double unit_factor = pull_conversion_factor_internal2userinput(&pcrd->params);
  99
 100     fprintf(out, "\t%g", pcrd->spatialData.value*unit_factor);
 101
 102     if (bPrintRefValue)
 103     {
 104         fprintf(out, "\t%g", pcrd->value_ref*unit_factor);
 105     }
 106
 107     if (bPrintComponents)
 108     {
 109         pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr01);
 110         if (pcrd->params.ngroup >= 4)
 111         {
 112             pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr23);
 113         }
 114         if (pcrd->params.ngroup >= 6)
 115         {
 116             pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr45);
 117         }
 118     }
 119 }
 120
 121 static void pull_print_x(FILE *out, struct pull_t *pull, double t)
 122 {
 123     fprintf(out, "%.4f", t);
 124
 125     for (size_t c = 0; c < pull->coord.size(); c++)
 126     {
 127         pull_coord_work_t *pcrd;
 128
 129         pcrd = &pull->coord[c];
 130
 131         pull_print_coord_dr(out, pcrd,
 132                             pull->params.bPrintRefValue && pcrd->params.eType != epullEXTERNAL,
 133                             pull->params.bPrintComp);
 134
 135         if (pull->params.bPrintCOM)
 136         {
 137             if (pcrd->params.eGeom == epullgCYL)
 138             {
 139                 pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]);
 140             }
 141             else
 142             {
 143                 pull_print_group_x(out, pcrd->params.dim,
 144                                    &pull->group[pcrd->params.group[0]]);
 145             }
 146             for (int g = 1; g < pcrd->params.ngroup; g++)
 147             {
 148                 pull_print_group_x(out, pcrd->params.dim, &pull->group[pcrd->params.group[g]]);
 149             }
 150         }
 151     }
 152     fprintf(out, "\n");
 153 }
 154
 155 static void pull_print_f(FILE *out, const pull_t *pull, double t)
 156 {
 157     fprintf(out, "%.4f", t);
 158
 159     for (const pull_coord_work_t &coord : pull->coord)
 160     {
 161         fprintf(out, "\t%g", coord.scalarForce);
 162     }
 163     fprintf(out, "\n");
 164 }
 165
 166 void pull_print_output(struct pull_t *pull, int64_t step, double time)
 167 {
 168     GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "pull_print_output called before all external pull potentials have been applied");
 169
 170     if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0))
 171     {
 172         pull_print_x(pull->out_x, pull, time);
 173     }
 174
 175     if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0))
 176     {
 177         pull_print_f(pull->out_f, pull, time);
 178     }
 179 }
 180
 181 static void set_legend_for_coord_components(const pull_coord_work_t *pcrd, int coord_index, char **setname, int *nsets_ptr)
 182 {
 183     /*  Loop over the distance vectors and print their components. Each vector is made up of two consecutive groups. */
 184     for (int g = 0; g < pcrd->params.ngroup; g += 2)
 185     {
 186         /* Loop over the components */
 187         for (int m  = 0; m < DIM; m++)
 188         {
 189             if (pcrd->params.dim[m])
 190             {
 191                 char legend[STRLEN];
 192
 193                 if (g == 0 && pcrd->params.ngroup <= 2)
 194                 {
 195                     /*  For the simplest case we print a simplified legend without group indices, just the cooordinate index
 196                         and which dimensional component it is. */
 197                     sprintf(legend, "%d d%c", coord_index + 1, 'X' + m);
 198                 }
 199                 else
 200                 {
 201                     /* Otherwise, print also the group indices (relative to the coordinate, not the global ones). */
 202                     sprintf(legend, "%d g %d-%d d%c", coord_index + 1, g + 1, g + 2, 'X' + m);
 203                 }
 204
 205                 setname[*nsets_ptr] = gmx_strdup(legend);
 206                 (*nsets_ptr)++;
 207             }
 208         }
 209     }
 210 }
 211
 212 static FILE *open_pull_out(const char *fn, struct pull_t *pull,
 213                            const gmx_output_env_t *oenv,
 214                            gmx_bool bCoord,
 215                            const ContinuationOptions &continuationOptions)
 216 {
 217     FILE  *fp;
 218     int    nsets, m;
 219     char **setname, buf[50];
 220
 221     if (continuationOptions.appendFiles)
 222     {
 223         fp = gmx_fio_fopen(fn, "a+");
 224     }
 225     else
 226     {
 227         fp = gmx_fio_fopen(fn, "w+");
 228         if (bCoord)
 229         {
 230             sprintf(buf, "Position (nm%s)", pull->bAngle ? ", deg" : "");
 231             xvgr_header(fp, "Pull COM",  "Time (ps)", buf,
 232                         exvggtXNY, oenv);
 233         }
 234         else
 235         {
 236             sprintf(buf, "Force (kJ/mol/nm%s)", pull->bAngle ? ", kJ/mol/rad" : "");
 237             xvgr_header(fp, "Pull force", "Time (ps)", buf,
 238                         exvggtXNY, oenv);
 239         }
 240
 241         /* With default mdp options only the actual coordinate value is printed (1),
 242          * but optionally the reference value (+ 1),
 243          * the group COMs for all the groups (+ ngroups_max*DIM)
 244          * and the components of the distance vectors can be printed (+ (ngroups_max/2)*DIM).
 245          */
 246         snew(setname, pull->coord.size()*(1 + 1 + c_pullCoordNgroupMax*DIM + c_pullCoordNgroupMax/2*DIM));
 247
 248         nsets = 0;
 249         for (size_t c = 0; c < pull->coord.size(); c++)
 250         {
 251             if (bCoord)
 252             {
 253                 /* The order of this legend should match the order of printing
 254                  * the data in print_pull_x above.
 255                  */
 256
 257                 /* The pull coord distance */
 258                 sprintf(buf, "%zu", c+1);
 259                 setname[nsets] = gmx_strdup(buf);
 260                 nsets++;
 261                 if (pull->params.bPrintRefValue &&
 262                     pull->coord[c].params.eType != epullEXTERNAL)
 263                 {
 264                     sprintf(buf, "%zu ref", c+1);
 265                     setname[nsets] = gmx_strdup(buf);
 266                     nsets++;
 267                 }
 268                 if (pull->params.bPrintComp)
 269                 {
 270                     set_legend_for_coord_components(&pull->coord[c], c, setname, &nsets);
 271                 }
 272
 273                 if (pull->params.bPrintCOM)
 274                 {
 275                     for (int g = 0; g < pull->coord[c].params.ngroup; g++)
 276                     {
 277                         /* Legend for reference group position */
 278                         for (m = 0; m < DIM; m++)
 279                         {
 280                             if (pull->coord[c].params.dim[m])
 281                             {
 282                                 sprintf(buf, "%zu g %d %c", c+1, g + 1, 'X'+m);
 283                                 setname[nsets] = gmx_strdup(buf);
 284                                 nsets++;
 285                             }
 286                         }
 287                     }
 288                 }
 289             }
 290             else
 291             {
 292                 /* For the pull force we always only use one scalar */
 293                 sprintf(buf, "%zu", c+1);
 294                 setname[nsets] = gmx_strdup(buf);
 295                 nsets++;
 296             }
 297         }
 298         if (nsets > 1)
 299         {
 300             xvgr_legend(fp, nsets, setname, oenv);
 301         }
 302         for (int c = 0; c < nsets; c++)
 303         {
 304             sfree(setname[c]);
 305         }
 306         sfree(setname);
 307     }
 308
 309     return fp;
 310 }
 311
 312 void init_pull_output_files(pull_t                    *pull,
 313                             int                        nfile,
 314                             const t_filenm             fnm[],
 315                             const gmx_output_env_t    *oenv,
 316                             const ContinuationOptions &continuationOptions)
 317 {
 318     /* Check for px and pf filename collision, if we are writing
 319        both files */
 320     std::string px_filename, pf_filename;
 321     std::string px_appended, pf_appended;
 322     try
 323     {
 324         px_filename  = std::string(opt2fn("-px", nfile, fnm));
 325         pf_filename  = std::string(opt2fn("-pf", nfile, fnm));
 326     }
 327     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 328
 329     if ((pull->params.nstxout != 0) &&
 330         (pull->params.nstfout != 0) &&
 331         (px_filename == pf_filename))
 332     {
 333         if (!opt2bSet("-px", nfile, fnm) && !opt2bSet("-pf", nfile, fnm))
 334         {
 335             /* We are writing both pull files but neither set directly. */
 336             try
 337             {
 338                 px_appended   = append_before_extension(px_filename, "_pullx");
 339                 pf_appended   = append_before_extension(pf_filename, "_pullf");
 340             }
 341             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 342             pull->out_x = open_pull_out(px_appended.c_str(), pull, oenv,
 343                                         TRUE, continuationOptions);
 344             pull->out_f = open_pull_out(pf_appended.c_str(), pull, oenv,
 345                                         FALSE, continuationOptions);
 346             return;
 347         }
 348         else
 349         {
 350             /* If at least one of -px and -pf is set but the filenames are identical: */
 351             gmx_fatal(FARGS, "Identical pull_x and pull_f output filenames %s",
 352                       px_filename.c_str());
 353         }
 354     }
 355     if (pull->params.nstxout != 0)
 356     {
 357         pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
 358                                     TRUE, continuationOptions);
 359     }
 360     if (pull->params.nstfout != 0)
 361     {
 362         pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
 363                                     FALSE, continuationOptions);
 364     }
 365 }