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