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