3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
70 real tauc,dtauc,S2,dS2;
75 void read_shifts(FILE *fp,int nx,real shiftx[],int ny,real shifty[])
80 for(i=0; (i<nx); i++) {
85 /*for(i=0; (i<ny); i++) {
92 complex c_sqr(complex c)
96 cs.re = c.re*c.re-c.im*c.im;
97 cs.im = 2.0*c.re*c.im;
102 complex c_add(complex c,complex d)
112 complex calc_ylm(int m,rvec rij,real r2,real r_3,real r_6)
114 static gmx_bool bFirst=TRUE;
115 static real y0,y1,y2;
116 real x,y,z,xsq,ysq,rxy,r1,cphi,sphi,cth,sth,fac;
120 y0 = sqrt(5/(4*M_PI));
121 y1 = sqrt(45/(24*M_PI));
122 y2 = sqrt(45/(96*M_PI));
140 /* Now calculate the spherical harmonics */
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;
156 cs.re = y0*(3*cth*cth-1);
166 void myfunc(real x,real a[],real *y,real dyda[],int na)
170 * y = a1 + (1-a1) exp(-a2 x)
172 * where in real life a1 is S^2 and a2 = 1/tauc
181 *y = S2 + (1-S2)*eee;
183 dyda[2] = x*tau1*tau1*(1-a[1])*eee;
186 void fit_one(gmx_bool bVerbose,
187 int nframes,real x[],real y[],real dy[],real ftol,
188 real *S2,real *dS2,real *tauc,real *dtauc)
190 void mrqmin(real x[],real y[],real sig[],int ndata,real a[],
191 int ma,int lista[],int mfit,real **covar,real **alpha,
193 void (*funcs)(real x,real a[],real *y,real dyda[],int na),
196 real *a,**covar,**alpha;
197 real chisq,ochisq,alamda;
199 int i,j,ma,mfit,*lista;
206 for(i=0; (i<ma+1); i++) {
212 a[1] = 0.99; /* S^2 */
213 a[2] = 0.1; /* tauc */
214 alamda = -1; /* Starting value */
219 mrqmin(x-1,y-1,dy-1,nframes,a,ma,lista,mfit,covar,alpha,
220 &chisq,myfunc,&alamda);
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]);
225 bCont = (((ochisq - chisq) > ftol*chisq) ||
226 ((ochisq == chisq)));
227 } while (bCont && (alamda != 0.0) && (j < 50));
229 fprintf(stderr,"\n");
231 /* Now get the covariance matrix out */
233 mrqmin(x-1,y-1,dy-1,nframes,a,ma,lista,mfit,covar,alpha,
234 &chisq,myfunc,&alamda);
237 *dS2 = sqrt(covar[1][1]);
239 *dtauc = sqrt(covar[2][2]);
241 for(i=0; (i<ma+1); i++) {
251 void calc_tauc(gmx_bool bVerbose,int npair,t_pair pair[],real dt,
252 int nframes,t_sij spec[],real **corr)
257 real S2,S22,tauc,fac;
263 for(i=0; (i<nframes); i++)
266 fprintf(stderr,"Fitting correlation function to Lipari&Szabo function\n");
267 fac=1.0/((real)nframes);
268 for(i=0; (i<npair); i++) {
270 /* Use Levenberg-Marquardt method to fit */
271 for(j=0; (j<nframes); j++)
274 nframes,x,corr[i],dy,ftol,
275 &(spec[i].S2),&(spec[i].dS2),
276 &(spec[i].tauc),&(spec[i].dtauc));
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));
290 void calc_aver(FILE *fp,int nframes,int npair,t_pair pair[],t_sij *spec,
297 md_6 = pow(maxdist,-6.0);
300 for(i=0; (i<npair); i++) {
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;
311 fprintf(fp," %8.3f+i%8.3f",c1.re,c1.im);
313 fprintf(fp," %8.3f-i%8.3f",c1.re,-c1.im);
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);
324 void plot_spectrum(char *noefn,int npair,t_pair pair[],t_sij *spec,real taum)
328 t_rgb rlo = { 1,0,0 },rhi = {1,1,1};
329 real Sijmax,Sijmin,pow6,pow3,pp3,pp6,ppy,tauc;
336 fp=xvgropen(noefn,"Cross Relaxation","Pair Index","\\8s\\4\\sij\\N");
337 for(i=0; (i<npair); i++) {
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;
342 Sijmax=max(Sijmax,sij.re);
343 Sijmin=min(Sijmin,sij.re);
344 fprintf(fp,"%5d %10g\n",i,sij.re);
347 fprintf(stderr,"Sijmin: %g, Sijmax: %g\n",Sijmin,Sijmax);
348 out=ffopen("spec.out","w");
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++) {
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);
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);
375 void spectrum(gmx_bool bVerbose,
376 char *trj,char *shifts,gmx_bool bAbInitio,
377 char *corrfn,char *noefn,
378 int maxframes,gmx_bool bFour,gmx_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)
384 int i,j,m,ii,jj,natoms,status,nframes;
388 real r2,r6,r_3,r_6,tauc;
392 gmx_rmpbc_t gpbc=NULL;
396 fprintf(stderr,"There is no kill like overkill! Going to malloc %d bytes\n",
397 npair*maxframes*sizeof(corr[0][0]));
399 for(i=0; (i<npair); i++)
400 snew(corr[i],maxframes);
402 natoms = read_first_x(&status,trj,&t0,&x,box);
404 gmx_fatal(FARGS,"Not enough atoms in trajectory");
405 gpbc = gmx_rmpbc_init(idef,ePBC,natoms,box);
408 if (nframes >= maxframes) {
409 fprintf(stderr,"\nThere are more than the %d frames you told me!",
415 fprintf(stderr,"\rframe: %d",nframes);
416 gmx_rmpbc(gpbc,box,x,x);
418 do_fit(natoms,w_rls,xp,x);
420 for(i=0; (i<npair); i++) {
423 rvec_sub(x[ii],x[jj],dx);
424 copy_rvec(dx,corr[i][nframes]);
428 r_3 = gmx_invsqrt(r6);
430 spec[i].rij_3 += r_3;
431 spec[i].rij_6 += r_6;
432 for(m=0; (m<5); m++) {
433 spec[i].Ylm[m] = c_add(spec[i].Ylm[m],
434 calc_ylm(m-2,dx,r2,r_3,r_6));
438 } while (read_next_x(status,&t,natoms,x,box));
441 fprintf(stderr,"\n");
443 gmx_rmpbc_done(gpbc);
445 fp=ffopen("ylm.out","w");
446 calc_aver(fp,nframes,npair,pair,spec,maxdist);
449 /* Select out the pairs that have to be correlated */
451 for(i=j=0; (i<npair); i++) {
453 Corr[j] = &(corr[i][0][0]);
457 fprintf(stderr,"There are %d NOEs in your simulation\n",j);
459 dt = (t1-t0)/(nframes-1);
462 do_autocorr(corrfn,"Correlation Function for Interproton Vectors",
463 nframes,j,Corr,dt,eacP2,nrestart,FALSE,FALSE,bFour,TRUE);
465 calc_tauc(bVerbose,npair,pair,dt,nframes/2,spec,(real **)corr);
467 plot_spectrum(noefn,npair,pair,spec,taum);
470 int gmx_relax(int argc,char *argv[])
472 const char *desc[] = {
473 "g_noe calculates a NOE spectrum"
478 int i,j,k,natoms,nprot,*prot_ind;
481 atom_id *ind_fit,*all_at;
490 { efTRX, "-f", NULL, ffREAD },
491 { efTPX, "-s", NULL, ffREAD },
492 { efNDX, NULL, NULL, ffREAD },
493 { efDAT, "-d", "shifts", ffREAD },
494 { efOUT, "-o","spec", ffWRITE },
495 { efXVG, "-corr", "rij-corr", ffWRITE },
496 { efXVG, "-noe", "noesy", ffWRITE }
498 #define NFILE asize(fnm)
499 static real taum = 0.0, maxdist = 0.6;
500 static int nlevels = 15;
501 static int nrestart = 1;
502 static int maxframes = 100;
503 static gmx_bool bFFT = TRUE,bFit = TRUE, bVerbose = TRUE;
505 { "-taum", FALSE, etREAL, &taum,
506 "Rotational correlation time for your molecule. It is obligatory to pass this option" },
507 { "-maxdist", FALSE, etREAL, &maxdist,
508 "Maximum distance to be plotted" },
509 { "-nlevels", FALSE, etINT, &nlevels,
510 "Number of levels for plotting" },
511 { "-nframes", FALSE, etINT, &maxframes,
512 "Number of frames in your trajectory. Will stop analysis after this" },
513 { "-fft", FALSE, etBOOL, &bFFT,
514 "Use FFT for correlation function" },
515 { "-nrestart", FALSE, etINT, &nrestart,
516 "Number of frames between starting point for computation of ACF without FFT" },
517 { "-fit", FALSE, etBOOL, &bFit,
518 "Do an optimal superposition on reference structure in tpx file" },
519 { "-v", FALSE, etBOOL, &bVerbose,
520 "Tell you what I am about to do" }
523 CopyRight(stderr,argv[0]);
524 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
525 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
527 gmx_fatal(FARGS,"Please give me a sensible taum!\n");
530 fprintf(stderr,"Warning: too many levels, setting to %d\n",nlevels);
533 top = read_top(ftp2fn(efTPX,NFILE,fnm));
534 natoms = top->atoms.nr;
536 read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,NULL,box,
537 &natoms,xp,NULL,NULL,NULL);
539 /* Determine the number of protons, and their index numbers
540 * by checking the mass
543 snew(prot_ind,natoms);
544 for(i=0; (i<natoms); i++)
545 if (top->atoms.atom[i].m < 2) {
546 prot_ind[nprot++] = i;
548 fprintf(stderr,"There %d protons in your topology\n",nprot);
549 snew(pair,(nprot*(nprot-1)/2));
550 for(i=k=0; (i<nprot); i++) {
551 for(j=i+1; (j<nprot); j++,k++) {
552 pair[k].ai = prot_ind[i];
553 pair[k].aj = prot_ind[j];
558 fprintf(stderr,"Select group for root least squares fit\n");
559 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&ifit,&ind_fit,&gn_fit);
562 gmx_fatal(FARGS,"Need >= 3 points to fit!\n");
564 /* Make an array with weights for fitting */
566 for(i=0; (i<ifit); i++)
567 w_rls[ind_fit[i]]=top->atoms.atom[ind_fit[i]].m;
569 /* Prepare reference frame */
571 for(j=0; (j<natoms); j++)
573 rm_pbc(&(top->idef),natoms,box,xp,xp);
574 reset_x(ifit,ind_fit,natoms,all_at,xp,w_rls);
578 ftp2fn(efTRX,NFILE,fnm),ftp2fn(efDAT,NFILE,fnm),
579 ftp2bSet(efDAT,NFILE,fnm),opt2fn("-corr",NFILE,fnm),
580 opt2fn("-noe",NFILE,fnm),
581 maxframes,bFFT,bFit,nrestart,
582 k,pair,natoms,shifts,
583 taum,maxdist,w_rls,xp,&(top->idef));