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