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