src/tools/correct.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32 static char *SRCID_correct_c = "$Id$";
  33 #include <math.h>
  34 #include "assert.h"
  35 #include "vec.h"
  36 #include "typedefs.h"
  37 #include "macros.h"
  38 #include "smalloc.h"
  39 #include "disco.h"
  40 #include "xtcio.h"
  41 #include "main.h"
  42 #include "random.h"
  43 #include "bondf.h"
  44 #include "names.h"
  45 #include "physics.h"
  46
  47 #ifdef debug_gmx
  48 #undef debug_gmx
  49 #endif
  50
  51 #define debug_gmx() do { FILE *fp=debug ? debug : (stdlog ? stdlog : stderr);\
  52 fprintf(fp,"PID=%d, %s  %d\n",gmx_cpu_id(),__FILE__,__LINE__); fflush(fp); } while (0)
  53
  54 #define NOT32
  55 #ifdef NOT32
  56 #define myrand(seed) rando(seed)
  57 #else
  58 static unsigned int IA=16807;
  59 static unsigned int IM=2147483647;
  60 static unsigned int IQ=127773;
  61 static unsigned int IR=2836;
  62 static unsigned int MASK=123459876;
  63 static float AM;
  64
  65 void init_rand()
  66 {
  67   AM = (1.0/IM);
  68 }
  69
  70 float myrand(int *seed)
  71 {
  72   int   k;
  73   float ans;
  74
  75   *seed ^= MASK;
  76   k      = (*seed)/IQ;
  77   *seed  = IA*(*seed-k*IQ)-IR*k;
  78   if (*seed < 0) *seed += IM;
  79   ans    = AM*(*seed);
  80   *seed ^= MASK;
  81
  82   return ans;
  83 }
  84 #endif
  85
  86 real mygauss(int *seed,real lb,real ub)
  87 {
  88   int  i;
  89   real myg=0.0;
  90
  91   for(i=0; (i<8); i++)
  92     myg += myrand(seed);
  93
  94   return lb+0.125*(ub-lb);
  95 }
  96
  97 void randomize_list(int n,int ip[],int *seed)
  98 {
  99   int i,iran,temp;
 100
 101   for(i=0; (i<n); i++) {
 102     iran     = ((int) (myrand(seed)*n)) % n;
 103     temp     = ip[i];
 104     ip[i]    = ip[iran];
 105     ip[iran] = temp;
 106   }
 107 }
 108
 109 void reset_list(int n,int ip[])
 110 {
 111   int i;
 112
 113   for(i=0; (i<n); i++)
 114     ip[i]    = i;
 115 }
 116
 117 int update_list(t_dist d[],int tag[],int natom,bool bViol[],int ip[])
 118 {
 119   int i,j,j0,j1,nchk;
 120
 121   nchk = 0;
 122   for(i=0; (i<natom); i++) {
 123     j0 = tag[i];
 124     j1 = tag[i+1];
 125     if (bViol[i]) {
 126       for(j=j0; (j<j1); j++) {
 127         if (d[j].ai != i)
 128           fatal_error(0,"Tags or distances inconsistent: "
 129                       "i=%d, j=%d, ai=%d, aj=%d, j0=%d, j1=%d",
 130                       i,j,d[j].ai,d[j].aj,j0,j1);
 131         ip[nchk] = j;
 132         nchk++;
 133       }
 134     }
 135     else {
 136       for(j=j0; (j<j1); j++) {
 137         if (d[j].ai != i)
 138           fatal_error(0,"Tags or distances inconsistent: "
 139                       "i=%d, j=%d, ai=%d, aj=%d, j0=%d, j1=%d",
 140                       i,j,d[j].ai,d[j].aj,j0,j1);
 141         if (bViol[d[j].aj]) {
 142           ip[nchk] = j;
 143           nchk++;
 144         }
 145       }
 146     }
 147   }
 148   return nchk;
 149 }
 150
 151 bool mirror_xx(int n1,int n2,int n3,int n4,rvec x[],bool bViol[],real xi)
 152 {
 153   /* Check the impropers, and fix if necessary */
 154   /* adapt the coordinates of atom n4 such that it is mirrorred
 155    * with respect to the plane defined by n1,n2,n3 (n1 origin)
 156    */
 157   rvec v1,v2,x4old,n,d14;
 158   real lambda,nsq,dn;
 159   int  m;
 160
 161   copy_rvec(x[n4],x4old);
 162
 163   /* v1 and v2 span the plane */
 164   rvec_sub(x[n2],x[n1],v1);
 165   rvec_sub(x[n3],x[n1],v2);
 166   rvec_sub(x[n1],x[n4],d14);
 167
 168   /* Take the normal vector to the plane */
 169   oprod(v1,v2,n);
 170   nsq = iprod(n,n);
 171
 172   /* Check for x[n1], x[n2], x[n3] on a line, if so come back next iteration */
 173   if (nsq == 0.0)
 174     return FALSE;
 175
 176   /* Find the point in which the vector parallel to n thru x[n4]
 177    * intersects with the plane:
 178    * by solving (x[n4]+lambda*n) . n = x[n1] . n
 179    * Compute lambda
 180    */
 181   lambda = iprod(d14,n)/nsq;
 182
 183   /* Check on the improper: If it is really far off than we should not
 184    * try to correct if, as this may only make it worse. The limits below
 185    * indicate a range that when mirrored will become a normal improper
 186    * again.
 187    */
 188   if ((xi < -10) && (xi > -60)) {
 189     /* The fourth atom is at the wrong side of the plane, fix it. */
 190
 191     /* Mirror all four atoms around the plane parallel to 123, going thru
 192      * the center of mass of the four atoms. It is easy to see that there
 193      * is no center-of-mass displacement this way.
 194      */
 195     for(m=0; (m<DIM); m++) {
 196       dn = 0.5*lambda*n[m];
 197       x[n1][m] -= dn;
 198       x[n2][m] -= dn;
 199       x[n3][m] -= dn;
 200       x[n4][m] += 3*dn;
 201     }
 202     if (debug) {
 203       fprintf(debug,"improper: %d %d %d %d, xi = %g deg.\n",
 204               n1+1,n2+1,n3+1,n4+1,xi);
 205       fprintf(debug,"lambda = %g, n = (%g,%g,%g)\n",lambda,n[XX],n[YY],n[ZZ]);
 206       fprintf(debug,"x4old = (%g,%g,%g) new = (%g,%g,%g)\n",
 207               x4old[XX],x4old[YY],x4old[ZZ],x[n4][XX],x[n4][YY],x[n4][ZZ]);
 208       bViol[n1] = TRUE;
 209       bViol[n2] = TRUE;
 210       bViol[n3] = TRUE;
 211       bViol[n4] = TRUE;
 212     }
 213     return TRUE;
 214   }
 215   return FALSE;
 216 }
 217
 218 void compute_dihs(FILE *log,int nq,t_quadruple q[],rvec x[],real phi[],
 219                   matrix box)
 220 {
 221   int    i;
 222   rvec   r_ij,r_kj,r_kl,m,n;
 223   real   cos_phi,sign;
 224
 225   for(i=0; (i<nq); i++)
 226     phi[i] = RAD2DEG*dih_angle(box,
 227                                x[q[i].ai],x[q[i].aj],x[q[i].ak],x[q[i].al],
 228                                r_ij,r_kj,r_kl,m,n,&cos_phi,&sign);
 229 }
 230
 231 int check_impropers(FILE *log,t_correct *c,int natom,rvec x[],matrix box)
 232
 233 {
 234   int  i,nmirror,nwrong;
 235   bool bFlip = FALSE;
 236
 237   nmirror = 0;
 238   if (c->bChiral) {
 239     compute_dihs(log,c->nimp,c->imp,x,c->idih,box);
 240     nwrong = 0;
 241     for(i=0; (i<c->nimp); i++)
 242       if (c->idih[i] < 0)
 243         nwrong++;
 244     if (nwrong > c->nimp/2) {
 245       if (debug)
 246         fprintf(debug,"Flipping entire molecule\n");
 247       for(i=0; (i<natom); i++)
 248         svmul(-1.0,x[i],x[i]);
 249     }
 250     for(i=0; (i<c->nimp); i++)
 251       if (mirror_xx(c->imp[i].ai,c->imp[i].aj,c->imp[i].ak,c->imp[i].al,
 252                     x,c->bViol,c->idih[i]))
 253         nmirror++;
 254     if (debug && (nmirror > 0))
 255       fprintf(debug,"mirrored %d improper dihedrals\n",nmirror);
 256   }
 257   if (bFlip)
 258     nmirror = c->nimp - nmirror;
 259
 260   return nmirror;
 261 }
 262
 263 int check_cis_trans(FILE *log,t_correct *c,rvec x[],matrix box)
 264 {
 265   int  i,nswap;
 266   rvec temp;
 267
 268   nswap = 0;
 269   if (c->bPep) {
 270     compute_dihs(log,c->npep,c->pepbond,x,c->omega,box);
 271     for(i=0; (i<c->npep); i++) {
 272       /* If one of the peptide bonds is cis, then swap H and Ca */
 273       if ((c->omega[i] < -90) || (c->omega[i] > 90)) {
 274         copy_rvec(x[c->pepbond[i].al],temp);
 275         copy_rvec(x[c->pepbond[i].am],x[c->pepbond[i].al]);
 276         copy_rvec(temp,x[c->pepbond[i].am]);
 277         /* These are now violated! */
 278         c->bViol[c->pepbond[i].al] = TRUE;
 279         c->bViol[c->pepbond[i].am] = TRUE;
 280         nswap++;
 281       }
 282     }
 283     if (debug && (nswap > 0))
 284       fprintf(debug,"swapped %d peptide bonds\n",nswap);
 285   }
 286   return nswap;
 287 }
 288
 289 bool do_1steep(int  ai,      int aj,       rvec dx,
 290                real len,     real lb,      real ub,     real wi,
 291                real *dev,    real *maxdev, real *ener, rvec fptr[],
 292                bool bViol[])
 293 {
 294   int  m;
 295   real ndev,wj,guess,ddx,delta;
 296
 297   /* Violated bounds ? */
 298   if ((len < lb) || (len > ub)) {
 299     /* Update counters */
 300     if (len < lb) {
 301       ndev  = lb-len;
 302       guess = lb;
 303     }
 304     else {
 305       ndev  = len-ub;
 306       guess = ub;
 307     }
 308
 309     *dev    += ndev;
 310     *ener   += (ndev*ndev);
 311     *maxdev  = max(*maxdev,ndev);
 312
 313     /* Set violation flag */
 314     bViol[ai] = TRUE;
 315     bViol[aj] = TRUE;
 316
 317     /* Correct non-bondeds only every N steps */
 318     wj  = 1.0-wi;
 319
 320     /* Calculate the force, that is the deviation from the nearest
 321      * bound. This is the steepest descent minimizaton part.
 322      */
 323     delta = (len-guess)/len;
 324     for(m=0; (m<DIM); m++) {
 325       ddx = delta*dx[m];
 326       fptr[ai][m] += wi*ddx;
 327       fptr[aj][m] -= wj*ddx;
 328     }
 329     return TRUE;
 330   }
 331   else
 332     return FALSE;
 333 }
 334
 335 bool do_1shake(int  cons_type,
 336                int  ai,      int aj,       rvec dx,
 337                real len,     real lb,      real ub,     real wi,
 338                real *dev,    real *maxdev, real *ener,  rvec xptr[],
 339                bool bViol[], int  *seed,   bool bLowerOnly)
 340 {
 341   int  m;
 342   real ndev,delta,ddx,wj,guess;
 343
 344   /* Violated bounds ? */
 345   if ((len < lb) || ((cons_type != edcNONBOND) && (len > ub))) {
 346     /* Update counters */
 347     if (cons_type == edcNONBOND)
 348       ndev = lb-len;
 349     else {
 350       if (len < lb)
 351         ndev = lb-len;
 352       else
 353         ndev = len-ub;
 354     }
 355     *dev    += ndev;
 356     *ener   += (ndev*ndev);
 357     *maxdev  = max(*maxdev,ndev);
 358
 359     /* Set violation flag */
 360     bViol[ai] = TRUE;
 361     bViol[aj] = TRUE;
 362
 363     wj  = 1.0-wi;
 364     if (cons_type == edcBOND)
 365       guess = mygauss(seed,lb,ub);
 366     else
 367       guess = lb + myrand(seed)*(ub-lb);
 368
 369     delta = (len-guess)/len;
 370     for(m=0; (m<DIM); m++) {
 371       ddx = delta*dx[m];
 372       xptr[ai][m] += wi*ddx;
 373       xptr[aj][m] -= wj*ddx;
 374     }
 375     return TRUE;
 376   }
 377   else
 378     return FALSE;
 379 }
 380
 381 int shake_coords(FILE *log,bool bVerbose,
 382                  int nstruct,int natom,rvec xref[],
 383                  rvec x[],int *seed,matrix box,t_correct *c,int *niter)
 384 {
 385   static bool bFirst=TRUE;
 386   static rvec *force[2],*xtry[2],*xwr;
 387   int    nit,nchk,low,oldlow;
 388   int    i,j,k,l,m,ai,aj,status,nprint;
 389   int    nviol[edcNR],nvtot,nswap,nmirror,cons_type;
 390   bool   bShort,bConverged,bNB,bPrint,bLowDev,bExplicit,bMyNB,bLowerOnly;
 391   bool   bVNow;
 392   bool   *bViol;
 393   real   lb,ub,lbfac,wi,len,dev,maxdev,len2;
 394   real   step,maxforce=0,rstep,step0;
 395   real   ener[2];
 396   int    cur=0;
 397 #define next 1-cur
 398   t_dist *d;
 399   int    *ip;
 400   rvec   dx,xcm,*xptr,*fptr;
 401   char   buf[32];
 402
 403   /* Initiate local variables */
 404   bViol      = c->bViol;
 405   d          = c->d;
 406   ip         = c->ip;
 407   nchk       = c->ndist;
 408   low        = c->ndist;
 409   oldlow     = c->ndist;
 410   bExplicit  = c->bExplicit;
 411   bConverged = FALSE;
 412   bLowDev    = FALSE;
 413   bLowerOnly = c->bLowerOnly;
 414   step0      = 0.1;
 415   step       = step0;
 416   nprint     = 0;
 417   ener[cur]  = 1e30;
 418   status     = 0;
 419
 420   debug_gmx();
 421   if (bFirst) {
 422     snew(force[0],natom);
 423     snew(force[1],natom);
 424     snew(xtry[0],natom);
 425     snew(xtry[1],natom);
 426     snew(xwr,natom);
 427 #ifndef NOT32
 428     init_rand();
 429 #endif
 430     bFirst = FALSE;
 431   }
 432   if (c->bDump) {
 433     sprintf(buf,"test%d.xtc",nstruct);
 434     status     = open_xtc(buf,"w");
 435     center_in_box(natom,x,box,xwr);
 436     write_xtc(status,natom,0,0,box,xwr,1000.0);
 437   }
 438   for(i=0; (i<natom); i++) {
 439     copy_rvec(x[i],xtry[cur][i]);
 440     copy_rvec(x[i],xtry[next][i]);
 441   }
 442   xptr = xtry[next];
 443   debug_gmx();
 444
 445   if (c->bRanlistFirst)
 446     randomize_list(nchk,ip,seed);
 447
 448   for(nit=0; ((nit < c->maxnit) && !bConverged); nit++) {
 449     /* Reset violation counters */
 450     for(i=0; (i<natom); i++)
 451       bViol[i] = FALSE;
 452
 453     /* Set progress indicators */
 454     nviol[edcNONBOND] = 0;
 455     nviol[edcBOND]    = 0;
 456     nviol[edcDISRE]   = 0;
 457     dev        = 0.0;
 458     maxdev     = 0.0;
 459     ener[next] = 0.0;
 460
 461     /* Nonbonded growing factor */
 462     lbfac      = (c->ngrow == 0) ? 1.0 : (min(1.0,(real)nit/(real)c->ngrow));
 463
 464     /* Set step dependent flags */
 465     bNB        = (c->nbcheck != 0) && (nit % c->nbcheck)  == 0;
 466     bPrint     = (c->nstprint!= 0) && (nit % c->nstprint) == 0;
 467     bExplicit  = bLowDev && c->bExplicit;
 468
 469     if (c->nstranlist && ((nit % c->nstranlist) == 0))
 470       randomize_list(nchk,ip,seed);
 471
 472     xptr = xtry[next];
 473     if (bExplicit) {
 474       fptr = force[next];
 475       for(i=0; (i<natom); i++)
 476         clear_rvec(fptr[i]);
 477     }
 478     else
 479       fptr = NULL;
 480
 481     /* Loop over distance list */
 482     for(l=0; (l<nchk); l++) {
 483       /* Unpack lists */
 484       k         = ip[l];
 485       ai        = d[k].ai;
 486       aj        = d[k].aj;
 487       lb        = d[k].lb;
 488       ub        = d[k].ub;
 489       wi        = d[k].wi;
 490       cons_type = d[k].cons_type;
 491
 492       if (cons_type == edcNONBOND) {
 493         lb *= lbfac;
 494         bMyNB = bNB;
 495       }
 496       else
 497         bMyNB = FALSE;
 498
 499       if (bMyNB || (cons_type != edcNONBOND)) {
 500         /* Compute distance */
 501         len2 = 0.0;
 502         for(m=0; (m<DIM); m++) {
 503           dx[m] = xptr[aj][m] - xptr[ai][m];
 504           len2 += dx[m]*dx[m];
 505         }
 506         len = sqrt(len2);
 507
 508         if (bExplicit)
 509           bVNow = do_1steep(ai,aj,dx,len,lb,ub,wi,&dev,&maxdev,&ener[next],
 510                             fptr,bViol);
 511         else {
 512           bVNow = do_1shake(cons_type,ai,aj,
 513                             dx,len,lb,ub,wi,&dev,&maxdev,&ener[next],
 514                             xptr,bViol,seed,bLowerOnly);
 515 #ifdef HEAVY
 516           if (bVNow && c->bDump && debug) {
 517             center_in_box(natom,xptr,box,xwr);
 518             write_xtc(status,natom,nit,(real)l,box,xwr,1000.0);
 519           }
 520 #endif
 521         }
 522         if (bVNow)
 523           nviol[cons_type]++;
 524       }
 525     }
 526     if (bExplicit) {
 527       if (ener[next] < ener[cur]) {
 528         cur = next;
 529         maxforce = 0;
 530         for(i=0; (i<natom); i++)
 531           for(m=0; (m<DIM); m++)
 532             maxforce = max(maxforce,fabs(fptr[i][m]));
 533         /* Increase step again */
 534         step  = min(step0,2*step);
 535         rstep = step/maxforce;
 536         for(i=0; (i<natom); i++)
 537           for(m=0; (m<DIM); m++)
 538             xtry[next][i][m] = xtry[cur][i][m] + fptr[i][m]*rstep;
 539         xptr = xtry[next];
 540       }
 541       else {
 542         /* Decrease step size and try again */
 543         step  = step/2;
 544         rstep = step/maxforce;
 545         for(i=0; (i<natom); i++)
 546           for(m=0; (m<DIM); m++)
 547             xtry[next][i][m] = xtry[cur][i][m] + fptr[i][m]*rstep;
 548       }
 549     }
 550     if (c->bDump && bPrint) {
 551       center_in_box(natom,xptr,box,xwr);
 552       write_xtc(status,natom,nit,(real)nit+1,box,xwr,1000.0);
 553     }
 554     /* Check Impropers */
 555     nmirror = check_impropers(bVerbose ? log : NULL,c,natom,xptr,box);
 556
 557     /* Check peptide bonds */
 558     nswap   = check_cis_trans(bVerbose ? log : NULL,c,xptr,box);
 559
 560     nvtot   = (nviol[edcNONBOND] + nviol[edcBOND] + nviol[edcDISRE] +
 561                nmirror + nswap);
 562     if (bPrint) {
 563       if ((nprint % 20) == 0)
 564         fprintf(stderr,"\n%5s%7s%7s%7s%7s%7s%7s%7s %10s %10s\n",
 565                 "nit","nchk","nimp","npep",
 566                 edc_names[edcNONBOND],edc_names[edcBOND],edc_names[edcDISRE],
 567                 "nvtot","maxdev",
 568                 bExplicit ? "ener A^2" : "<dev>");
 569       fprintf(stderr,"%5d%7d%7d%7d%7d%7d%7d%7d %10.4e %10.4e\n",
 570               nit,nchk,nmirror,nswap,
 571               nviol[edcNONBOND],nviol[edcBOND],nviol[edcDISRE],nvtot,maxdev,
 572               bExplicit ? ener[cur] : dev/natom);
 573       nprint++;
 574     }
 575     bLowDev    = (bLowDev || (dev/natom < c->lodev));
 576
 577     /* Update neighbour list */
 578 #define UPDL
 579 #ifdef UPDL
 580     nchk = update_list(c->d,c->tag,natom,bViol,ip);
 581
 582     bConverged = (nchk == 0);
 583 #else
 584     reset_list(nchk,ip);
 585
 586     bConverged = (nvtot == 0);
 587 #endif
 588   }
 589   if (c->bDump)
 590     close_xtc(status);
 591
 592   /* Copy final structure back to input array */
 593   for(i=0; (i<natom); i++)
 594     copy_rvec(xptr[i],x[i]);
 595
 596   fprintf(log,"struct: %5d, box: %8.3f  %8.3f  %8.3f\n",
 597           nstruct,box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
 598   fprintf(log,"Converged: %3s nit: %5d  <dev>: %8.4f  maxdev: %8.4f\n",
 599           yesno_names[bConverged],nit,dev/natom,maxdev);
 600
 601   *niter = nit;
 602
 603   return nvtot;
 604 }
 605
 606 int quick_check(FILE *log,int natom,rvec x[],matrix box,t_correct *c)
 607 {
 608   int  i,j,k,m,ai,aj,nviol;
 609   real lb,ub,len,len2,dx;
 610
 611   if (log)
 612     fprintf(log,"Double checking final structure\n");
 613   nviol = 0;
 614   for(k=0; (k<c->ndist); k++) {
 615     /* Unpack lists */
 616     ai    = c->d[k].ai;
 617     aj    = c->d[k].aj;
 618     lb    = c->d[k].lb;
 619     ub    = c->d[k].ub;
 620
 621     /* Compute distance */
 622     len2 = 0.0;
 623     for(m=0; (m<DIM); m++) {
 624       dx    = x[aj][m] - x[ai][m];
 625       len2 += dx*dx;
 626     }
 627     len = sqrt(len2);
 628     if ((len < lb) || ((ub > 0) && (len > ub))) {
 629       if (debug)
 630         fprintf(debug,"VIOL: ai=%4d aj=%4d lb=%8.4f len=%8.4f ub=%8.4f\n",
 631                 ai,aj,lb,len,ub);
 632       nviol++;
 633     }
 634   }
 635   /* Check Impropers */
 636   /* nmirror = check_impropers(log,c,natom,x,box);*/
 637
 638   /* Check peptide bonds */
 639   /* nswap   = check_cis_trans(log,c,x,box);*/
 640
 641   return nviol;
 642 }