src/programs/mdrun/xutils.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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <assert.h>
  40
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "strdb.h"
  44 #include "string2.h"
  45 #include "xmdrun.h"
  46 #include "vec.h"
  47 #include "genalg.h"
  48 #include "random.h"
  49 #include "macros.h"
  50
  51 real mol_dipole(int k0, int k1, rvec x[], real q[])
  52 {
  53     int  k, m;
  54     rvec mu;
  55
  56     clear_rvec(mu);
  57     for (k = k0; (k < k1); k++)
  58     {
  59         for (m = 0; (m < DIM); m++)
  60         {
  61             mu[m] += q[k]*x[k][m];
  62         }
  63     }
  64     return norm(mu); /* Dipole moment of this molecule in e nm */
  65 }
  66
  67 real calc_mu_aver(t_commrec *cr, rvec x[], real q[], rvec mu,
  68                   t_block *mols, t_mdatoms *md, int gnx, atom_id grpindex[])
  69 {
  70     int     i, start, end;
  71     real    mu_ave;
  72
  73     start = md->start;
  74     end   = md->homenr + start;
  75
  76     /*
  77        clear_rvec(mu);
  78        for(i=start; (i<end); i++)
  79        for(m=0; (m<DIM); m++)
  80         mu[m] += q[i]*x[i][m];
  81        if (PAR(cr)) {
  82        gmx_sum(DIM,mu,cr);
  83        }
  84      */
  85     /* I guess we have to parallelise this one! */
  86
  87     if (gnx > 0)
  88     {
  89         mu_ave = 0.0;
  90         for (i = 0; (i < gnx); i++)
  91         {
  92             int gi = grpindex[i];
  93             mu_ave += mol_dipole(mols->index[gi], mols->index[gi+1], x, q);
  94         }
  95
  96         return(mu_ave/gnx);
  97     }
  98     else
  99     {
 100         return 0;
 101     }
 102 }
 103
 104 /* Lots of global variables! Yummy... */
 105 static t_ffscan ff;
 106
 107 void set_ffvars(t_ffscan *fff)
 108 {
 109     ff = *fff;
 110 }
 111
 112 real cost(tensor P, real MSF, real E)
 113 {
 114     return (ff.fac_msf*MSF+ff.fac_epot*sqr(E-ff.epot)+ff.fac_pres*
 115             (sqr(P[XX][XX]-ff.pres)+sqr(P[YY][YY]-ff.pres)+sqr(P[ZZ][ZZ]-ff.pres)));
 116
 117 }
 118
 119 static const char     *esenm[eseNR] = { "SIG", "EPS", "BHAMA", "BHAMB", "BHAMC", "CELlX", "CELLY", "CELLZ" };
 120 static int             nparm = 0, *param_val = NULL;
 121 static t_range        *range = NULL;
 122 static t_genalg       *ga    = NULL;
 123 static rvec            scale = { 1, 1, 1 };
 124
 125 static void init_range(t_range *r, int np, int atype, int ptype,
 126                        real rmin, real rmax)
 127 {
 128     if (rmin > rmax)
 129     {
 130         gmx_fatal(FARGS, "rmin (%f) > rmax (%f)", rmin, rmax);
 131     }
 132     if (np <= 0)
 133     {
 134         gmx_fatal(FARGS, "np (%d) should be > 0", np);
 135     }
 136     if ((rmax > rmin) && (np <= 1))
 137     {
 138         gmx_fatal(FARGS, "If rmax > rmin, np should be > 1");
 139     }
 140     if ((ptype < 0) || (ptype >= eseNR))
 141     {
 142         gmx_fatal(FARGS, "ptype (%d) should be < %d", ptype, eseNR);
 143     }
 144     r->np    = np;
 145     r->atype = atype;
 146     r->ptype = ptype;
 147     r->rmin  = rmin;
 148     r->rmax  = rmax;
 149     r->rval  = rmin;
 150     r->dr    = r->rmax - r->rmin;
 151 }
 152
 153 static t_range *read_range(const char *db, int *nrange)
 154 {
 155     int       nlines, nr, np, i;
 156     char    **lines;
 157     t_range  *range;
 158     int       atype, ptype;
 159     double    rmin, rmax;
 160
 161     nlines = get_file(db, &lines);
 162     snew(range, nlines);
 163
 164     nr = 0;
 165     for (i = 0; (i < nlines); i++)
 166     {
 167         strip_comment(lines[i]);
 168         if (sscanf(lines[i], "%d%d%d%lf%lf", &np, &atype, &ptype, &rmin, &rmax) == 5)
 169         {
 170             if (ff.bLogEps && (ptype == eseEPSILON) && (rmin <= 0))
 171             {
 172                 gmx_fatal(FARGS, "When using logarithmic epsilon increments the minimum"
 173                           "value must be > 0");
 174             }
 175             init_range(&range[nr], np, atype, ptype, rmin, rmax);
 176             nr++;
 177         }
 178     }
 179     fprintf(stderr, "found %d variables to iterate over\n", nr);
 180
 181     *nrange = nr;
 182
 183     for (nr = 0; (nr < nlines); nr++)
 184     {
 185         sfree(lines[nr]);
 186     }
 187     sfree(lines);
 188
 189     return range;
 190 }
 191
 192 static real value_range(t_range *r, int n)
 193 {
 194     real logrmin, logrmax;
 195
 196     if ((n < 0) || (n > r->np))
 197     {
 198         gmx_fatal(FARGS, "Value (%d) out of range for value_range (max %d)", n, r->np);
 199     }
 200
 201     if (r->np == 1)
 202     {
 203         r->rval = r->rmin;
 204     }
 205     else
 206     {
 207         if ((r->ptype == eseEPSILON) && ff.bLogEps)
 208         {
 209             logrmin = log(r->rmin);
 210             logrmax = log(r->rmax);
 211             r->rval = exp(logrmin + (n*(logrmax-logrmin))/(r->np-1));
 212         }
 213         else
 214         {
 215             r->rval = r->rmin+(n*(r->dr))/(r->np-1);
 216         }
 217     }
 218     return r->rval;
 219 }
 220
 221 real value_rand(t_range *r, int *seed)
 222 {
 223     real logrmin, logrmax;
 224     real mr;
 225
 226     if (r->np == 1)
 227     {
 228         r->rval = r->rmin;
 229     }
 230     else
 231     {
 232         mr = rando(seed);
 233         if ((r->ptype == eseEPSILON) && ff.bLogEps)
 234         {
 235             logrmin = log(r->rmin);
 236             logrmax = log(r->rmax);
 237             r->rval = exp(logrmin + mr*(logrmax-logrmin));
 238         }
 239         else
 240         {
 241             r->rval = r->rmin + mr*(r->rmax-r->rmin);
 242         }
 243     }
 244     if (debug)
 245     {
 246         fprintf(debug, "type: %s, value: %g\n", esenm[r->ptype], r->rval);
 247     }
 248     return r->rval;
 249 }
 250
 251 static void update_ff(t_forcerec *fr, int nparm, t_range range[], int param_val[])
 252 {
 253     static double *sigma = NULL, *eps = NULL, *c6 = NULL, *cn = NULL, *bhama = NULL, *bhamb = NULL, *bhamc = NULL;
 254     real           val, *nbfp;
 255     int            i, j, atnr;
 256
 257     atnr = fr->ntype;
 258     nbfp = fr->nbfp;
 259
 260     if (fr->bBHAM)
 261     {
 262         if (bhama == NULL)
 263         {
 264             assert(bhamb == NULL && bhamc == NULL);
 265             snew(bhama, atnr);
 266             snew(bhamb, atnr);
 267             snew(bhamc, atnr);
 268         }
 269     }
 270     else
 271     {
 272         if (sigma == NULL)
 273         {
 274             assert(eps == NULL && c6 == NULL && cn == NULL);
 275             snew(sigma, atnr);
 276             snew(eps, atnr);
 277             snew(c6, atnr);
 278             snew(cn, atnr);
 279         }
 280     }
 281     /* Get current values for everything */
 282     for (i = 0; (i < nparm); i++)
 283     {
 284         if (ga)
 285         {
 286             val = range[i].rval;
 287         }
 288         else
 289         {
 290             val = value_range(&range[i], param_val[i]);
 291         }
 292         if (debug)
 293         {
 294             fprintf(debug, "val = %g\n", val);
 295         }
 296         switch (range[i].ptype)
 297         {
 298             case eseSIGMA:
 299                 assert(!fr->bBHAM);
 300                 sigma[range[i].atype] = val;
 301                 break;
 302             case eseEPSILON:
 303                 assert(!fr->bBHAM);
 304                 eps[range[i].atype] = val;
 305                 break;
 306             case eseBHAMA:
 307                 assert(fr->bBHAM);
 308                 bhama[range[i].atype] = val;
 309                 break;
 310             case eseBHAMB:
 311                 assert(fr->bBHAM);
 312                 bhamb[range[i].atype] = val;
 313                 break;
 314             case eseBHAMC:
 315                 assert(fr->bBHAM);
 316                 bhamc[range[i].atype] = val;
 317                 break;
 318             case eseCELLX:
 319                 scale[XX] = val;
 320                 break;
 321             case eseCELLY:
 322                 scale[YY] = val;
 323                 break;
 324             case eseCELLZ:
 325                 scale[ZZ] = val;
 326                 break;
 327             default:
 328                 gmx_fatal(FARGS, "Unknown ptype");
 329         }
 330     }
 331     if (fr->bBHAM)
 332     {
 333         for (i = 0; (i < atnr); i++)
 334         {
 335             for (j = 0; (j <= i); j++)
 336             {
 337                 BHAMA(nbfp, atnr, i, j) = BHAMA(nbfp, atnr, j, i) = sqrt(bhama[i]*bhama[j]);
 338                 BHAMB(nbfp, atnr, i, j) = BHAMB(nbfp, atnr, j, i) = sqrt(bhamb[i]*bhamb[j]);
 339                 BHAMC(nbfp, atnr, i, j) = BHAMC(nbfp, atnr, j, i) = sqrt(bhamc[i]*bhamc[j]);
 340             }
 341         }
 342     }
 343     else
 344     {
 345         /* Now build a new matrix */
 346         for (i = 0; (i < atnr); i++)
 347         {
 348             c6[i] = 4*eps[i]*pow(sigma[i], 6.0);
 349             cn[i] = 4*eps[i]*pow(sigma[i], ff.npow);
 350         }
 351         for (i = 0; (i < atnr); i++)
 352         {
 353             for (j = 0; (j <= i); j++)
 354             {
 355                 C6(nbfp, atnr, i, j)  = C6(nbfp, atnr, j, i)  = sqrt(c6[i]*c6[j]);
 356                 C12(nbfp, atnr, i, j) = C12(nbfp, atnr, j, i) = sqrt(cn[i]*cn[j]);
 357             }
 358         }
 359     }
 360
 361     if (debug)
 362     {
 363         if (!fr->bBHAM)
 364         {
 365             for (i = 0; (i < atnr); i++)
 366             {
 367                 fprintf(debug, "atnr = %2d  sigma = %8.4f  eps = %8.4f\n", i, sigma[i], eps[i]);
 368             }
 369         }
 370         for (i = 0; (i < atnr); i++)
 371         {
 372             for (j = 0; (j < atnr); j++)
 373             {
 374                 if (fr->bBHAM)
 375                 {
 376                     fprintf(debug, "i: %2d  j: %2d  A:  %10.5e  B:  %10.5e  C:  %10.5e\n", i, j,
 377                             BHAMA(nbfp, atnr, i, j), BHAMB(nbfp, atnr, i, j), BHAMC(nbfp, atnr, i, j));
 378                 }
 379                 else
 380                 {
 381                     fprintf(debug, "i: %2d  j: %2d  c6:  %10.5e  cn:  %10.5e\n", i, j,
 382                             C6(nbfp, atnr, i, j), C12(nbfp, atnr, i, j));
 383                 }
 384             }
 385         }
 386     }
 387 }
 388
 389 static void scale_box(int natoms, rvec x[], matrix box)
 390 {
 391     int i, m;
 392
 393     if ((scale[XX] != 1.0) ||   (scale[YY] != 1.0) ||   (scale[ZZ] != 1.0))
 394     {
 395         if (debug)
 396         {
 397             fprintf(debug, "scale = %8.4f  %8.4f  %8.4f\n",
 398                     scale[XX], scale[YY], scale[ZZ]);
 399         }
 400         for (m = 0; (m < DIM); m++)
 401         {
 402             box[m][m] *= scale[m];
 403         }
 404         for (i = 0; (i < natoms); i++)
 405         {
 406             for (m = 0; (m < DIM); m++)
 407             {
 408                 x[i][m] *= scale[m];
 409             }
 410         }
 411     }
 412 }
 413
 414 gmx_bool update_forcefield(FILE *fplog,
 415                            int nfile, const t_filenm fnm[], t_forcerec *fr,
 416                            int natoms, rvec x[], matrix box)
 417 {
 418     static int ntry, ntried;
 419     int        i, j;
 420     gmx_bool   bDone;
 421
 422     /* First time around we have to read the parameters */
 423     if (nparm == 0)
 424     {
 425         range = read_range(ftp2fn(efDAT, nfile, fnm), &nparm);
 426         if (nparm == 0)
 427         {
 428             gmx_fatal(FARGS, "No correct parameter info in %s", ftp2fn(efDAT, nfile, fnm));
 429         }
 430         snew(param_val, nparm);
 431
 432         if (opt2bSet("-ga", nfile, fnm))
 433         {
 434             /* Genetic algorithm time */
 435             ga = init_ga(fplog, opt2fn("-ga", nfile, fnm), nparm, range);
 436         }
 437         else
 438         {
 439             /* Determine the grid size */
 440             ntry = 1;
 441             for (i = 0; (i < nparm); i++)
 442             {
 443                 ntry *= range[i].np;
 444             }
 445             ntried = 0;
 446
 447             fprintf(fplog, "Going to try %d different combinations of %d parameters\n",
 448                     ntry, nparm);
 449         }
 450     }
 451     if (ga)
 452     {
 453         update_ga(fplog, range, ga);
 454     }
 455     else
 456     {
 457         /* Increment the counter
 458          * Non-trivial, since this is nparm nested loops in principle
 459          */
 460         for (i = 0; (i < nparm); i++)
 461         {
 462             if (param_val[i] < (range[i].np-1))
 463             {
 464                 param_val[i]++;
 465                 for (j = 0; (j < i); j++)
 466                 {
 467                     param_val[j] = 0;
 468                 }
 469                 ntried++;
 470                 break;
 471             }
 472         }
 473         if (i == nparm)
 474         {
 475             fprintf(fplog, "Finished with %d out of %d iterations\n", ntried+1, ntry);
 476             return TRUE;
 477         }
 478     }
 479
 480     /* Now do the real updating */
 481     update_ff(fr, nparm, range, param_val);
 482
 483     /* Update box and coordinates if necessary */
 484     scale_box(natoms, x, box);
 485
 486     return FALSE;
 487 }
 488
 489 static void print_range(FILE *fp, tensor P, real MSF, real energy)
 490 {
 491     int  i;
 492
 493     fprintf(fp, "%8.3f  %8.3f  %8.3f  %8.3f",
 494             cost(P, MSF, energy), trace(P)/3, MSF, energy);
 495     for (i = 0; (i < nparm); i++)
 496     {
 497         fprintf(fp, " %s %10g", esenm[range[i].ptype], range[i].rval);
 498     }
 499     fprintf(fp, " FF\n");
 500     fflush(fp);
 501 }
 502
 503 static real msf(int n, rvec f1[], rvec f2[])
 504 {
 505     int  i, j;
 506     rvec ff2;
 507     real msf1 = 0;
 508
 509     for (i = 0; (i < n); )
 510     {
 511         clear_rvec(ff2);
 512         for (j = 0; ((j < ff.molsize) && (i < n)); j++, i++)
 513         {
 514             rvec_inc(ff2, f1[i]);
 515             if (f2)
 516             {
 517                 rvec_inc(ff2, f2[i]);
 518             }
 519         }
 520         msf1 += iprod(ff2, ff2);
 521     }
 522
 523     return msf1/n;
 524 }
 525
 526 static void print_grid(FILE *fp, real ener[], int natoms, rvec f[], rvec fshake[],
 527                        rvec x[], t_block *mols, real mass[], tensor pres)
 528 {
 529     static gmx_bool    bFirst = TRUE;
 530     static const char *desc[] = {
 531         "------------------------------------------------------------------------",
 532         "In the output from the forcefield scan we have the potential energy,",
 533         "then the root mean square force on the atoms, and finally the parameters",
 534         "in the order they appear in the input file.",
 535         "------------------------------------------------------------------------"
 536     };
 537     real               msf1;
 538     int                i;
 539
 540     if (bFirst)
 541     {
 542         for (i = 0; (i < asize(desc)); i++)
 543         {
 544             fprintf(fp, "%s\n", desc[i]);
 545         }
 546         fflush(fp);
 547         bFirst = FALSE;
 548     }
 549     if ((ff.tol == 0) || (fabs(ener[F_EPOT]/ff.nmol-ff.epot) < ff.tol))
 550     {
 551         msf1 = msf(natoms, f, fshake);
 552         if ((ff.f_max == 0) || (msf1 < sqr(ff.f_max)))
 553         {
 554             print_range(fp, pres, msf1, ener[F_EPOT]/ff.nmol);
 555         }
 556     }
 557 }
 558
 559 gmx_bool print_forcefield(FILE *fp, real ener[], int natoms, rvec f[], rvec fshake[],
 560                           rvec x[], t_block *mols, real mass[], tensor pres)
 561 {
 562     real msf1;
 563
 564     if (ga)
 565     {
 566         msf1 = msf(natoms, f, fshake);
 567         if (debug)
 568         {
 569             fprintf(fp, "Pressure: %12g, RMSF: %12g, Energy-Epot: %12g, cost: %12g\n",
 570                     ener[F_PRES], sqrt(msf1), ener[F_EPOT]/ff.nmol-ff.epot,
 571                     cost(pres, msf1, ener[F_EPOT]/ff.nmol));
 572         }
 573         if (print_ga(fp, ga, msf1, pres, scale, (ener[F_EPOT]/ff.nmol), range, ff.tol))
 574         {
 575             return TRUE;
 576         }
 577         fflush(fp);
 578     }
 579     else
 580     {
 581         print_grid(fp, ener, natoms, f, fshake, x, mols, mass, pres);
 582     }
 583     return FALSE;
 584 }