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(
 216                     out, pull->params, pcrd.params, *pcrdHistory, pcrdHistory->valueRef, numValuesInSum);
 217         }
 218         else
 219         {
 220             pull_print_coord_dr(
 221                     out, pull->params, pcrd.params, pcrd.spatialData, pcrd.value_ref, 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,
 253                                     "\t%g",
 254                                     pull->coordForceHistory->pullGroupSums[pcrd.params.group[0]].x[m]
 255                                             / numValuesInSum);
 256                         }
 257                         else
 258                         {
 259                             fprintf(out, "\t%g", pull->group[pcrd.params.group[0]].x[m]);
 260                         }
 261                     }
 262                 }
 263             }
 264             for (int g = 1; g < pcrd.params.ngroup; g++)
 265             {
 266                 for (int m = 0; m < DIM; m++)
 267                 {
 268                     if (pcrd.params.dim[m])
 269                     {
 270                         if (pull->bXOutAverage)
 271                         {
 272                             fprintf(out,
 273                                     "\t%g",
 274                                     pull->coordForceHistory->pullGroupSums[pcrd.params.group[g]].x[m]
 275                                             / numValuesInSum);
 276                         }
 277                         else
 278                         {
 279                             fprintf(out, "\t%g", pull->group[pcrd.params.group[g]].x[m]);
 280                         }
 281                     }
 282                 }
 283             }
 284         }
 285     }
 286     fprintf(out, "\n");
 287
 288     if (pull->bXOutAverage)
 289     {
 290         pullResetHistory(pull->coordForceHistory, TRUE, FALSE);
 291     }
 292 }
 293
 294 static void pull_print_f(FILE* out, const pull_t* pull, double t)
 295 {
 296     fprintf(out, "%.4f", t);
 297
 298     if (pull->bFOutAverage)
 299     {
 300         for (size_t c = 0; c < pull->coord.size(); c++)
 301         {
 302             fprintf(out,
 303                     "\t%g",
 304                     pull->coordForceHistory->pullCoordinateSums[c].scalarForce
 305                             / pull->coordForceHistory->numValuesInFSum);
 306         }
 307     }
 308     else
 309     {
 310         for (const pull_coord_work_t& coord : pull->coord)
 311         {
 312             fprintf(out, "\t%g", coord.scalarForce);
 313         }
 314     }
 315     fprintf(out, "\n");
 316
 317     if (pull->bFOutAverage)
 318     {
 319         pullResetHistory(pull->coordForceHistory, FALSE, TRUE);
 320     }
 321 }
 322
 323 void pull_print_output(struct pull_t* pull, int64_t step, double time)
 324 {
 325     GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
 326                "pull_print_output called before all external pull potentials have been applied");
 327
 328     if (pull->params.nstxout != 0)
 329     {
 330         /* Do not update the average if the number of observations already equal (or are
 331          * higher than) what should be in each average output. This can happen when
 332          * appending to a file from a checkpoint, which would otherwise include the
 333          * last value twice.*/
 334         if (pull->bXOutAverage && !pull->coord.empty()
 335             && pull->coordForceHistory->numValuesInXSum < pull->params.nstxout)
 336         {
 337             addToPullxHistory(pull);
 338         }
 339         if (step % pull->params.nstxout == 0)
 340         {
 341             pull_print_x(pull->out_x, pull, time);
 342         }
 343     }
 344
 345     if (pull->params.nstfout != 0)
 346     {
 347         /* Do not update the average if the number of observations already equal (or are
 348          * higher than) what should be in each average output. This can happen when
 349          * appending to a file from a checkpoint, which would otherwise include the
 350          * last value twice.*/
 351         if (pull->bFOutAverage && !pull->coord.empty()
 352             && pull->coordForceHistory->numValuesInFSum < pull->params.nstfout)
 353         {
 354             addToPullfHistory(pull);
 355         }
 356         if (step % pull->params.nstfout == 0)
 357         {
 358             pull_print_f(pull->out_f, pull, time);
 359         }
 360     }
 361 }
 362
 363 static void set_legend_for_coord_components(const pull_coord_work_t* pcrd,
 364                                             int                      coord_index,
 365                                             char**                   setname,
 366                                             int*                     nsets_ptr)
 367 {
 368     /*  Loop over the distance vectors and print their components. Each vector is made up of two consecutive groups. */
 369     for (int g = 0; g < pcrd->params.ngroup; g += 2)
 370     {
 371         /* Loop over the components */
 372         for (int m = 0; m < DIM; m++)
 373         {
 374             if (pcrd->params.dim[m])
 375             {
 376                 char legend[STRLEN];
 377
 378                 if (g == 0 && pcrd->params.ngroup <= 2)
 379                 {
 380                     /*  For the simplest case we print a simplified legend without group indices,
 381                        just the cooordinate index and which dimensional component it is. */
 382                     sprintf(legend, "%d d%c", coord_index + 1, 'X' + m);
 383                 }
 384                 else
 385                 {
 386                     /* Otherwise, print also the group indices (relative to the coordinate, not the global ones). */
 387                     sprintf(legend, "%d g %d-%d d%c", coord_index + 1, g + 1, g + 2, 'X' + m);
 388                 }
 389
 390                 setname[*nsets_ptr] = gmx_strdup(legend);
 391                 (*nsets_ptr)++;
 392             }
 393         }
 394     }
 395 }
 396
 397 static FILE* open_pull_out(const char*             fn,
 398                            struct pull_t*          pull,
 399                            const gmx_output_env_t* oenv,
 400                            gmx_bool                bCoord,
 401                            const bool              restartWithAppending)
 402 {
 403     FILE*  fp;
 404     int    nsets, m;
 405     char **setname, buf[50];
 406
 407     if (restartWithAppending)
 408     {
 409         fp = gmx_fio_fopen(fn, "a+");
 410     }
 411     else
 412     {
 413         fp = gmx_fio_fopen(fn, "w+");
 414         if (bCoord)
 415         {
 416             sprintf(buf, "Position (nm%s)", pull->bAngle ? ", deg" : "");
 417             if (pull->bXOutAverage)
 418             {
 419                 xvgr_header(fp, "Pull Average COM", "Time (ps)", buf, exvggtXNY, oenv);
 420             }
 421             else
 422             {
 423                 xvgr_header(fp, "Pull COM", "Time (ps)", buf, exvggtXNY, oenv);
 424             }
 425         }
 426         else
 427         {
 428             sprintf(buf, "Force (kJ/mol/nm%s)", pull->bAngle ? ", kJ/mol/rad" : "");
 429             if (pull->bFOutAverage)
 430             {
 431                 xvgr_header(fp, "Pull Average force", "Time (ps)", buf, exvggtXNY, oenv);
 432             }
 433             else
 434             {
 435                 xvgr_header(fp, "Pull force", "Time (ps)", buf, exvggtXNY, oenv);
 436             }
 437         }
 438
 439         /* With default mdp options only the actual coordinate value is printed (1),
 440          * but optionally the reference value (+ 1),
 441          * the group COMs for all the groups (+ ngroups_max*DIM)
 442          * and the components of the distance vectors can be printed (+ (ngroups_max/2)*DIM).
 443          */
 444         snew(setname,
 445              pull->coord.size() * (1 + 1 + c_pullCoordNgroupMax * DIM + c_pullCoordNgroupMax / 2 * DIM));
 446
 447         nsets = 0;
 448         for (size_t c = 0; c < pull->coord.size(); c++)
 449         {
 450             if (bCoord)
 451             {
 452                 /* The order of this legend should match the order of printing
 453                  * the data in print_pull_x above.
 454                  */
 455
 456                 /* The pull coord distance */
 457                 sprintf(buf, "%zu", c + 1);
 458                 setname[nsets] = gmx_strdup(buf);
 459                 nsets++;
 460                 if (pull->params.bPrintRefValue && pull->coord[c].params.eType != epullEXTERNAL)
 461                 {
 462                     sprintf(buf, "%zu ref", c + 1);
 463                     setname[nsets] = gmx_strdup(buf);
 464                     nsets++;
 465                 }
 466                 if (pull->params.bPrintComp)
 467                 {
 468                     set_legend_for_coord_components(&pull->coord[c], c, setname, &nsets);
 469                 }
 470
 471                 if (pull->params.bPrintCOM)
 472                 {
 473                     for (int g = 0; g < pull->coord[c].params.ngroup; g++)
 474                     {
 475                         /* Legend for reference group position */
 476                         for (m = 0; m < DIM; m++)
 477                         {
 478                             if (pull->coord[c].params.dim[m])
 479                             {
 480                                 sprintf(buf, "%zu g %d %c", c + 1, g + 1, 'X' + m);
 481                                 setname[nsets] = gmx_strdup(buf);
 482                                 nsets++;
 483                             }
 484                         }
 485                     }
 486                 }
 487             }
 488             else
 489             {
 490                 /* For the pull force we always only use one scalar */
 491                 sprintf(buf, "%zu", c + 1);
 492                 setname[nsets] = gmx_strdup(buf);
 493                 nsets++;
 494             }
 495         }
 496         if (nsets > 1)
 497         {
 498             xvgr_legend(fp, nsets, setname, oenv);
 499         }
 500         for (int c = 0; c < nsets; c++)
 501         {
 502             sfree(setname[c]);
 503         }
 504         sfree(setname);
 505     }
 506
 507     return fp;
 508 }
 509
 510 void init_pull_output_files(pull_t*                     pull,
 511                             int                         nfile,
 512                             const t_filenm              fnm[],
 513                             const gmx_output_env_t*     oenv,
 514                             const gmx::StartingBehavior startingBehavior)
 515 {
 516     /* Check for px and pf filename collision, if we are writing
 517        both files */
 518     std::string px_filename, pf_filename;
 519     std::string px_appended, pf_appended;
 520     try
 521     {
 522         px_filename = std::string(opt2fn("-px", nfile, fnm));
 523         pf_filename = std::string(opt2fn("-pf", nfile, fnm));
 524     }
 525     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 526
 527     bool restartWithAppending = startingBehavior == gmx::StartingBehavior::RestartWithAppending;
 528     if ((pull->params.nstxout != 0) && (pull->params.nstfout != 0) && (px_filename == pf_filename))
 529     {
 530         if (!opt2bSet("-px", nfile, fnm) && !opt2bSet("-pf", nfile, fnm))
 531         {
 532             /* We are writing both pull files but neither set directly. */
 533             try
 534             {
 535                 px_appended = append_before_extension(px_filename, "_pullx");
 536                 pf_appended = append_before_extension(pf_filename, "_pullf");
 537             }
 538             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 539             pull->out_x = open_pull_out(px_appended.c_str(), pull, oenv, TRUE, restartWithAppending);
 540             pull->out_f = open_pull_out(pf_appended.c_str(), pull, oenv, FALSE, restartWithAppending);
 541             return;
 542         }
 543         else
 544         {
 545             /* If at least one of -px and -pf is set but the filenames are identical: */
 546             gmx_fatal(FARGS, "Identical pull_x and pull_f output filenames %s", px_filename.c_str());
 547         }
 548     }
 549     if (pull->params.nstxout != 0)
 550     {
 551         pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, restartWithAppending);
 552     }
 553     if (pull->params.nstfout != 0)
 554     {
 555         pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv, FALSE, restartWithAppending);
 556     }
 557 }
 558
 559 void initPullHistory(pull_t* pull, ObservablesHistory* observablesHistory)
 560 {
 561     GMX_RELEASE_ASSERT(pull, "Need a valid pull object");
 562
 563     if (observablesHistory == nullptr)
 564     {
 565         pull->coordForceHistory = nullptr;
 566         return;
 567     }
 568     /* If pull->coordForceHistory is already set we are starting from a checkpoint. Do not reset it. */
 569     if (observablesHistory->pullHistory == nullptr)
 570     {
 571         observablesHistory->pullHistory          = std::make_unique<PullHistory>();
 572         pull->coordForceHistory                  = observablesHistory->pullHistory.get();
 573         pull->coordForceHistory->numValuesInXSum = 0;
 574         pull->coordForceHistory->numValuesInFSum = 0;
 575         pull->coordForceHistory->pullCoordinateSums.resize(pull->coord.size());
 576         pull->coordForceHistory->pullGroupSums.resize(pull->group.size());
 577     }
 578     else
 579     {
 580         pull->coordForceHistory = observablesHistory->pullHistory.get();
 581     }
 582 }