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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "output.h"
  41
  42 #include <cstdio>
  43
  44 #include <memory>
  45
  46 #include "gromacs/commandline/filenm.h"
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/xvgr.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/mdrunutility/handlerestart.h"
  51 #include "gromacs/mdtypes/mdrunoptions.h"
  52 #include "gromacs/mdtypes/observableshistory.h"
  53 #include "gromacs/mdtypes/pullhistory.h"
  54 #include "gromacs/pulling/pull.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/exceptions.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/gmxassert.h"
  59
  60 #include "pull_internal.h"
  61
  62 static std::string append_before_extension(const std::string& pathname, const std::string& to_append)
  63 {
  64     /* Appends to_append before last '.' in pathname */
  65     size_t extPos = pathname.find_last_of('.');
  66     if (extPos == std::string::npos)
  67     {
  68         return pathname + to_append;
  69     }
  70     else
  71     {
  72         return pathname.substr(0, extPos) + to_append + 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
 120         PullCoordinateHistory& pcrdHistory = pull->coordForceHistory->pullCoordinateSums[c];
 121
 122         pcrdHistory.scalarForce += pcrd.scalarForce;
 123     }
 124 }
 125
 126 static void pullResetHistory(PullHistory* history, bool resetXHistory, bool resetFHistory)
 127 {
 128     if (resetXHistory)
 129     {
 130         history->numValuesInXSum = 0;
 131
 132         for (PullCoordinateHistory& pcrdHistory : history->pullCoordinateSums)
 133         {
 134             pcrdHistory.value    = 0;
 135             pcrdHistory.valueRef = 0;
 136
 137             clear_dvec(pcrdHistory.dr01);
 138             clear_dvec(pcrdHistory.dr23);
 139             clear_dvec(pcrdHistory.dr45);
 140             clear_dvec(pcrdHistory.dynaX);
 141         }
 142
 143         for (PullGroupHistory& pgrpHistory : history->pullGroupSums)
 144         {
 145             clear_dvec(pgrpHistory.x);
 146         }
 147     }
 148     if (resetFHistory)
 149     {
 150         history->numValuesInFSum = 0;
 151         for (PullCoordinateHistory& pcrdHistory : history->pullCoordinateSums)
 152         {
 153             pcrdHistory.scalarForce = 0;
 154         }
 155     }
 156 }
 157
 158 static void pull_print_coord_dr_components(FILE* out, const ivec dim, const dvec dr, const int numValuesInSum)
 159 {
 160     for (int m = 0; m < DIM; m++)
 161     {
 162         if (dim[m])
 163         {
 164             fprintf(out, "\t%g", dr[m] / numValuesInSum);
 165         }
 166     }
 167 }
 168
 169 template<typename T>
 170 static void pull_print_coord_dr(FILE*                out,
 171                                 const pull_params_t& pullParams,
 172                                 const t_pull_coord&  coordParams,
 173                                 const T&             pcrdData,
 174                                 double               referenceValue,
 175                                 const int            numValuesInSum)
 176 {
 177     const double unit_factor = pull_conversion_factor_internal2userinput(&coordParams);
 178
 179     fprintf(out, "\t%g", pcrdData.value * unit_factor / numValuesInSum);
 180
 181     if (pullParams.bPrintRefValue && coordParams.eType != epullEXTERNAL)
 182     {
 183         fprintf(out, "\t%g", referenceValue * unit_factor / numValuesInSum);
 184     }
 185
 186     if (pullParams.bPrintComp)
 187     {
 188         pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr01, numValuesInSum);
 189         if (coordParams.ngroup >= 4)
 190         {
 191             pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr23, numValuesInSum);
 192         }
 193         if (coordParams.ngroup >= 6)
 194         {
 195             pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr45, numValuesInSum);
 196         }
 197     }
 198 }
 199
 200 static void pull_print_x(FILE* out, pull_t* pull, double t)
 201 {
 202     fprintf(out, "%.4f", t);
 203
 204     for (size_t c = 0; c < pull->coord.size(); c++)
 205     {
 206         const pull_coord_work_t&     pcrd           = pull->coord[c];
 207         int                          numValuesInSum = 1;
 208         const PullCoordinateHistory* pcrdHistory    = nullptr;
 209
 210         if (pull->bXOutAverage)
 211         {
 212             pcrdHistory = &pull->coordForceHistory->pullCoordinateSums[c];
 213
 214             numValuesInSum = pull->coordForceHistory->numValuesInXSum;
 215             pull_print_coord_dr(out, pull->params, pcrd.params, *pcrdHistory, pcrdHistory->valueRef,
 216                                 numValuesInSum);
 217         }
 218         else
 219         {
 220             pull_print_coord_dr(out, pull->params, pcrd.params, pcrd.spatialData, pcrd.value_ref,
 221                                 numValuesInSum);
 222         }
 223
 224         if (pull->params.bPrintCOM)
 225         {
 226             if (pcrd.params.eGeom == epullgCYL)
 227             {
 228                 for (int m = 0; m < DIM; m++)
 229                 {
 230                     if (pcrd.params.dim[m])
 231                     {
 232                         /* This equates to if (pull->bXOutAverage) */
 233                         if (pcrdHistory)
 234                         {
 235                             fprintf(out, "\t%g", pcrdHistory->dynaX[m] / numValuesInSum);
 236                         }
 237                         else
 238                         {
 239                             fprintf(out, "\t%g", pull->dyna[c].x[m]);
 240                         }
 241                     }
 242                 }
 243             }
 244             else
 245             {
 246                 for (int m = 0; m < DIM; m++)
 247                 {
 248                     if (pcrd.params.dim[m])
 249                     {
 250                         if (pull->bXOutAverage)
 251                         {
 252                             fprintf(out, "\t%g",
 253                                     pull->coordForceHistory->pullGroupSums[pcrd.params.group[0]].x[m]
 254                                             / numValuesInSum);
 255                         }
 256                         else
 257                         {
 258                             fprintf(out, "\t%g", pull->group[pcrd.params.group[0]].x[m]);
 259                         }
 260                     }
 261                 }
 262             }
 263             for (int g = 1; g < pcrd.params.ngroup; g++)
 264             {
 265                 for (int m = 0; m < DIM; m++)
 266                 {
 267                     if (pcrd.params.dim[m])
 268                     {
 269                         if (pull->bXOutAverage)
 270                         {
 271                             fprintf(out, "\t%g",
 272                                     pull->coordForceHistory->pullGroupSums[pcrd.params.group[g]].x[m]
 273                                             / 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",
 301                     pull->coordForceHistory->pullCoordinateSums[c].scalarForce
 302                             / pull->coordForceHistory->numValuesInFSum);
 303         }
 304     }
 305     else
 306     {
 307         for (const pull_coord_work_t& coord : pull->coord)
 308         {
 309             fprintf(out, "\t%g", coord.scalarForce);
 310         }
 311     }
 312     fprintf(out, "\n");
 313
 314     if (pull->bFOutAverage)
 315     {
 316         pullResetHistory(pull->coordForceHistory, FALSE, TRUE);
 317     }
 318 }
 319
 320 void pull_print_output(struct pull_t* pull, int64_t step, double time)
 321 {
 322     GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
 323                "pull_print_output called before all external pull potentials have been applied");
 324
 325     if (pull->params.nstxout != 0)
 326     {
 327         /* Do not update the average if the number of observations already equal (or are
 328          * higher than) what should be in each average output. This can happen when
 329          * appending to a file from a checkpoint, which would otherwise include the
 330          * last value twice.*/
 331         if (pull->bXOutAverage && !pull->coord.empty()
 332             && pull->coordForceHistory->numValuesInXSum < pull->params.nstxout)
 333         {
 334             addToPullxHistory(pull);
 335         }
 336         if (step % pull->params.nstxout == 0)
 337         {
 338             pull_print_x(pull->out_x, pull, time);
 339         }
 340     }
 341
 342     if (pull->params.nstfout != 0)
 343     {
 344         /* Do not update the average if the number of observations already equal (or are
 345          * higher than) what should be in each average output. This can happen when
 346          * appending to a file from a checkpoint, which would otherwise include the
 347          * last value twice.*/
 348         if (pull->bFOutAverage && !pull->coord.empty()
 349             && pull->coordForceHistory->numValuesInFSum < pull->params.nstfout)
 350         {
 351             addToPullfHistory(pull);
 352         }
 353         if (step % pull->params.nstfout == 0)
 354         {
 355             pull_print_f(pull->out_f, pull, time);
 356         }
 357     }
 358 }
 359
 360 static void set_legend_for_coord_components(const pull_coord_work_t* pcrd,
 361                                             int                      coord_index,
 362                                             char**                   setname,
 363                                             int*                     nsets_ptr)
 364 {
 365     /*  Loop over the distance vectors and print their components. Each vector is made up of two consecutive groups. */
 366     for (int g = 0; g < pcrd->params.ngroup; g += 2)
 367     {
 368         /* Loop over the components */
 369         for (int m = 0; m < DIM; m++)
 370         {
 371             if (pcrd->params.dim[m])
 372             {
 373                 char legend[STRLEN];
 374
 375                 if (g == 0 && pcrd->params.ngroup <= 2)
 376                 {
 377                     /*  For the simplest case we print a simplified legend without group indices,
 378                        just the cooordinate index and which dimensional component it is. */
 379                     sprintf(legend, "%d d%c", coord_index + 1, 'X' + m);
 380                 }
 381                 else
 382                 {
 383                     /* Otherwise, print also the group indices (relative to the coordinate, not the global ones). */
 384                     sprintf(legend, "%d g %d-%d d%c", coord_index + 1, g + 1, g + 2, 'X' + m);
 385                 }
 386
 387                 setname[*nsets_ptr] = gmx_strdup(legend);
 388                 (*nsets_ptr)++;
 389             }
 390         }
 391     }
 392 }
 393
 394 static FILE* open_pull_out(const char*             fn,
 395                            struct pull_t*          pull,
 396                            const gmx_output_env_t* oenv,
 397                            gmx_bool                bCoord,
 398                            const bool              restartWithAppending)
 399 {
 400     FILE*  fp;
 401     int    nsets, m;
 402     char **setname, buf[50];
 403
 404     if (restartWithAppending)
 405     {
 406         fp = gmx_fio_fopen(fn, "a+");
 407     }
 408     else
 409     {
 410         fp = gmx_fio_fopen(fn, "w+");
 411         if (bCoord)
 412         {
 413             sprintf(buf, "Position (nm%s)", pull->bAngle ? ", deg" : "");
 414             if (pull->bXOutAverage)
 415             {
 416                 xvgr_header(fp, "Pull Average COM", "Time (ps)", buf, exvggtXNY, oenv);
 417             }
 418             else
 419             {
 420                 xvgr_header(fp, "Pull COM", "Time (ps)", buf, exvggtXNY, oenv);
 421             }
 422         }
 423         else
 424         {
 425             sprintf(buf, "Force (kJ/mol/nm%s)", pull->bAngle ? ", kJ/mol/rad" : "");
 426             if (pull->bFOutAverage)
 427             {
 428                 xvgr_header(fp, "Pull Average force", "Time (ps)", buf, exvggtXNY, oenv);
 429             }
 430             else
 431             {
 432                 xvgr_header(fp, "Pull force", "Time (ps)", buf, exvggtXNY, oenv);
 433             }
 434         }
 435
 436         /* With default mdp options only the actual coordinate value is printed (1),
 437          * but optionally the reference value (+ 1),
 438          * the group COMs for all the groups (+ ngroups_max*DIM)
 439          * and the components of the distance vectors can be printed (+ (ngroups_max/2)*DIM).
 440          */
 441         snew(setname, pull->coord.size()
 442                               * (1 + 1 + c_pullCoordNgroupMax * DIM + c_pullCoordNgroupMax / 2 * DIM));
 443
 444         nsets = 0;
 445         for (size_t c = 0; c < pull->coord.size(); c++)
 446         {
 447             if (bCoord)
 448             {
 449                 /* The order of this legend should match the order of printing
 450                  * the data in print_pull_x above.
 451                  */
 452
 453                 /* The pull coord distance */
 454                 sprintf(buf, "%zu", c + 1);
 455                 setname[nsets] = gmx_strdup(buf);
 456                 nsets++;
 457                 if (pull->params.bPrintRefValue && pull->coord[c].params.eType != epullEXTERNAL)
 458                 {
 459                     sprintf(buf, "%zu ref", c + 1);
 460                     setname[nsets] = gmx_strdup(buf);
 461                     nsets++;
 462                 }
 463                 if (pull->params.bPrintComp)
 464                 {
 465                     set_legend_for_coord_components(&pull->coord[c], c, setname, &nsets);
 466                 }
 467
 468                 if (pull->params.bPrintCOM)
 469                 {
 470                     for (int g = 0; g < pull->coord[c].params.ngroup; g++)
 471                     {
 472                         /* Legend for reference group position */
 473                         for (m = 0; m < DIM; m++)
 474                         {
 475                             if (pull->coord[c].params.dim[m])
 476                             {
 477                                 sprintf(buf, "%zu g %d %c", c + 1, g + 1, 'X' + m);
 478                                 setname[nsets] = gmx_strdup(buf);
 479                                 nsets++;
 480                             }
 481                         }
 482                     }
 483                 }
 484             }
 485             else
 486             {
 487                 /* For the pull force we always only use one scalar */
 488                 sprintf(buf, "%zu", c + 1);
 489                 setname[nsets] = gmx_strdup(buf);
 490                 nsets++;
 491             }
 492         }
 493         if (nsets > 1)
 494         {
 495             xvgr_legend(fp, nsets, setname, oenv);
 496         }
 497         for (int c = 0; c < nsets; c++)
 498         {
 499             sfree(setname[c]);
 500         }
 501         sfree(setname);
 502     }
 503
 504     return fp;
 505 }
 506
 507 void init_pull_output_files(pull_t*                     pull,
 508                             int                         nfile,
 509                             const t_filenm              fnm[],
 510                             const gmx_output_env_t*     oenv,
 511                             const gmx::StartingBehavior startingBehavior)
 512 {
 513     /* Check for px and pf filename collision, if we are writing
 514        both files */
 515     std::string px_filename, pf_filename;
 516     std::string px_appended, pf_appended;
 517     try
 518     {
 519         px_filename = std::string(opt2fn("-px", nfile, fnm));
 520         pf_filename = std::string(opt2fn("-pf", nfile, fnm));
 521     }
 522     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 523
 524     bool restartWithAppending = startingBehavior == gmx::StartingBehavior::RestartWithAppending;
 525     if ((pull->params.nstxout != 0) && (pull->params.nstfout != 0) && (px_filename == pf_filename))
 526     {
 527         if (!opt2bSet("-px", nfile, fnm) && !opt2bSet("-pf", nfile, fnm))
 528         {
 529             /* We are writing both pull files but neither set directly. */
 530             try
 531             {
 532                 px_appended = append_before_extension(px_filename, "_pullx");
 533                 pf_appended = append_before_extension(pf_filename, "_pullf");
 534             }
 535             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 536             pull->out_x = open_pull_out(px_appended.c_str(), pull, oenv, TRUE, restartWithAppending);
 537             pull->out_f = open_pull_out(pf_appended.c_str(), pull, oenv, FALSE, restartWithAppending);
 538             return;
 539         }
 540         else
 541         {
 542             /* If at least one of -px and -pf is set but the filenames are identical: */
 543             gmx_fatal(FARGS, "Identical pull_x and pull_f output filenames %s", px_filename.c_str());
 544         }
 545     }
 546     if (pull->params.nstxout != 0)
 547     {
 548         pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, restartWithAppending);
 549     }
 550     if (pull->params.nstfout != 0)
 551     {
 552         pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv, FALSE, restartWithAppending);
 553     }
 554 }
 555
 556 void initPullHistory(pull_t* pull, ObservablesHistory* observablesHistory)
 557 {
 558     GMX_RELEASE_ASSERT(pull, "Need a valid pull object");
 559
 560     if (observablesHistory == nullptr)
 561     {
 562         pull->coordForceHistory = nullptr;
 563         return;
 564     }
 565     /* If pull->coordForceHistory is already set we are starting from a checkpoint. Do not reset it. */
 566     if (observablesHistory->pullHistory == nullptr)
 567     {
 568         observablesHistory->pullHistory          = std::make_unique<PullHistory>();
 569         pull->coordForceHistory                  = observablesHistory->pullHistory.get();
 570         pull->coordForceHistory->numValuesInXSum = 0;
 571         pull->coordForceHistory->numValuesInFSum = 0;
 572         pull->coordForceHistory->pullCoordinateSums.resize(pull->coord.size());
 573         pull->coordForceHistory->pullGroupSums.resize(pull->group.size());
 574     }
 575     else
 576     {
 577         pull->coordForceHistory = observablesHistory->pullHistory.get();
 578     }
 579 }