src/gromacs/gmxana/gmx_genion.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "copyrite.h"
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "sysstuff.h"
  44 #include "confio.h"
  45 #include "statutil.h"
  46 #include "pbc.h"
  47 #include "force.h"
  48 #include "gmx_fatal.h"
  49 #include "futil.h"
  50 #include "maths.h"
  51 #include "macros.h"
  52 #include "physics.h"
  53 #include "vec.h"
  54 #include "tpxio.h"
  55 #include "mdrun.h"
  56 #include "calcpot.h"
  57 #include "main.h"
  58 #include "random.h"
  59 #include "index.h"
  60 #include "mtop_util.h"
  61 #include "gmx_ana.h"
  62
  63 static int greatest_common_divisor(int p, int q)
  64 {
  65     int tmp;
  66     while (q != 0)
  67     {
  68         tmp = q;
  69         q = p % q;
  70         p = tmp;
  71     }
  72     return p;
  73 }
  74
  75 static void insert_ion(int nsa, int *nwater,
  76                        gmx_bool bSet[], int repl[], atom_id index[],
  77                        real pot[], rvec x[], t_pbc *pbc,
  78                        int sign, int q, const char *ionname,
  79                        t_mdatoms *mdatoms,
  80                        real rmin, gmx_bool bRandom, int *seed)
  81 {
  82     int             i, ii, ei, owater, wlast, m, nw;
  83     real            extr_e, poti, rmin2;
  84     rvec            xei, dx;
  85     gmx_bool        bSub = FALSE;
  86     gmx_large_int_t maxrand;
  87
  88     ei       = -1;
  89     nw       = *nwater;
  90     maxrand  = nw;
  91     maxrand *= 1000;
  92     if (bRandom)
  93     {
  94         do
  95         {
  96             ei = nw*rando(seed);
  97             maxrand--;
  98         }
  99         while (bSet[ei] && (maxrand > 0));
 100         if (bSet[ei])
 101         {
 102             gmx_fatal(FARGS, "No more replaceable solvent!");
 103         }
 104     }
 105     else
 106     {
 107         extr_e = 0;
 108         for (i = 0; (i < nw); i++)
 109         {
 110             if (!bSet[i])
 111             {
 112                 ii   = index[nsa*i];
 113                 poti = pot[ii];
 114                 if (q > 0)
 115                 {
 116                     if ((poti <= extr_e) || !bSub)
 117                     {
 118                         extr_e = poti;
 119                         ei     = i;
 120                         bSub   = TRUE;
 121                     }
 122                 }
 123                 else
 124                 {
 125                     if ((poti >= extr_e) || !bSub)
 126                     {
 127                         extr_e = poti;
 128                         ei     = i;
 129                         bSub   = TRUE;
 130                     }
 131                 }
 132             }
 133         }
 134         if (ei == -1)
 135         {
 136             gmx_fatal(FARGS, "No more replaceable solvent!");
 137         }
 138     }
 139     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
 140             ei, index[nsa*ei], ionname);
 141
 142     /* Replace solvent molecule charges with ion charge */
 143     bSet[ei] = TRUE;
 144     repl[ei] = sign;
 145     mdatoms->chargeA[index[nsa*ei]] = q;
 146     for (i = 1; i < nsa; i++)
 147     {
 148         mdatoms->chargeA[index[nsa*ei+i]] = 0;
 149     }
 150
 151     /* Mark all solvent molecules within rmin as unavailable for substitution */
 152     if (rmin > 0)
 153     {
 154         rmin2 = rmin*rmin;
 155         for (i = 0; (i < nw); i++)
 156         {
 157             if (!bSet[i])
 158             {
 159                 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
 160                 if (iprod(dx, dx) < rmin2)
 161                 {
 162                     bSet[i] = TRUE;
 163                 }
 164             }
 165         }
 166     }
 167 }
 168
 169 static char *aname(const char *mname)
 170 {
 171     char *str;
 172     int   i;
 173
 174     str = strdup(mname);
 175     i   = strlen(str)-1;
 176     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 177     {
 178         str[i] = '\0';
 179         i--;
 180     }
 181
 182     return str;
 183 }
 184
 185 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
 186                t_atoms *atoms, rvec x[],
 187                const char *p_name, const char *n_name)
 188 {
 189     int    i, j, k, r, np, nn, starta, startr, npi, nni;
 190     rvec  *xt;
 191     char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
 192
 193     snew(xt, atoms->nr);
 194
 195     /* Put all the solvent in front and count the added ions */
 196     np = 0;
 197     nn = 0;
 198     j  = index[0];
 199     for (i = 0; i < nw; i++)
 200     {
 201         r = repl[i];
 202         if (r == 0)
 203         {
 204             for (k = 0; k < nsa; k++)
 205             {
 206                 copy_rvec(x[index[nsa*i+k]], xt[j++]);
 207             }
 208         }
 209         else if (r > 0)
 210         {
 211             np++;
 212         }
 213         else if (r < 0)
 214         {
 215             nn++;
 216         }
 217     }
 218
 219     if (np+nn > 0)
 220     {
 221         /* Put the positive and negative ions at the end */
 222         starta = index[nsa*(nw - np - nn)];
 223         startr = atoms->atom[starta].resind;
 224
 225         if (np)
 226         {
 227             snew(pptr, 1);
 228             pptr[0] = strdup(p_name);
 229             snew(paptr, 1);
 230             paptr[0] = aname(p_name);
 231         }
 232         if (nn)
 233         {
 234             snew(nptr, 1);
 235             nptr[0] = strdup(n_name);
 236             snew(naptr, 1);
 237             naptr[0] = aname(n_name);
 238         }
 239         npi = 0;
 240         nni = 0;
 241         for (i = 0; i < nw; i++)
 242         {
 243             r = repl[i];
 244             if (r > 0)
 245             {
 246                 j = starta+npi;
 247                 k = startr+npi;
 248                 copy_rvec(x[index[nsa*i]], xt[j]);
 249                 atoms->atomname[j]     = paptr;
 250                 atoms->atom[j].resind  = k;
 251                 atoms->resinfo[k].name = pptr;
 252                 npi++;
 253             }
 254             else if (r < 0)
 255             {
 256                 j = starta+np+nni;
 257                 k = startr+np+nni;
 258                 copy_rvec(x[index[nsa*i]], xt[j]);
 259                 atoms->atomname[j]     = naptr;
 260                 atoms->atom[j].resind  = k;
 261                 atoms->resinfo[k].name = nptr;
 262                 nni++;
 263             }
 264         }
 265         for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
 266         {
 267             j                  = i-(nsa-1)*(np+nn);
 268             atoms->atomname[j] = atoms->atomname[i];
 269             atoms->atom[j]     = atoms->atom[i];
 270             copy_rvec(x[i], xt[j]);
 271         }
 272         atoms->nr -= (nsa-1)*(np+nn);
 273
 274         /* Copy the new positions back */
 275         for (i = index[0]; i < atoms->nr; i++)
 276         {
 277             copy_rvec(xt[i], x[i]);
 278         }
 279         sfree(xt);
 280     }
 281 }
 282
 283 static void update_topol(const char *topinout, int p_num, int n_num,
 284                          const char *p_name, const char *n_name, char *grpname)
 285 {
 286 #define TEMP_FILENM "temp.top"
 287     FILE    *fpin, *fpout;
 288     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 289     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 290     gmx_bool bMolecules;
 291
 292     printf("\nProcessing topology\n");
 293     fpin  = ffopen(topinout, "r");
 294     fpout = ffopen(TEMP_FILENM, "w");
 295
 296     line       = 0;
 297     bMolecules = FALSE;
 298     nmol_line  = 0;
 299     sol_line   = -1;
 300     nsol_last  = -1;
 301     while (fgets(buf, STRLEN, fpin))
 302     {
 303         line++;
 304         strcpy(buf2, buf);
 305         if ((temp = strchr(buf2, '\n')) != NULL)
 306         {
 307             temp[0] = '\0';
 308         }
 309         ltrim(buf2);
 310         if (buf2[0] == '[')
 311         {
 312             buf2[0] = ' ';
 313             if ((temp = strchr(buf2, '\n')) != NULL)
 314             {
 315                 temp[0] = '\0';
 316             }
 317             rtrim(buf2);
 318             if (buf2[strlen(buf2)-1] == ']')
 319             {
 320                 buf2[strlen(buf2)-1] = '\0';
 321                 ltrim(buf2);
 322                 rtrim(buf2);
 323                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 324             }
 325             fprintf(fpout, "%s", buf);
 326         }
 327         else if (!bMolecules)
 328         {
 329             fprintf(fpout, "%s", buf);
 330         }
 331         else
 332         {
 333             /* Check if this is a line with solvent molecules */
 334             sscanf(buf, "%s", buf2);
 335             if (gmx_strcasecmp(buf2, grpname) == 0)
 336             {
 337                 sol_line = nmol_line;
 338                 sscanf(buf, "%*s %d", &nsol_last);
 339             }
 340             /* Store this molecules section line */
 341             srenew(mol_line, nmol_line+1);
 342             mol_line[nmol_line] = strdup(buf);
 343             nmol_line++;
 344         }
 345     }
 346     ffclose(fpin);
 347
 348     if (sol_line == -1)
 349     {
 350         ffclose(fpout);
 351         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 352     }
 353     if (nsol_last < p_num+n_num)
 354     {
 355         ffclose(fpout);
 356         gmx_fatal(FARGS, "The last entry for moleculetype '%s' in the [ molecules ] section of file '%s' has less solvent molecules (%d) than were replaced (%d)", grpname, topinout, nsol_last, p_num+n_num);
 357     }
 358
 359     /* Print all the molecule entries */
 360     for (i = 0; i < nmol_line; i++)
 361     {
 362         if (i != sol_line)
 363         {
 364             fprintf(fpout, "%s", mol_line[i]);
 365         }
 366         else
 367         {
 368             printf("Replacing %d solute molecules in topology file (%s) "
 369                    " by %d %s and %d %s ions.\n",
 370                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 371             nsol_last -= p_num + n_num;
 372             if (nsol_last > 0)
 373             {
 374                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 375             }
 376             if (p_num > 0)
 377             {
 378                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 379             }
 380             if (n_num > 0)
 381             {
 382                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 383             }
 384         }
 385     }
 386     ffclose(fpout);
 387     /* use ffopen to generate backup of topinout */
 388     fpout = ffopen(topinout, "w");
 389     ffclose(fpout);
 390     rename(TEMP_FILENM, topinout);
 391 #undef TEMP_FILENM
 392 }
 393
 394 int gmx_genion(int argc, char *argv[])
 395 {
 396     const char        *desc[] = {
 397         "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
 398         "the position of the first atoms with the most favorable electrostatic",
 399         "potential or at random. The potential is calculated on all atoms, using",
 400         "normal GROMACS particle-based methods (in contrast to other methods",
 401         "based on solving the Poisson-Boltzmann equation).",
 402         "The potential is recalculated after every ion insertion.",
 403         "If specified in the run input file, a reaction field, shift function",
 404         "or user function can be used. For the user function a table file",
 405         "can be specified with the option [TT]-table[tt].",
 406         "The group of solvent molecules should be continuous and all molecules",
 407         "should have the same number of atoms.",
 408         "The user should add the ion molecules to the topology file or use",
 409         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 410         "The ion molecule type, residue and atom names in all force fields",
 411         "are the capitalized element names without sign. This molecule name",
 412         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 413         "[TT][molecules][tt] section of your topology updated accordingly,",
 414         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 415         "[PAR]Ions which can have multiple charge states get the multiplicity",
 416         "added, without sign, for the uncommon states only.[PAR]",
 417         "With the option [TT]-pot[tt] the potential can be written as B-factors",
 418         "in a [TT].pdb[tt] file (for visualisation using e.g. Rasmol).",
 419         "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
 420         "with the [TT]-scale[tt] option.[PAR]",
 421         "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
 422     };
 423     const char        *bugs[] = {
 424         "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
 425         "If you specify a salt concentration existing ions are not taken into account. In effect you therefore specify the amount of salt to be added."
 426     };
 427     static int         p_num   = 0, n_num = 0, p_q = 1, n_q = -1;
 428     static const char *p_name  = "NA", *n_name = "CL";
 429     static real        rmin    = 0.6, scale = 0.001, conc = 0;
 430     static int         seed    = 1993;
 431     static gmx_bool    bRandom = TRUE, bNeutral = FALSE;
 432     static t_pargs     pa[]    = {
 433         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 434         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
 435         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 436         { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 437         { "-nname", FALSE, etSTR,  {&n_name}, "Name of the negative ion"      },
 438         { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 439         { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 440         { "-random", FALSE, etBOOL, {&bRandom}, "Use random placement of ions instead of based on potential. The rmin option should still work" },
 441         { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" },
 442         { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
 443         { "-conc",  FALSE, etREAL, {&conc},
 444           "Specify salt concentration (mol/liter). This will add sufficient ions to reach up to the specified concentration as computed from the volume of the cell in the input [TT].tpr[tt] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
 445         { "-neutral", FALSE, etBOOL, {&bNeutral}, "This option will add enough ions to neutralize the system. These ions are added on top of those specified with [TT]-np[tt]/[TT]-nn[tt] or [TT]-conc[tt]. "}
 446     };
 447     gmx_mtop_t        *mtop;
 448     gmx_localtop_t    *top;
 449     t_inputrec         inputrec;
 450     t_commrec         *cr;
 451     t_mdatoms         *mdatoms;
 452     gmx_enerdata_t     enerd;
 453     t_graph           *graph;
 454     t_forcerec        *fr;
 455     rvec              *x, *v;
 456     real              *pot, vol, qtot;
 457     matrix             box;
 458     t_atoms            atoms;
 459     t_pbc              pbc;
 460     int               *repl;
 461     atom_id           *index;
 462     char              *grpname;
 463     gmx_bool          *bSet, bPDB;
 464     int                i, nw, nwa, nsa, nsalt, iqtot;
 465     FILE              *fplog;
 466     output_env_t       oenv;
 467     t_filenm           fnm[] = {
 468         { efTPX, NULL,  NULL,      ffREAD  },
 469         { efXVG, "-table", "table", ffOPTRD },
 470         { efNDX, NULL,  NULL,      ffOPTRD },
 471         { efSTO, "-o",  NULL,      ffWRITE },
 472         { efLOG, "-g",  "genion",  ffWRITE },
 473         { efPDB, "-pot", "pot",    ffOPTWR },
 474         { efTOP, "-p",  "topol",   ffOPTRW }
 475     };
 476 #define NFILE asize(fnm)
 477
 478     parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 479                       asize(desc), desc, asize(bugs), bugs, &oenv);
 480     bPDB = ftp2bSet(efPDB, NFILE, fnm);
 481     if (bRandom && bPDB)
 482     {
 483         fprintf(stderr, "Not computing potential with random option!\n");
 484         bPDB = FALSE;
 485     }
 486
 487     /* Check input for something sensible */
 488     if ((p_num < 0) || (n_num < 0))
 489     {
 490         gmx_fatal(FARGS, "Negative number of ions to add?");
 491     }
 492
 493     snew(mtop, 1);
 494     snew(top, 1);
 495     fplog = init_calcpot(ftp2fn(efLOG, NFILE, fnm), ftp2fn(efTPX, NFILE, fnm),
 496                          opt2fn("-table", NFILE, fnm), mtop, top, &inputrec, &cr,
 497                          &graph, &mdatoms, &fr, &enerd, &pot, box, &x, oenv);
 498
 499     atoms = gmx_mtop_global_atoms(mtop);
 500
 501     qtot = 0;
 502     for (i = 0; (i < atoms.nr); i++)
 503     {
 504         qtot += atoms.atom[i].q;
 505     }
 506     iqtot = gmx_nint(qtot);
 507
 508
 509     if (conc > 0)
 510     {
 511         /* Compute number of ions to be added */
 512         vol = det(box);
 513         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 514         p_num = abs(nsalt*n_q);
 515         n_num = abs(nsalt*p_q);
 516     }
 517     if (bNeutral)
 518     {
 519         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 520
 521         /* Check if the system is neutralizable
 522          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 523         int gcd = greatest_common_divisor(n_q, p_q);
 524         if ((qdelta % gcd) != 0)
 525         {
 526             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 527                     " -pq %d.\n", n_q, p_q);
 528         }
 529
 530         while (qdelta != 0)
 531         {
 532             while (qdelta < 0)
 533             {
 534                 p_num++;
 535                 qdelta += p_q;
 536             }
 537             while (qdelta > 0)
 538             {
 539                 n_num++;
 540                 qdelta += n_q;
 541             }
 542         }
 543     }
 544
 545     if ((p_num == 0) && (n_num == 0))
 546     {
 547         if (!bPDB)
 548         {
 549             fprintf(stderr, "No ions to add and no potential to calculate.\n");
 550             exit(0);
 551         }
 552         nw  = 0;
 553         nsa = 0; /* to keep gcc happy */
 554     }
 555     else
 556     {
 557         printf("Will try to add %d %s ions and %d %s ions.\n",
 558                p_num, p_name, n_num, n_name);
 559         printf("Select a continuous group of solvent molecules\n");
 560         get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 561         for (i = 1; i < nwa; i++)
 562         {
 563             if (index[i] != index[i-1]+1)
 564             {
 565                 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
 566                           "index[%d]=%d, index[%d]=%d",
 567                           grpname, i, index[i-1]+1, i+1, index[i]+1);
 568             }
 569         }
 570         nsa = 1;
 571         while ((nsa < nwa) &&
 572                (atoms.atom[index[nsa]].resind ==
 573                 atoms.atom[index[nsa-1]].resind))
 574         {
 575             nsa++;
 576         }
 577         if (nwa % nsa)
 578         {
 579             gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
 580                       nwa, nsa);
 581         }
 582         nw = nwa/nsa;
 583         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 584         if (p_num+n_num > nw)
 585         {
 586             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 587         }
 588     }
 589
 590     if (opt2bSet("-p", NFILE, fnm))
 591     {
 592         update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 593     }
 594
 595     snew(bSet, nw);
 596     snew(repl, nw);
 597
 598     snew(v, atoms.nr);
 599     snew(atoms.pdbinfo, atoms.nr);
 600
 601     set_pbc(&pbc, inputrec.ePBC, box);
 602
 603     /* Now loop over the ions that have to be placed */
 604     do
 605     {
 606         if (!bRandom)
 607         {
 608             calc_pot(fplog, cr, mtop, &inputrec, top, x, fr, &enerd, mdatoms, pot, box, graph);
 609             if (bPDB || debug)
 610             {
 611                 char buf[STRLEN];
 612
 613                 if (debug)
 614                 {
 615                     sprintf(buf, "%d_%s", p_num+n_num, ftp2fn(efPDB, NFILE, fnm));
 616                 }
 617                 else
 618                 {
 619                     strcpy(buf, ftp2fn(efPDB, NFILE, fnm));
 620                 }
 621                 for (i = 0; (i < atoms.nr); i++)
 622                 {
 623                     atoms.pdbinfo[i].bfac = pot[i]*scale;
 624                 }
 625                 write_sto_conf(buf, "Potential calculated by genion",
 626                                &atoms, x, v, inputrec.ePBC, box);
 627                 bPDB = FALSE;
 628             }
 629         }
 630         if ((p_num > 0) && (p_num >= n_num))
 631         {
 632             insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
 633                        1, p_q, p_name, mdatoms, rmin, bRandom, &seed);
 634             p_num--;
 635         }
 636         else if (n_num > 0)
 637         {
 638             insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
 639                        -1, n_q, n_name, mdatoms, rmin, bRandom, &seed);
 640             n_num--;
 641         }
 642     }
 643     while (p_num+n_num > 0);
 644     fprintf(stderr, "\n");
 645
 646     if (nw)
 647     {
 648         sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 649     }
 650
 651     sfree(atoms.pdbinfo);
 652     atoms.pdbinfo = NULL;
 653     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *mtop->name, &atoms, x, NULL,
 654                    inputrec.ePBC, box);
 655
 656     thanx(stderr);
 657
 658     gmx_log_close(fplog);
 659
 660     return 0;
 661 }