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