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