src/tools/g_relax.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_relax_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <stdlib.h>
  33 #include "sysstuff.h"
  34 #include "string.h"
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "vec.h"
  39 #include "xvgr.h"
  40 #include "pbc.h"
  41 #include "copyrite.h"
  42 #include "futil.h"
  43 #include "statutil.h"
  44 #include "rdgroup.h"
  45 #include "gstat.h"
  46 #include "tpxio.h"
  47
  48 typedef struct {
  49   real re,im;
  50 } complex;
  51
  52 typedef struct {
  53   int ai,aj;
  54 } t_pair;
  55
  56 typedef struct {
  57   real    rij_3;
  58   real    rij_6;
  59   bool    bNOE;
  60   real    tauc,dtauc,S2,dS2;
  61   complex y2;
  62   complex Ylm[5];
  63 } t_sij;
  64
  65 void read_shifts(FILE *fp,int nx,real shiftx[],int ny,real shifty[])
  66 {
  67   int    i;
  68   double d;
  69
  70   for(i=0; (i<nx); i++) {
  71     fscanf(fp,"%lf",&d);
  72     shiftx[i]=d;
  73     shifty[i]=d;
  74   }
  75   /*for(i=0; (i<ny); i++) {
  76     fscanf(fp,"%lf",&d);
  77     shifty[i]=d;
  78     }
  79     */
  80 }
  81
  82 complex c_sqr(complex c)
  83 {
  84   complex cs;
  85
  86   cs.re = c.re*c.re-c.im*c.im;
  87   cs.im = 2.0*c.re*c.im;
  88
  89   return cs;
  90 }
  91
  92 complex c_add(complex c,complex d)
  93 {
  94   complex cs;
  95
  96   cs.re = c.re+d.re;
  97   cs.im = c.im+d.im;
  98
  99   return cs;
 100 }
 101
 102 complex calc_ylm(int m,rvec rij,real r2,real r_3,real r_6)
 103 {
 104   static  bool bFirst=TRUE;
 105   static  real y0,y1,y2;
 106   real    x,y,z,xsq,ysq,rxy,r1,cphi,sphi,cth,sth,fac;
 107   complex cs;
 108
 109   if (bFirst) {
 110     y0 = sqrt(5/(4*M_PI));
 111     y1 = sqrt(45/(24*M_PI));
 112     y2 = sqrt(45/(96*M_PI));
 113     bFirst = FALSE;
 114   }
 115   r1  = sqrt(r2);
 116   x   = rij[XX];
 117   y   = rij[YY];
 118   z   = rij[ZZ];
 119   xsq = x*x;
 120   ysq = y*y;
 121   rxy = sqrt(xsq+ysq);
 122   cphi= x/rxy;
 123   sphi= y/rxy;
 124   cth = z/r1;
 125   if (cphi != 0.0)
 126     sth = x/(r1*cphi);
 127   else
 128     sth = y/(r1*sphi);
 129
 130   /* Now calculate the spherical harmonics */
 131   switch(m) {
 132   case -2:
 133   case  2:
 134     fac    = y2*sth*sth;
 135     cs.re  = fac*(cphi*cphi-sphi*sphi);
 136     /* Use the index as a prefactor... check your formulas */
 137     cs.im  = m*fac*cphi*sphi;
 138     break;
 139   case -1:
 140   case  1:
 141     fac    = y1*cth*sth;
 142     cs.re  = -m*fac*cphi;
 143     cs.im  = -fac*sphi;
 144     break;
 145   case 0:
 146     cs.re = y0*(3*cth*cth-1);
 147     cs.im = 0;
 148     break;
 149   }
 150   cs.re *= r_3;
 151   cs.im *= r_3;
 152
 153   return cs;
 154 }
 155
 156 void myfunc(real x,real a[],real *y,real dyda[],int na)
 157 {
 158   /* Fit to function
 159    *
 160    * y = a1 + (1-a1) exp(-a2 x)
 161    *
 162    * where in real life a1 is S^2 and a2 = 1/tauc
 163    *
 164    */
 165
 166   real eee,S2,tau1;
 167
 168   S2      = a[1];
 169   tau1    = 1.0/a[2];
 170   eee     = exp(-x*tau1);
 171   *y      = S2 + (1-S2)*eee;
 172   dyda[1] = 1 - eee;
 173   dyda[2] = x*tau1*tau1*(1-a[1])*eee;
 174 }
 175
 176 void fit_one(bool bVerbose,
 177              int nframes,real x[],real y[],real dy[],real ftol,
 178              real *S2,real *dS2,real *tauc,real *dtauc)
 179 {
 180   void mrqmin(real x[],real y[],real sig[],int ndata,real a[],
 181               int ma,int lista[],int mfit,real **covar,real **alpha,
 182               real *chisq,
 183               void (*funcs)(real x,real a[],real *y,real dyda[],int na),
 184               real *alamda);
 185
 186   real *a,**covar,**alpha;
 187   real chisq,ochisq,alamda;
 188   bool bCont;
 189   int  i,j,ma,mfit,*lista;
 190
 191   ma=mfit=2;
 192   snew(a,ma+1);
 193   snew(covar,ma+1);
 194   snew(alpha,ma+1);
 195   snew(lista,ma+1);
 196   for(i=0; (i<ma+1); i++) {
 197     lista[i] = i;
 198     snew(covar[i],ma+1);
 199     snew(alpha[i],ma+1);
 200   }
 201
 202   a[1]   = 0.99;  /* S^2              */
 203   a[2]   = 0.1;   /* tauc             */
 204   alamda = -1;    /* Starting value   */
 205   chisq  = 1e12;
 206   j      = 0;
 207   do {
 208     ochisq = chisq;
 209     mrqmin(x-1,y-1,dy-1,nframes,a,ma,lista,mfit,covar,alpha,
 210            &chisq,myfunc,&alamda);
 211     if (bVerbose)
 212       fprintf(stderr,"\rFitting %d chisq=%g, alamda=%g, tau=%g, S^2=%g\t\t\n",
 213               j,chisq,alamda,1.0/a[2],a[1]);
 214     j++;
 215     bCont = (((ochisq - chisq) > ftol*chisq) ||
 216              ((ochisq == chisq)));
 217   } while (bCont && (alamda != 0.0) && (j < 50));
 218   if (bVerbose)
 219     fprintf(stderr,"\n");
 220
 221   /* Now get the covariance matrix out */
 222   alamda = 0;
 223   mrqmin(x-1,y-1,dy-1,nframes,a,ma,lista,mfit,covar,alpha,
 224          &chisq,myfunc,&alamda);
 225
 226   *S2    = a[1];
 227   *dS2   = sqrt(covar[1][1]);
 228   *tauc  = a[2];
 229   *dtauc = sqrt(covar[2][2]);
 230
 231   for(i=0; (i<ma+1); i++) {
 232     sfree(covar[i]);
 233     sfree(alpha[i]);
 234   }
 235   sfree(a);
 236   sfree(covar);
 237   sfree(alpha);
 238   sfree(lista);
 239 }
 240
 241 void calc_tauc(bool bVerbose,int npair,t_pair pair[],real dt,
 242                int nframes,t_sij spec[],real **corr)
 243 {
 244   FILE *fp;
 245   char buf[32];
 246   int  i,j,k,n;
 247   real S2,S22,tauc,fac;
 248   real *x,*dy;
 249   real ftol = 1e-3;
 250
 251   snew(x,nframes);
 252   snew(dy,nframes);
 253   for(i=0; (i<nframes); i++)
 254     x[i]  = i*dt;
 255
 256   fprintf(stderr,"Fitting correlation function to Lipari&Szabo function\n");
 257   fac=1.0/((real)nframes);
 258   for(i=0; (i<npair); i++) {
 259     if (spec[i].bNOE) {
 260       /* Use Levenberg-Marquardt method to fit */
 261       for(j=0; (j<nframes); j++)
 262         dy[j] = fac;
 263       fit_one(bVerbose,
 264               nframes,x,corr[i],dy,ftol,
 265               &(spec[i].S2),&(spec[i].dS2),
 266               &(spec[i].tauc),&(spec[i].dtauc));
 267       if (bVerbose) {
 268         sprintf(buf,"test%d.xvg",i);
 269         fp = ffopen(buf,"w");
 270         for(j=0; (j<nframes); j++) {
 271           fprintf(fp,"%10g  %10g  %10g\n",j*dt,corr[i][j],
 272                   spec[i].S2 + (1-spec[i].S2)*exp(-j*dt/spec[i].tauc));
 273         }
 274         fclose(fp);
 275       }
 276     }
 277   }
 278 }
 279
 280 void calc_aver(FILE *fp,int nframes,int npair,t_pair pair[],t_sij *spec,
 281                real maxdist)
 282 {
 283   int     i,j,m;
 284   real    nf_1,fac,md_6;
 285   complex c1,c2,dc2;
 286
 287   md_6 = pow(maxdist,-6.0);
 288   fac  = 4*M_PI/5;
 289   nf_1 = 1.0/nframes;
 290   for(i=0; (i<npair); i++) {
 291     c2.re = 0;
 292     c2.im = 0;
 293     fprintf(fp,"%5d  %5d",pair[i].ai,pair[i].aj);
 294     for(m=0; (m<5); m++) {
 295       c1.re  = spec[i].Ylm[m].re*nf_1;
 296       c1.im  = spec[i].Ylm[m].im*nf_1;
 297       dc2    = c_sqr(c1);
 298       c2     = c_add(dc2,c2);
 299
 300       if (c1.im > 0)
 301         fprintf(fp,"  %8.3f+i%8.3f",c1.re,c1.im);
 302       else
 303         fprintf(fp,"  %8.3f-i%8.3f",c1.re,-c1.im);
 304     }
 305     fprintf(fp,"\n");
 306     spec[i].rij_3 *= nf_1;
 307     spec[i].rij_6 *= nf_1;
 308     spec[i].y2.re  = fac*c2.re;
 309     spec[i].y2.im  = fac*c2.im;
 310     spec[i].bNOE   = (spec[i].rij_6 > md_6);
 311   }
 312 }
 313
 314 void plot_spectrum(char *noefn,int npair,t_pair pair[],t_sij *spec,real taum)
 315 {
 316   FILE    *fp,*out;
 317   int     i,j,m;
 318   t_rgb   rlo = { 1,0,0 },rhi = {1,1,1};
 319   real    Sijmax,Sijmin,pow6,pow3,pp3,pp6,ppy,tauc;
 320   real    *Sij;
 321   complex sij;
 322
 323   snew(Sij,npair);
 324   Sijmax = -1000.0;
 325   Sijmin =  1000.0;
 326   fp=xvgropen(noefn,"Cross Relaxation","Pair Index","\\8s\\4\\sij\\N");
 327   for(i=0; (i<npair); i++) {
 328     tauc   = spec[i].tauc;
 329     sij.re = -0.4*((taum-tauc)*spec[i].y2.re + tauc*spec[i].rij_6);
 330     sij.im = -0.4* (taum-tauc)*spec[i].y2.im;
 331     Sij[i]=sij.re;
 332     Sijmax=max(Sijmax,sij.re);
 333     Sijmin=min(Sijmin,sij.re);
 334     fprintf(fp,"%5d  %10g\n",i,sij.re);
 335   }
 336   fclose(fp);
 337   fprintf(stderr,"Sijmin: %g, Sijmax: %g\n",Sijmin,Sijmax);
 338   out=ffopen("spec.out","w");
 339   pow6 = -1.0/6.0;
 340   pow3 = -1.0/3.0;
 341   fprintf(out,"%5s  %5s  %8s  %8s  %8s  %8s  %8s  %8s  %8s  %8s  %8s\n",
 342           "at i","at j","S^2","Sig S^2","tauc","Sig tauc",
 343           "<rij6>","<rij3>","<ylm>","rij3-6","ylm-rij6");
 344   for(i=0; (i<npair); i++) {
 345     if (spec[i].bNOE) {
 346       pp6 = pow(spec[i].rij_6,pow6);
 347       pp3 = pow(spec[i].rij_3,pow3);
 348       if (spec[i].y2.re < 0)
 349         ppy = -pow(-spec[i].y2.re,pow6);
 350       else
 351         ppy = pow(spec[i].y2.re,pow6);
 352       fprintf(out,"%5d  %5d  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f\n",
 353               pair[i].ai,pair[i].aj,
 354               spec[i].S2,spec[i].dS2,spec[i].tauc,spec[i].dtauc,
 355               pp6,pp3,ppy,pp3-pp6,ppy-pp6);
 356     }
 357   }
 358   fclose(out);
 359
 360   sfree(Sij);
 361
 362   xvgr_file(noefn,NULL);
 363 }
 364
 365 void spectrum(bool bVerbose,
 366               char *trj,char *shifts,bool bAbInitio,
 367               char *corrfn,char *noefn,
 368               int maxframes,bool bFour,bool bFit,int nrestart,
 369               int npair,t_pair pair[],int nat,real chem_shifts[],
 370               real taum,real maxdist,
 371               real w_rls[],rvec xp[],t_idef *idef)
 372 {
 373   FILE   *fp;
 374   int    i,j,m,ii,jj,natoms,status,nframes;
 375   rvec   *x,dx;
 376   matrix box;
 377   real   t0,t1,t,dt;
 378   real   r2,r6,r_3,r_6,tauc;
 379   rvec   **corr;
 380   real   **Corr;
 381   t_sij  *spec;
 382
 383   snew(spec,npair);
 384
 385   fprintf(stderr,"There is no kill like overkill! Going to malloc %d bytes\n",
 386           npair*maxframes*sizeof(corr[0][0]));
 387   snew(corr,npair);
 388   for(i=0; (i<npair); i++)
 389     snew(corr[i],maxframes);
 390   nframes = 0;
 391   natoms  = read_first_x(&status,trj,&t0,&x,box);
 392   if (natoms > nat)
 393     fatal_error(0,"Not enough atoms in trajectory");
 394   do {
 395     if (nframes >= maxframes) {
 396       fprintf(stderr,"\nThere are more than the %d frames you told me!",
 397               maxframes);
 398       break;
 399     }
 400     t1 = t;
 401     if (bVerbose)
 402       fprintf(stderr,"\rframe: %d",nframes);
 403     rm_pbc(idef,natoms,box,x,x);
 404     if (bFit)
 405       do_fit(natoms,w_rls,xp,x);
 406
 407     for(i=0; (i<npair); i++) {
 408       ii = pair[i].ai;
 409       jj = pair[i].aj;
 410       rvec_sub(x[ii],x[jj],dx);
 411       copy_rvec(dx,corr[i][nframes]);
 412
 413       r2  = iprod(dx,dx);
 414       r6  = r2*r2*r2;
 415       r_3 = invsqrt(r6);
 416       r_6 = r_3*r_3;
 417       spec[i].rij_3 += r_3;
 418       spec[i].rij_6 += r_6;
 419       for(m=0; (m<5); m++) {
 420         spec[i].Ylm[m] = c_add(spec[i].Ylm[m],
 421                                calc_ylm(m-2,dx,r2,r_3,r_6));
 422       }
 423     }
 424     nframes++;
 425   } while (read_next_x(status,&t,natoms,x,box));
 426   close_trj(status);
 427   if (bVerbose)
 428     fprintf(stderr,"\n");
 429
 430   fp=ffopen("ylm.out","w");
 431   calc_aver(fp,nframes,npair,pair,spec,maxdist);
 432   fclose(fp);
 433
 434   /* Select out the pairs that have to be correlated */
 435   snew(Corr,npair);
 436   for(i=j=0; (i<npair); i++) {
 437     if (spec[i].bNOE) {
 438       Corr[j] = &(corr[i][0][0]);
 439       j++;
 440     }
 441   }
 442   fprintf(stderr,"There are %d NOEs in your simulation\n",j);
 443   if (nframes > 1)
 444     dt = (t1-t0)/(nframes-1);
 445   else
 446     dt = 1;
 447   do_autocorr(corrfn,"Correlation Function for Interproton Vectors",
 448               nframes,j,Corr,dt,eacP2,nrestart,FALSE,FALSE,bFour,TRUE);
 449
 450   calc_tauc(bVerbose,npair,pair,dt,nframes/2,spec,(real **)corr);
 451
 452   plot_spectrum(noefn,npair,pair,spec,taum);
 453 }
 454
 455 int main(int argc,char *argv[])
 456 {
 457   static char *desc[] = {
 458     "g_noe calculates a NOE spectrum"
 459   };
 460
 461   int        status;
 462   t_topology *top;
 463   int        i,j,k,natoms,nprot,*prot_ind;
 464   int        ifit;
 465   char       *gn_fit;
 466   atom_id    *ind_fit,*all_at;
 467   real       *w_rls;
 468   rvec       *xp;
 469   t_pair     *pair;
 470   matrix     box;
 471   int        step,nre;
 472   real       t,lambda;
 473   real       *shifts=NULL;
 474   t_filenm   fnm[] = {
 475     { efTRX, "-f", NULL,     ffREAD },
 476     { efTPX, "-s", NULL,     ffREAD },
 477     { efNDX, NULL, NULL,     ffREAD },
 478     { efDAT, "-d", "shifts", ffREAD },
 479     { efOUT, "-o","spec",    ffWRITE },
 480     { efXVG, "-corr", "rij-corr", ffWRITE },
 481     { efXVG, "-noe", "noesy", ffWRITE }
 482   };
 483 #define NFILE asize(fnm)
 484   static real taum      = 0.0, maxdist = 0.6;
 485   static int  nlevels   = 15;
 486   static int  nrestart  = 1;
 487   static int  maxframes = 100;
 488   static bool bFFT      = TRUE,bFit = TRUE, bVerbose = TRUE;
 489   t_pargs pa[] = {
 490     { "-taum",     FALSE, etREAL, &taum,
 491       "Rotational correlation time for your molecule. It is obligatory to pass this option" },
 492     { "-maxdist",  FALSE, etREAL, &maxdist,
 493       "Maximum distance to be plotted" },
 494     { "-nlevels",  FALSE, etINT,  &nlevels,
 495       "Number of levels for plotting" },
 496     { "-nframes", FALSE, etINT,  &maxframes,
 497       "Number of frames in your trajectory. Will stop analysis after this" },
 498     { "-fft",      FALSE, etBOOL, &bFFT,
 499       "Use FFT for correlation function" },
 500     { "-nrestart", FALSE, etINT,  &nrestart,
 501       "Number of frames between starting point for computation of ACF without FFT" },
 502     { "-fit",      FALSE, etBOOL, &bFit,
 503       "Do an optimal superposition on reference structure in tpx file" },
 504     { "-v",        FALSE, etBOOL, &bVerbose,
 505       "Tell you what I am about to do" }
 506   };
 507
 508   CopyRight(stderr,argv[0]);
 509   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 510                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 511   if (taum <= 0)
 512     fatal_error(0,"Please give me a sensible taum!\n");
 513   if (nlevels > 50) {
 514     nlevels = 50;
 515     fprintf(stderr,"Warning: too many levels, setting to %d\n",nlevels);
 516   }
 517
 518   top    = read_top(ftp2fn(efTPX,NFILE,fnm));
 519   natoms = top->atoms.nr;
 520   snew(xp,natoms);
 521   read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,NULL,box,
 522            &natoms,xp,NULL,NULL,NULL);
 523
 524   /* Determine the number of protons, and their index numbers
 525    * by checking the mass
 526    */
 527   nprot  = 0;
 528   snew(prot_ind,natoms);
 529   for(i=0; (i<natoms); i++)
 530     if (top->atoms.atom[i].m  < 2) {
 531       prot_ind[nprot++] = i;
 532     }
 533   fprintf(stderr,"There %d protons in your topology\n",nprot);
 534   snew(pair,(nprot*(nprot-1)/2));
 535   for(i=k=0; (i<nprot); i++) {
 536     for(j=i+1; (j<nprot); j++,k++) {
 537       pair[k].ai = prot_ind[i];
 538       pair[k].aj = prot_ind[j];
 539     }
 540   }
 541   sfree(prot_ind);
 542
 543   fprintf(stderr,"Select group for root least squares fit\n");
 544   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&ifit,&ind_fit,&gn_fit);
 545
 546   if (ifit < 3)
 547     fatal_error(0,"Need >= 3 points to fit!\n");
 548
 549   /* Make an array with weights for fitting */
 550   snew(w_rls,natoms);
 551   for(i=0; (i<ifit); i++)
 552     w_rls[ind_fit[i]]=top->atoms.atom[ind_fit[i]].m;
 553
 554   /* Prepare reference frame */
 555   snew(all_at,natoms);
 556   for(j=0; (j<natoms); j++)
 557     all_at[j]=j;
 558   rm_pbc(&(top->idef),natoms,box,xp,xp);
 559   reset_x(ifit,ind_fit,natoms,all_at,xp,w_rls);
 560   sfree(all_at);
 561
 562   spectrum(bVerbose,
 563            ftp2fn(efTRX,NFILE,fnm),ftp2fn(efDAT,NFILE,fnm),
 564            ftp2bSet(efDAT,NFILE,fnm),opt2fn("-corr",NFILE,fnm),
 565            opt2fn("-noe",NFILE,fnm),
 566            maxframes,bFFT,bFit,nrestart,
 567            k,pair,natoms,shifts,
 568            taum,maxdist,w_rls,xp,&(top->idef));
 569
 570   thanx(stdout);
 571
 572   return 0;
 573 }
 574