src/kernel/xutils.c

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