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,2019, 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 <memory>
  44
  45 #include "gromacs/commandline/filenm.h"
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/mdtypes/mdrunoptions.h"
  50 #include "gromacs/mdtypes/observableshistory.h"
  51 #include "gromacs/mdtypes/pullhistory.h"
  52 #include "gromacs/pulling/pull.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/exceptions.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/gmxassert.h"
  57
  58 #include "pull_internal.h"
  59
  60 static std::string append_before_extension(const std::string &pathname,
  61                                            const std::string &to_append)
  62 {
  63     /* Appends to_append before last '.' in pathname */
  64     size_t extPos = pathname.find_last_of('.');
  65     if (extPos == std::string::npos)
  66     {
  67         return pathname + to_append;
  68     }
  69     else
  70     {
  71         return pathname.substr(0, extPos) + to_append +
  72                pathname.substr(extPos, std::string::npos);
  73     }
  74 }
  75
  76 static void addToPullxHistory(pull_t *pull)
  77 {
  78     GMX_RELEASE_ASSERT(pull->coordForceHistory, "Pull history does not exist.");
  79     pull->coordForceHistory->numValuesInXSum++;
  80     for (size_t c = 0; c < pull->coord.size(); c++)
  81     {
  82         const pull_coord_work_t &pcrd        = pull->coord[c];
  83         PullCoordinateHistory   &pcrdHistory = pull->coordForceHistory->pullCoordinateSums[c];
  84
  85         pcrdHistory.value    += pcrd.spatialData.value;
  86         pcrdHistory.valueRef += pcrd.value_ref;
  87
  88         for (int m = 0; m < DIM; m++)
  89         {
  90             pcrdHistory.dr01[m] += pcrd.spatialData.dr01[m];
  91             pcrdHistory.dr23[m] += pcrd.spatialData.dr23[m];
  92             pcrdHistory.dr45[m] += pcrd.spatialData.dr45[m];
  93         }
  94         if (pcrd.params.eGeom == epullgCYL)
  95         {
  96             for (int m = 0; m < DIM; m++)
  97             {
  98                 pcrdHistory.dynaX[m] += pull->dyna[c].x[m];
  99             }
 100         }
 101     }
 102     for (size_t g = 0; g < pull->group.size(); g++)
 103     {
 104         PullGroupHistory &pgrpHistory = pull->coordForceHistory->pullGroupSums[g];
 105         for (int m = 0; m < DIM; m++)
 106         {
 107             pgrpHistory.x[m] += pull->group[g].x[m];
 108         }
 109     }
 110 }
 111
 112 static void addToPullfHistory(pull_t *pull)
 113 {
 114     GMX_RELEASE_ASSERT(pull->coordForceHistory, "Pull history does not exist.");
 115     pull->coordForceHistory->numValuesInFSum++;
 116     for (size_t c = 0; c < pull->coord.size(); c++)
 117     {
 118         const pull_coord_work_t &pcrd        = pull->coord[c];;
 119         PullCoordinateHistory   &pcrdHistory = pull->coordForceHistory->pullCoordinateSums[c];
 120
 121         pcrdHistory.scalarForce += pcrd.scalarForce;
 122     }
 123 }
 124
 125 static void pullResetHistory(PullHistory *history,
 126                              bool         resetXHistory,
 127                              bool         resetFHistory)
 128 {
 129     if (resetXHistory)
 130     {
 131         history->numValuesInXSum = 0;
 132
 133         for (PullCoordinateHistory &pcrdHistory : history->pullCoordinateSums)
 134         {
 135             pcrdHistory.value        = 0;
 136             pcrdHistory.valueRef     = 0;
 137
 138             clear_dvec(pcrdHistory.dr01);
 139             clear_dvec(pcrdHistory.dr23);
 140             clear_dvec(pcrdHistory.dr45);
 141             clear_dvec(pcrdHistory.dynaX);
 142         }
 143
 144         for (PullGroupHistory &pgrpHistory : history->pullGroupSums)
 145         {
 146             clear_dvec(pgrpHistory.x);
 147         }
 148     }
 149     if (resetFHistory)
 150     {
 151         history->numValuesInFSum = 0;
 152         for (PullCoordinateHistory &pcrdHistory : history->pullCoordinateSums)
 153         {
 154             pcrdHistory.scalarForce = 0;
 155         }
 156     }
 157 }
 158
 159 static void pull_print_coord_dr_components(FILE *out, const ivec dim, const dvec dr,
 160                                            const int numValuesInSum)
 161 {
 162     for (int m = 0; m < DIM; m++)
 163     {
 164         if (dim[m])
 165         {
 166             fprintf(out, "\t%g", dr[m] / numValuesInSum);
 167         }
 168     }
 169 }
 170
 171 template <typename T>
 172 static void pull_print_coord_dr(FILE                *out,
 173                                 const pull_params_t &pullParams,
 174                                 const t_pull_coord  &coordParams,
 175                                 const T             &pcrdData,
 176                                 double               referenceValue,
 177                                 const int            numValuesInSum)
 178 {
 179     const double unit_factor = pull_conversion_factor_internal2userinput(&coordParams);
 180
 181     fprintf(out, "\t%g", pcrdData.value*unit_factor/numValuesInSum);
 182
 183     if (pullParams.bPrintRefValue && coordParams.eType != epullEXTERNAL)
 184     {
 185         fprintf(out, "\t%g", referenceValue*unit_factor/numValuesInSum);
 186     }
 187
 188     if (pullParams.bPrintComp)
 189     {
 190         pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr01, numValuesInSum);
 191         if (coordParams.ngroup >= 4)
 192         {
 193             pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr23, numValuesInSum);
 194         }
 195         if (coordParams.ngroup >= 6)
 196         {
 197             pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr45, numValuesInSum);
 198         }
 199     }
 200 }
 201
 202 static void pull_print_x(FILE *out, pull_t *pull, double t)
 203 {
 204     fprintf(out, "%.4f", t);
 205
 206     for (size_t c = 0; c < pull->coord.size(); c++)
 207     {
 208         const pull_coord_work_t     &pcrd           = pull->coord[c];
 209         int                          numValuesInSum = 1;
 210         const PullCoordinateHistory *pcrdHistory    = nullptr;
 211
 212         if (pull->bXOutAverage)
 213         {
 214             pcrdHistory = &pull->coordForceHistory->pullCoordinateSums[c];
 215
 216             numValuesInSum = pull->coordForceHistory->numValuesInXSum;
 217             pull_print_coord_dr(out, pull->params, pcrd.params,
 218                                 *pcrdHistory, pcrdHistory->valueRef,
 219                                 numValuesInSum);
 220         }
 221         else
 222         {
 223             pull_print_coord_dr(out, pull->params, pcrd.params,
 224                                 pcrd.spatialData, pcrd.value_ref,
 225                                 numValuesInSum);
 226         }
 227
 228         if (pull->params.bPrintCOM)
 229         {
 230             if (pcrd.params.eGeom == epullgCYL)
 231             {
 232                 for (int m = 0; m < DIM; m++)
 233                 {
 234                     if (pcrd.params.dim[m])
 235                     {
 236                         /* This equates to if (pull->bXOutAverage) */
 237                         if (pcrdHistory)
 238                         {
 239                             fprintf(out, "\t%g", pcrdHistory->dynaX[m] / numValuesInSum);
 240                         }
 241                         else
 242                         {
 243                             fprintf(out, "\t%g", pull->dyna[c].x[m]);
 244                         }
 245                     }
 246                 }
 247             }
 248             else
 249             {
 250                 for (int m = 0; m < DIM; m++)
 251                 {
 252                     if (pcrd.params.dim[m])
 253                     {
 254                         if (pull->bXOutAverage)
 255                         {
 256                             fprintf(out, "\t%g", pull->coordForceHistory->pullGroupSums[pcrd.params.group[0]].x[m] / numValuesInSum);
 257                         }
 258                         else
 259                         {
 260                             fprintf(out, "\t%g", pull->group[pcrd.params.group[0]].x[m]);
 261                         }
 262                     }
 263                 }
 264             }
 265             for (int g = 1; g < pcrd.params.ngroup; g++)
 266             {
 267                 for (int m = 0; m < DIM; m++)
 268                 {
 269                     if (pcrd.params.dim[m])
 270                     {
 271                         if (pull->bXOutAverage)
 272                         {
 273                             fprintf(out, "\t%g", pull->coordForceHistory->pullGroupSums[pcrd.params.group[g]].x[m] / numValuesInSum);
 274                         }
 275                         else
 276                         {
 277                             fprintf(out, "\t%g", pull->group[pcrd.params.group[g]].x[m]);
 278                         }
 279                     }
 280                 }
 281             }
 282         }
 283     }
 284     fprintf(out, "\n");
 285
 286     if (pull->bXOutAverage)
 287     {
 288         pullResetHistory(pull->coordForceHistory, TRUE, FALSE);
 289     }
 290 }
 291
 292 static void pull_print_f(FILE *out, const pull_t *pull, double t)
 293 {
 294     fprintf(out, "%.4f", t);
 295
 296     if (pull->bFOutAverage)
 297     {
 298         for (size_t c = 0; c < pull->coord.size(); c++)
 299         {
 300             fprintf(out, "\t%g", pull->coordForceHistory->pullCoordinateSums[c].scalarForce / pull->coordForceHistory->numValuesInFSum);
 301         }
 302     }
 303     else
 304     {
 305         for (const pull_coord_work_t &coord : pull->coord)
 306         {
 307             fprintf(out, "\t%g", coord.scalarForce);
 308         }
 309     }
 310     fprintf(out, "\n");
 311
 312     if (pull->bFOutAverage)
 313     {
 314         pullResetHistory(pull->coordForceHistory, FALSE, TRUE);
 315     }
 316 }
 317
 318 void pull_print_output(struct pull_t *pull, int64_t step, double time)
 319 {
 320     GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "pull_print_output called before all external pull potentials have been applied");
 321
 322     if (pull->params.nstxout != 0)
 323     {
 324         /* Do not update the average if the number of observations already equal (or are
 325          * higher than) what should be in each average output. This can happen when
 326          * appending to a file from a checkpoint, which would otherwise include the
 327          * last value twice.*/
 328         if (pull->bXOutAverage && !pull->coord.empty() && pull->coordForceHistory->numValuesInXSum < pull->params.nstxout)
 329         {
 330             addToPullxHistory(pull);
 331         }
 332         if (step % pull->params.nstxout == 0)
 333         {
 334             pull_print_x(pull->out_x, pull, time);
 335         }
 336     }
 337
 338     if (pull->params.nstfout != 0)
 339     {
 340         /* Do not update the average if the number of observations already equal (or are
 341          * higher than) what should be in each average output. This can happen when
 342          * appending to a file from a checkpoint, which would otherwise include the
 343          * last value twice.*/
 344         if (pull->bFOutAverage && !pull->coord.empty() && pull->coordForceHistory->numValuesInFSum < pull->params.nstfout)
 345         {
 346             addToPullfHistory(pull);
 347         }
 348         if (step % pull->params.nstfout == 0)
 349         {
 350             pull_print_f(pull->out_f, pull, time);
 351         }
 352     }
 353 }
 354
 355 static void set_legend_for_coord_components(const pull_coord_work_t *pcrd, int coord_index, char **setname, int *nsets_ptr)
 356 {
 357     /*  Loop over the distance vectors and print their components. Each vector is made up of two consecutive groups. */
 358     for (int g = 0; g < pcrd->params.ngroup; g += 2)
 359     {
 360         /* Loop over the components */
 361         for (int m  = 0; m < DIM; m++)
 362         {
 363             if (pcrd->params.dim[m])
 364             {
 365                 char legend[STRLEN];
 366
 367                 if (g == 0 && pcrd->params.ngroup <= 2)
 368                 {
 369                     /*  For the simplest case we print a simplified legend without group indices, just the cooordinate index
 370                         and which dimensional component it is. */
 371                     sprintf(legend, "%d d%c", coord_index + 1, 'X' + m);
 372                 }
 373                 else
 374                 {
 375                     /* Otherwise, print also the group indices (relative to the coordinate, not the global ones). */
 376                     sprintf(legend, "%d g %d-%d d%c", coord_index + 1, g + 1, g + 2, 'X' + m);
 377                 }
 378
 379                 setname[*nsets_ptr] = gmx_strdup(legend);
 380                 (*nsets_ptr)++;
 381             }
 382         }
 383     }
 384 }
 385
 386 static FILE *open_pull_out(const char *fn, struct pull_t *pull,
 387                            const gmx_output_env_t *oenv,
 388                            gmx_bool bCoord,
 389                            const gmx::ContinuationOptions &continuationOptions)
 390 {
 391     FILE  *fp;
 392     int    nsets, m;
 393     char **setname, buf[50];
 394
 395     if (continuationOptions.appendFiles)
 396     {
 397         fp = gmx_fio_fopen(fn, "a+");
 398     }
 399     else
 400     {
 401         fp = gmx_fio_fopen(fn, "w+");
 402         if (bCoord)
 403         {
 404             sprintf(buf, "Position (nm%s)", pull->bAngle ? ", deg" : "");
 405             if (pull->bXOutAverage)
 406             {
 407                 xvgr_header(fp, "Pull Average COM",  "Time (ps)", buf,
 408                             exvggtXNY, oenv);
 409             }
 410             else
 411             {
 412                 xvgr_header(fp, "Pull COM",  "Time (ps)", buf,
 413                             exvggtXNY, oenv);
 414             }
 415         }
 416         else
 417         {
 418             sprintf(buf, "Force (kJ/mol/nm%s)", pull->bAngle ? ", kJ/mol/rad" : "");
 419             if (pull->bFOutAverage)
 420             {
 421                 xvgr_header(fp, "Pull Average force", "Time (ps)", buf,
 422                             exvggtXNY, oenv);
 423             }
 424             else
 425             {
 426                 xvgr_header(fp, "Pull force", "Time (ps)", buf,
 427                             exvggtXNY, oenv);
 428             }
 429         }
 430
 431         /* With default mdp options only the actual coordinate value is printed (1),
 432          * but optionally the reference value (+ 1),
 433          * the group COMs for all the groups (+ ngroups_max*DIM)
 434          * and the components of the distance vectors can be printed (+ (ngroups_max/2)*DIM).
 435          */
 436         snew(setname, pull->coord.size()*(1 + 1 + c_pullCoordNgroupMax*DIM + c_pullCoordNgroupMax/2*DIM));
 437
 438         nsets = 0;
 439         for (size_t c = 0; c < pull->coord.size(); c++)
 440         {
 441             if (bCoord)
 442             {
 443                 /* The order of this legend should match the order of printing
 444                  * the data in print_pull_x above.
 445                  */
 446
 447                 /* The pull coord distance */
 448                 sprintf(buf, "%zu", c+1);
 449                 setname[nsets] = gmx_strdup(buf);
 450                 nsets++;
 451                 if (pull->params.bPrintRefValue &&
 452                     pull->coord[c].params.eType != epullEXTERNAL)
 453                 {
 454                     sprintf(buf, "%zu ref", c+1);
 455                     setname[nsets] = gmx_strdup(buf);
 456                     nsets++;
 457                 }
 458                 if (pull->params.bPrintComp)
 459                 {
 460                     set_legend_for_coord_components(&pull->coord[c], c, setname, &nsets);
 461                 }
 462
 463                 if (pull->params.bPrintCOM)
 464                 {
 465                     for (int g = 0; g < pull->coord[c].params.ngroup; g++)
 466                     {
 467                         /* Legend for reference group position */
 468                         for (m = 0; m < DIM; m++)
 469                         {
 470                             if (pull->coord[c].params.dim[m])
 471                             {
 472                                 sprintf(buf, "%zu g %d %c", c+1, g + 1, 'X'+m);
 473                                 setname[nsets] = gmx_strdup(buf);
 474                                 nsets++;
 475                             }
 476                         }
 477                     }
 478                 }
 479             }
 480             else
 481             {
 482                 /* For the pull force we always only use one scalar */
 483                 sprintf(buf, "%zu", c+1);
 484                 setname[nsets] = gmx_strdup(buf);
 485                 nsets++;
 486             }
 487         }
 488         if (nsets > 1)
 489         {
 490             xvgr_legend(fp, nsets, setname, oenv);
 491         }
 492         for (int c = 0; c < nsets; c++)
 493         {
 494             sfree(setname[c]);
 495         }
 496         sfree(setname);
 497     }
 498
 499     return fp;
 500 }
 501
 502 void init_pull_output_files(pull_t                         *pull,
 503                             int                             nfile,
 504                             const t_filenm                  fnm[],
 505                             const gmx_output_env_t         *oenv,
 506                             const gmx::ContinuationOptions &continuationOptions)
 507 {
 508     /* Check for px and pf filename collision, if we are writing
 509        both files */
 510     std::string px_filename, pf_filename;
 511     std::string px_appended, pf_appended;
 512     try
 513     {
 514         px_filename  = std::string(opt2fn("-px", nfile, fnm));
 515         pf_filename  = std::string(opt2fn("-pf", nfile, fnm));
 516     }
 517     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 518
 519     if ((pull->params.nstxout != 0) &&
 520         (pull->params.nstfout != 0) &&
 521         (px_filename == pf_filename))
 522     {
 523         if (!opt2bSet("-px", nfile, fnm) && !opt2bSet("-pf", nfile, fnm))
 524         {
 525             /* We are writing both pull files but neither set directly. */
 526             try
 527             {
 528                 px_appended   = append_before_extension(px_filename, "_pullx");
 529                 pf_appended   = append_before_extension(pf_filename, "_pullf");
 530             }
 531             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 532             pull->out_x = open_pull_out(px_appended.c_str(), pull, oenv,
 533                                         TRUE, continuationOptions);
 534             pull->out_f = open_pull_out(pf_appended.c_str(), pull, oenv,
 535                                         FALSE, continuationOptions);
 536             return;
 537         }
 538         else
 539         {
 540             /* If at least one of -px and -pf is set but the filenames are identical: */
 541             gmx_fatal(FARGS, "Identical pull_x and pull_f output filenames %s",
 542                       px_filename.c_str());
 543         }
 544     }
 545     if (pull->params.nstxout != 0)
 546     {
 547         pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
 548                                     TRUE, continuationOptions);
 549     }
 550     if (pull->params.nstfout != 0)
 551     {
 552         pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
 553                                     FALSE, continuationOptions);
 554     }
 555 }
 556
 557 void initPullHistory(pull_t             *pull,
 558                      ObservablesHistory *observablesHistory)
 559 {
 560     GMX_RELEASE_ASSERT(pull, "Need a valid pull object");
 561
 562     if (observablesHistory == nullptr)
 563     {
 564         pull->coordForceHistory = nullptr;
 565         return;
 566     }
 567     /* If pull->coordForceHistory is already set we are starting from a checkpoint. Do not reset it. */
 568     if (observablesHistory->pullHistory == nullptr)
 569     {
 570         observablesHistory->pullHistory           = std::make_unique<PullHistory>();
 571         pull->coordForceHistory                   = observablesHistory->pullHistory.get();
 572         pull->coordForceHistory->numValuesInXSum  = 0;
 573         pull->coordForceHistory->numValuesInFSum  = 0;
 574         pull->coordForceHistory->pullCoordinateSums.resize(pull->coord.size());
 575         pull->coordForceHistory->pullGroupSums.resize(pull->group.size());
 576     }
 577     else
 578     {
 579         pull->coordForceHistory                   = observablesHistory->pullHistory.get();
 580     }
 581 }