?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added Verlet scheme and NxN non-bonded functionality
[gromacs.git] / src / tools / gmx_disre.c
blob6a8935c66758b4503c446e2805190bfebb974606
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 #include <math.h>
39 #include <string.h>
40
41 #include "typedefs.h"
42 #include "macros.h"
43 #include "copyrite.h"
44 #include "mshift.h"
45 #include "xvgr.h"
46 #include "vec.h"
47 #include "do_fit.h"
48 #include "confio.h"
49 #include "smalloc.h"
50 #include "nrnb.h"
51 #include "disre.h"
52 #include "statutil.h"
53 #include "force.h"
54 #include "gstat.h"
55 #include "main.h"
56 #include "pdbio.h"
57 #include "index.h"
58 #include "mdatoms.h"
59 #include "tpxio.h"
60 #include "mdrun.h"
61 #include "names.h"
62 #include "matio.h"
63 #include "mtop_util.h"
64 #include "gmx_ana.h"
65
66
67 typedef struct {
68 int n;
69 real v;
70 } t_toppop;
71
72 t_toppop *top=NULL;
73 int ntop=0;
74
75 typedef struct {
76 int nv,nframes;
77 real sumv,averv,maxv;
78 real *aver1,*aver2,*aver_3,*aver_6;
79 } t_dr_result;
80
81 static void init5(int n)
82 {
83 ntop=n;
84 snew(top,ntop);
85 }
86
87 static void reset5(void)
88 {
89 int i;
90
91 for(i=0; (i<ntop); i++) {
92 top[i].n=-1;
93 top[i].v= 0;
94 }
95 }
96
97 int tpcomp(const void *a,const void *b)
98 {
99 t_toppop *tpa;
100 t_toppop *tpb;
101
102 tpa=(t_toppop *)a;
103 tpb=(t_toppop *)b;
104
105 return (1e7*(tpb->v-tpa->v));
106 }
107
108 static void add5(int ndr,real viol)
109 {
110 int i,mini;
111
112 mini=0;
113 for(i=1; (i<ntop); i++)
114 if (top[i].v < top[mini].v)
115 mini=i;
116 if (viol > top[mini].v) {
117 top[mini].v=viol;
118 top[mini].n=ndr;
119 }
120 }
121
122 static void print5(FILE *fp)
123 {
124 int i;
125
126 qsort(top,ntop,sizeof(top[0]),tpcomp);
127 fprintf(fp,"Index:");
128 for(i=0; (i<ntop); i++)
129 fprintf(fp," %6d",top[i].n);
130 fprintf(fp,"\nViol: ");
131 for(i=0; (i<ntop); i++)
132 fprintf(fp," %6.3f",top[i].v);
133 fprintf(fp,"\n");
134 }
135
136 static void check_viol(FILE *log,t_commrec *cr,
137 t_ilist *disres,t_iparams forceparams[],
138 t_functype functype[],rvec x[],rvec f[],
139 t_forcerec *fr,t_pbc *pbc,t_graph *g,t_dr_result dr[],
140 int clust_id,int isize,atom_id index[],real vvindex[],
141 t_fcdata *fcd)
142 {
143 t_iatom *forceatoms;
144 int i,j,nat,n,type,nviol,ndr,label;
145 real ener,rt,mviol,tviol,viol,lam,dvdl,drt;
146 static gmx_bool bFirst=TRUE;
147
148 lam =0;
149 dvdl =0;
150 tviol=0;
151 nviol=0;
152 mviol=0;
153 ndr=0;
154 if (ntop)
155 reset5();
156 forceatoms=disres->iatoms;
157 for(j=0; (j<isize); j++) {
158 vvindex[j]=0;
159 }
160 nat = interaction_function[F_DISRES].nratoms+1;
161 for(i=0; (i<disres->nr); ) {
162 type = forceatoms[i];
163 n = 0;
164 label = forceparams[type].disres.label;
165 if (debug)
166 fprintf(debug,"DISRE: ndr = %d, label = %d i=%d, n =%d\n",
167 ndr,label,i,n);
168 if (ndr != label)
169 gmx_fatal(FARGS,"tpr inconsistency. ndr = %d, label = %d\n",ndr,label);
170 do {
171 n += nat;
172 } while (((i+n) < disres->nr) &&
173 (forceparams[forceatoms[i+n]].disres.label == label));
174
175 calc_disres_R_6(cr->ms,n,&forceatoms[i],forceparams,
176 (const rvec*)x,pbc,fcd,NULL);
177
178 if (fcd->disres.Rt_6[0] <= 0)
179 gmx_fatal(FARGS,"ndr = %d, rt_6 = %f",ndr,fcd->disres.Rt_6[0]);
180
181 rt = pow(fcd->disres.Rt_6[0],-1.0/6.0);
182 dr[clust_id].aver1[ndr] += rt;
183 dr[clust_id].aver2[ndr] += sqr(rt);
184 drt = pow(rt,-3.0);
185 dr[clust_id].aver_3[ndr] += drt;
186 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
187
188 ener=interaction_function[F_DISRES].ifunc(n,&forceatoms[i],forceparams,
189 (const rvec*)x,f,fr->fshift,
190 pbc,g,lam,&dvdl,NULL,fcd,NULL);
191 viol = fcd->disres.sumviol;
192
193 if (viol > 0) {
194 nviol++;
195 if (ntop)
196 add5(forceparams[type].disres.label,viol);
197 if (viol > mviol)
198 mviol = viol;
199 tviol += viol;
200 for(j=0; (j<isize); j++) {
201 if (index[j] == forceparams[type].disres.label)
202 vvindex[j]=pow(fcd->disres.Rt_6[0],-1.0/6.0);
203 }
204 }
205 ndr ++;
206 i += n;
207 }
208 dr[clust_id].nv = nviol;
209 dr[clust_id].maxv = mviol;
210 dr[clust_id].sumv = tviol;
211 dr[clust_id].averv= tviol/ndr;
212 dr[clust_id].nframes++;
213
214 if (bFirst) {
215 fprintf(stderr,"\nThere are %d restraints and %d pairs\n",
216 ndr,disres->nr/nat);
217 bFirst = FALSE;
218 }
219 if (ntop)
220 print5(log);
221 }
222
223 typedef struct {
224 int index;
225 gmx_bool bCore;
226 real up1,r,rT3,rT6,viol,violT3,violT6;
227 } t_dr_stats;
228
229 static int drs_comp(const void *a,const void *b)
230 {
231 t_dr_stats *da,*db;
232
233 da = (t_dr_stats *)a;
234 db = (t_dr_stats *)b;
235
236 if (da->viol > db->viol)
237 return -1;
238 else if (da->viol < db->viol)
239 return 1;
240 else
241 return 0;
242 }
243
244 static void dump_dump(FILE *log,int ndr,t_dr_stats drs[])
245 {
246 static const char *core[] = { "All restraints", "Core restraints" };
247 static const char *tp[] = { "linear", "third power", "sixth power" };
248 real viol_tot,viol_max,viol=0;
249 gmx_bool bCore;
250 int nviol,nrestr;
251 int i,kkk;
252
253 for(bCore = FALSE; (bCore <= TRUE); bCore++) {
254 for(kkk=0; (kkk<3); kkk++) {
255 viol_tot = 0;
256 viol_max = 0;
257 nviol = 0;
258 nrestr = 0;
259 for(i=0; (i<ndr); i++) {
260 if (!bCore || (bCore && drs[i].bCore)) {
261 switch (kkk) {
262 case 0:
263 viol = drs[i].viol;
264 break;
265 case 1:
266 viol = drs[i].violT3;
267 break;
268 case 2:
269 viol = drs[i].violT6;
270 break;
271 default:
272 gmx_incons("Dumping violations");
273 }
274 viol_max = max(viol_max,viol);
275 if (viol > 0)
276 nviol++;
277 viol_tot += viol;
278 nrestr++;
279 }
280 }
281 if ((nrestr > 0) || (bCore && (nrestr < ndr))) {
282 fprintf(log,"\n");
283 fprintf(log,"+++++++ %s ++++++++\n",core[bCore]);
284 fprintf(log,"+++++++ Using %s averaging: ++++++++\n",tp[kkk]);
285 fprintf(log,"Sum of violations: %8.3f nm\n",viol_tot);
286 if (nrestr > 0)
287 fprintf(log,"Average violation: %8.3f nm\n",viol_tot/nrestr);
288 fprintf(log,"Largest violation: %8.3f nm\n",viol_max);
289 fprintf(log,"Number of violated restraints: %d/%d\n",nviol,nrestr);
290 }
291 }
292 }
293 }
294
295 static void dump_viol(FILE *log,int ndr,t_dr_stats *drs,gmx_bool bLinear)
296 {
297 int i;
298
299 fprintf(log,"Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
300 for(i=0; (i<ndr); i++) {
301 if (bLinear && (drs[i].viol == 0))
302 break;
303 fprintf(log,"%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
304 drs[i].index,yesno_names[drs[i].bCore],
305 drs[i].up1,drs[i].r,drs[i].rT3,drs[i].rT6,
306 drs[i].viol,drs[i].violT3,drs[i].violT6);
307 }
308 }
309
310 static gmx_bool is_core(int i,int isize,atom_id index[])
311 {
312 int kk;
313 gmx_bool bIC = FALSE;
314
315 for(kk=0; !bIC && (kk<isize); kk++)
316 bIC = (index[kk] == i);
317
318 return bIC;
319 }
320
321 static void dump_stats(FILE *log,int nsteps,int ndr,t_ilist *disres,
322 t_iparams ip[],t_dr_result *dr,
323 int isize,atom_id index[],t_atoms *atoms)
324 {
325 int i,j,nra;
326 t_dr_stats *drs;
327
328 fprintf(log,"\n");
329 fprintf(log,"++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
330 snew(drs,ndr);
331 j = 0;
332 nra = interaction_function[F_DISRES].nratoms+1;
333 for(i=0; (i<ndr); i++) {
334 /* Search for the appropriate restraint */
335 for( ; (j<disres->nr); j+=nra) {
336 if (ip[disres->iatoms[j]].disres.label == i)
337 break;
338 }
339 drs[i].index = i;
340 drs[i].bCore = is_core(i,isize,index);
341 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
342 drs[i].r = dr->aver1[i]/nsteps;
343 drs[i].rT3 = pow(dr->aver_3[i]/nsteps,-1.0/3.0);
344 drs[i].rT6 = pow(dr->aver_6[i]/nsteps,-1.0/6.0);
345 drs[i].viol = max(0,drs[i].r-drs[i].up1);
346 drs[i].violT3 = max(0,drs[i].rT3-drs[i].up1);
347 drs[i].violT6 = max(0,drs[i].rT6-drs[i].up1);
348 if (atoms) {
349 int j1 = disres->iatoms[j+1];
350 int j2 = disres->iatoms[j+2];
351 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
352 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
353 }
354 }
355 dump_viol(log,ndr,drs,FALSE);
356
357 fprintf(log,"+++ Sorted by linear averaged violations: +++\n");
358 qsort(drs,ndr,sizeof(drs[0]),drs_comp);
359 dump_viol(log,ndr,drs,TRUE);
360
361 dump_dump(log,ndr,drs);
362
363 sfree(drs);
364 }
365
366 static void dump_clust_stats(FILE *fp,int ndr,t_ilist *disres,
367 t_iparams ip[],t_blocka *clust,t_dr_result dr[],
368 char *clust_name[],int isize,atom_id index[])
369 {
370 int i,j,k,nra,mmm=0;
371 double sumV,maxV,sumVT3,sumVT6,maxVT3,maxVT6;
372 t_dr_stats *drs;
373
374 fprintf(fp,"\n");
375 fprintf(fp,"++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
376 fprintf(fp,"Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
377
378 snew(drs,ndr);
379
380 for(k=0; (k<clust->nr); k++) {
381 if (dr[k].nframes == 0)
382 continue;
383 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
384 gmx_fatal(FARGS,"Inconsistency in cluster %s.\n"
385 "Found %d frames in trajectory rather than the expected %d\n",
386 clust_name[k],dr[k].nframes,
387 clust->index[k+1]-clust->index[k]);
388 if (!clust_name[k])
389 gmx_fatal(FARGS,"Inconsistency with cluster %d. Invalid name",k);
390 j = 0;
391 nra = interaction_function[F_DISRES].nratoms+1;
392 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
393 for(i=0; (i<ndr); i++) {
394 /* Search for the appropriate restraint */
395 for( ; (j<disres->nr); j+=nra) {
396 if (ip[disres->iatoms[j]].disres.label == i)
397 break;
398 }
399 drs[i].index = i;
400 drs[i].bCore = is_core(i,isize,index);
401 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
402 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
403 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
404 gmx_fatal(FARGS,"dr[%d].aver_3[%d] = %f",k,i,dr[k].aver_3[i]);
405 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes,-1.0/3.0);
406 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes,-1.0/6.0);
407 drs[i].viol = max(0,drs[i].r-drs[i].up1);
408 drs[i].violT3 = max(0,drs[i].rT3-drs[i].up1);
409 drs[i].violT6 = max(0,drs[i].rT6-drs[i].up1);
410 sumV += drs[i].viol;
411 sumVT3 += drs[i].violT3;
412 sumVT6 += drs[i].violT6;
413 maxV = max(maxV,drs[i].viol);
414 maxVT3 = max(maxVT3,drs[i].violT3);
415 maxVT6 = max(maxVT6,drs[i].violT6);
416 }
417 if (strcmp(clust_name[k],"1000") == 0) {
418 mmm ++;
419 }
420 fprintf(fp,"%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
421 clust_name[k],
422 dr[k].nframes,sumV,maxV,sumVT3,maxVT3,sumVT6,maxVT6);
423
424 }
425 fflush(fp);
426 sfree(drs);
427 }
428
429 static void init_dr_res(t_dr_result *dr,int ndr)
430 {
431 snew(dr->aver1,ndr+1);
432 snew(dr->aver2,ndr+1);
433 snew(dr->aver_3,ndr+1);
434 snew(dr->aver_6,ndr+1);
435 dr->nv = 0;
436 dr->nframes = 0;
437 dr->sumv = 0;
438 dr->maxv = 0;
439 dr->averv = 0;
440 }
441
442 static void dump_disre_matrix(const char *fn,t_dr_result *dr,int ndr,
443 int nsteps,t_idef *idef,gmx_mtop_t *mtop,
444 real max_dr,int nlevels,gmx_bool bThird)
445 {
446 FILE *fp;
447 int *resnr;
448 int n_res,a_offset,mb,mol,a;
449 t_atoms *atoms;
450 int iii,i,j,nra,nratoms,tp,ri,rj,index,nlabel,label;
451 atom_id ai,aj,*ptr;
452 real **matrix,*t_res,hi,*w_dr,rav,rviol;
453 t_rgb rlo = { 1, 1, 1 };
454 t_rgb rhi = { 0, 0, 0 };
455 if (fn == NULL)
456 return;
457
458 snew(resnr,mtop->natoms);
459 n_res = 0;
460 a_offset = 0;
461 for(mb=0; mb<mtop->nmolblock; mb++) {
462 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
463 for(mol=0; mol<mtop->molblock[mb].nmol; mol++) {
464 for(a=0; a<atoms->nr; a++) {
465 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
466 }
467 n_res += atoms->nres;
468 a_offset += atoms->nr;
469 }
470 }
471
472 snew(t_res,n_res);
473 for(i=0; (i<n_res); i++)
474 t_res[i] = i+1;
475 snew(matrix,n_res);
476 for(i=0; (i<n_res); i++)
477 snew(matrix[i],n_res);
478 nratoms = interaction_function[F_DISRES].nratoms;
479 nra = (idef->il[F_DISRES].nr/(nratoms+1));
480 snew(ptr,nra+1);
481 index = 0;
482 nlabel = 0;
483 ptr[0] = 0;
484 snew(w_dr,ndr);
485 for(i=0; (i<idef->il[F_DISRES].nr); i+=nratoms+1) {
486 tp = idef->il[F_DISRES].iatoms[i];
487 label = idef->iparams[tp].disres.label;
488
489 if (label != index) {
490 /* Set index pointer */
491 ptr[index+1] = i;
492 if (nlabel <= 0)
493 gmx_fatal(FARGS,"nlabel is %d, label = %d",nlabel,label);
494 if (index >= ndr)
495 gmx_fatal(FARGS,"ndr = %d, index = %d",ndr,index);
496 /* Update the weight */
497 w_dr[index] = 1.0/nlabel;
498 index = label;
499 nlabel = 1;
500 }
501 else {
502 nlabel++;
503 }
504 }
505 printf("nlabel = %d, index = %d, ndr = %d\n",nlabel,index,ndr);
506 hi = 0;
507 for(i=0; (i<ndr); i++) {
508 for(j=ptr[i]; (j<ptr[i+1]); j+=nratoms+1) {
509 tp = idef->il[F_DISRES].iatoms[j];
510 ai = idef->il[F_DISRES].iatoms[j+1];
511 aj = idef->il[F_DISRES].iatoms[j+2];
512
513 ri = resnr[ai];
514 rj = resnr[aj];
515 if (bThird)
516 rav = pow(dr->aver_3[i]/nsteps,-1.0/3.0);
517 else
518 rav = dr->aver1[i]/nsteps;
519 if (debug)
520 fprintf(debug,"DR %d, atoms %d, %d, distance %g\n",i,ai,aj,rav);
521 rviol = max(0,rav-idef->iparams[tp].disres.up1);
522 matrix[ri][rj] += w_dr[i]*rviol;
523 matrix[rj][ri] += w_dr[i]*rviol;
524 hi = max(hi,matrix[ri][rj]);
525 hi = max(hi,matrix[rj][ri]);
526 }
527 }
528
529 sfree(resnr);
530
531 if (max_dr > 0) {
532 if (hi > max_dr)
533 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n",max_dr,hi);
534 hi = max_dr;
535 }
536 printf("Highest level in the matrix will be %g\n",hi);
537 fp = ffopen(fn,"w");
538 write_xpm(fp,0,"Distance Violations","<V> (nm)","Residue","Residue",
539 n_res,n_res,t_res,t_res,matrix,0,hi,rlo,rhi,&nlevels);
540 ffclose(fp);
541 }
542
543 int gmx_disre(int argc,char *argv[])
544 {
545 const char *desc[] = {
546 "[TT]g_disre[tt] computes violations of distance restraints.",
547 "If necessary, all protons can be added to a protein molecule ",
548 "using the [TT]g_protonate[tt] program.[PAR]",
549 "The program always",
550 "computes the instantaneous violations rather than time-averaged,",
551 "because this analysis is done from a trajectory file afterwards",
552 "it does not make sense to use time averaging. However,",
553 "the time averaged values per restraint are given in the log file.[PAR]",
554 "An index file may be used to select specific restraints for",
555 "printing.[PAR]",
556 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
557 "amount of average violations.[PAR]",
558 "When the [TT]-c[tt] option is given, an index file will be read",
559 "containing the frames in your trajectory corresponding to the clusters",
560 "(defined in another manner) that you want to analyze. For these clusters",
561 "the program will compute average violations using the third power",
562 "averaging algorithm and print them in the log file."
563 };
564 static int ntop = 0;
565 static int nlevels = 20;
566 static real max_dr = 0;
567 static gmx_bool bThird = TRUE;
568 t_pargs pa[] = {
569 { "-ntop", FALSE, etINT, {&ntop},
570 "Number of large violations that are stored in the log file every step" },
571 { "-maxdr", FALSE, etREAL, {&max_dr},
572 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
573 { "-nlevels", FALSE, etINT, {&nlevels},
574 "Number of levels in the matrix output" },
575 { "-third", FALSE, etBOOL, {&bThird},
576 "Use inverse third power averaging or linear for matrix output" }
577 };
578
579 FILE *out=NULL,*aver=NULL,*numv=NULL,*maxxv=NULL,*xvg=NULL;
580 t_tpxheader header;
581 t_inputrec ir;
582 gmx_mtop_t mtop;
583 rvec *xtop;
584 gmx_localtop_t *top;
585 t_atoms *atoms=NULL;
586 t_forcerec *fr;
587 t_fcdata fcd;
588 t_nrnb nrnb;
589 t_commrec *cr;
590 t_graph *g;
591 int ntopatoms,natoms,i,j,kkk;
592 t_trxstatus *status;
593 real t;
594 rvec *x,*f,*xav=NULL;
595 matrix box;
596 gmx_bool bPDB;
597 int isize;
598 atom_id *index=NULL,*ind_fit=NULL;
599 char *grpname;
600 t_cluster_ndx *clust=NULL;
601 t_dr_result dr,*dr_clust=NULL;
602 char **leg;
603 real *vvindex=NULL,*w_rls=NULL;
604 t_mdatoms *mdatoms;
605 t_pbc pbc,*pbc_null;
606 int my_clust;
607 FILE *fplog;
608 output_env_t oenv;
609 gmx_rmpbc_t gpbc=NULL;
610
611 t_filenm fnm[] = {
612 { efTPX, NULL, NULL, ffREAD },
613 { efTRX, "-f", NULL, ffREAD },
614 { efXVG, "-ds", "drsum", ffWRITE },
615 { efXVG, "-da", "draver", ffWRITE },
616 { efXVG, "-dn", "drnum", ffWRITE },
617 { efXVG, "-dm", "drmax", ffWRITE },
618 { efXVG, "-dr", "restr", ffWRITE },
619 { efLOG, "-l", "disres", ffWRITE },
620 { efNDX, NULL, "viol", ffOPTRD },
621 { efPDB, "-q", "viol", ffOPTWR },
622 { efNDX, "-c", "clust", ffOPTRD },
623 { efXPM, "-x", "matrix", ffOPTWR }
624 };
625 #define NFILE asize(fnm)
626
627 cr = init_par(&argc,&argv);
628 CopyRight(stderr,argv[0]);
629 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
630 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
631
632 gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,FALSE,0,&fplog);
633
634 if (ntop)
635 init5(ntop);
636
637 read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header,FALSE,NULL,NULL);
638 snew(xtop,header.natoms);
639 read_tpx(ftp2fn(efTPX,NFILE,fnm),&ir,box,&ntopatoms,xtop,NULL,NULL,&mtop);
640 bPDB = opt2bSet("-q",NFILE,fnm);
641 if (bPDB) {
642 snew(xav,ntopatoms);
643 snew(ind_fit,ntopatoms);
644 snew(w_rls,ntopatoms);
645 for(kkk=0; (kkk<ntopatoms); kkk++) {
646 w_rls[kkk] = 1;
647 ind_fit[kkk] = kkk;
648 }
649
650 snew(atoms,1);
651 *atoms = gmx_mtop_global_atoms(&mtop);
652
653 if (atoms->pdbinfo == NULL) {
654 snew(atoms->pdbinfo,atoms->nr);
655 }
656 }
657
658 top = gmx_mtop_generate_local_top(&mtop,&ir);
659
660 g = NULL;
661 pbc_null = NULL;
662 if (ir.ePBC != epbcNONE) {
663 if (ir.bPeriodicMols)
664 pbc_null = &pbc;
665 else
666 g = mk_graph(fplog,&top->idef,0,mtop.natoms,FALSE,FALSE);
667 }
668
669 if (ftp2bSet(efNDX,NFILE,fnm)) {
670 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
671 xvg=xvgropen(opt2fn("-dr",NFILE,fnm),"Individual Restraints","Time (ps)",
672 "nm",oenv);
673 snew(vvindex,isize);
674 snew(leg,isize);
675 for(i=0; (i<isize); i++) {
676 index[i]++;
677 snew(leg[i],12);
678 sprintf(leg[i],"index %d",index[i]);
679 }
680 xvgr_legend(xvg,isize,(const char**)leg,oenv);
681 }
682 else
683 isize=0;
684
685 ir.dr_tau=0.0;
686 init_disres(fplog,&mtop,&ir,NULL,FALSE,&fcd,NULL);
687
688 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
689 snew(f,5*natoms);
690
691 init_dr_res(&dr,fcd.disres.nres);
692 if (opt2bSet("-c",NFILE,fnm)) {
693 clust = cluster_index(fplog,opt2fn("-c",NFILE,fnm));
694 snew(dr_clust,clust->clust->nr+1);
695 for(i=0; (i<=clust->clust->nr); i++)
696 init_dr_res(&dr_clust[i],fcd.disres.nres);
697 }
698 else {
699 out =xvgropen(opt2fn("-ds",NFILE,fnm),
700 "Sum of Violations","Time (ps)","nm",oenv);
701 aver=xvgropen(opt2fn("-da",NFILE,fnm),
702 "Average Violation","Time (ps)","nm",oenv);
703 numv=xvgropen(opt2fn("-dn",NFILE,fnm),
704 "# Violations","Time (ps)","#",oenv);
705 maxxv=xvgropen(opt2fn("-dm",NFILE,fnm),
706 "Largest Violation","Time (ps)","nm",oenv);
707 }
708
709 mdatoms = init_mdatoms(fplog,&mtop,ir.efep!=efepNO);
710 atoms2md(&mtop,&ir,0,NULL,0,mtop.natoms,mdatoms);
711 update_mdatoms(mdatoms,ir.fepvals->init_lambda);
712 fr = mk_forcerec();
713 fprintf(fplog,"Made forcerec\n");
714 init_forcerec(fplog,oenv,fr,NULL,&ir,&mtop,cr,box,FALSE,
715 NULL,NULL,NULL,NULL,NULL,FALSE,-1);
716 init_nrnb(&nrnb);
717 if (ir.ePBC != epbcNONE)
718 gpbc = gmx_rmpbc_init(&top->idef,ir.ePBC,natoms,box);
719
720 j=0;
721 do {
722 if (ir.ePBC != epbcNONE) {
723 if (ir.bPeriodicMols)
724 set_pbc(&pbc,ir.ePBC,box);
725 else
726 gmx_rmpbc(gpbc,natoms,box,x);
727 }
728
729 if (clust) {
730 if (j > clust->maxframe)
731 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file. t = %8f\n",t);
732 my_clust = clust->inv_clust[j];
733 range_check(my_clust,0,clust->clust->nr);
734 check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
735 top->idef.iparams,top->idef.functype,
736 x,f,fr,pbc_null,g,dr_clust,my_clust,isize,index,vvindex,&fcd);
737 }
738 else
739 check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
740 top->idef.iparams,top->idef.functype,
741 x,f,fr,pbc_null,g,&dr,0,isize,index,vvindex,&fcd);
742 if (bPDB) {
743 reset_x(atoms->nr,ind_fit,atoms->nr,NULL,x,w_rls);
744 do_fit(atoms->nr,w_rls,x,x);
745 if (j == 0) {
746 /* Store the first frame of the trajectory as 'characteristic'
747 * for colouring with violations.
748 */
749 for(kkk=0; (kkk<atoms->nr); kkk++)
750 copy_rvec(x[kkk],xav[kkk]);
751 }
752 }
753 if (!clust) {
754 if (isize > 0) {
755 fprintf(xvg,"%10g",t);
756 for(i=0; (i<isize); i++)
757 fprintf(xvg," %10g",vvindex[i]);
758 fprintf(xvg,"\n");
759 }
760 fprintf(out, "%10g %10g\n",t,dr.sumv);
761 fprintf(aver, "%10g %10g\n",t,dr.averv);
762 fprintf(maxxv,"%10g %10g\n",t,dr.maxv);
763 fprintf(numv, "%10g %10d\n",t,dr.nv);
764 }
765 j++;
766 } while (read_next_x(oenv,status,&t,natoms,x,box));
767 close_trj(status);
768 if (ir.ePBC != epbcNONE)
769 gmx_rmpbc_done(gpbc);
770
771 if (clust) {
772 dump_clust_stats(fplog,fcd.disres.nres,&(top->idef.il[F_DISRES]),
773 top->idef.iparams,clust->clust,dr_clust,
774 clust->grpname,isize,index);
775 }
776 else {
777 dump_stats(fplog,j,fcd.disres.nres,&(top->idef.il[F_DISRES]),
778 top->idef.iparams,&dr,isize,index,
779 bPDB ? atoms : NULL);
780 if (bPDB) {
781 write_sto_conf(opt2fn("-q",NFILE,fnm),
782 "Coloured by average violation in Angstrom",
783 atoms,xav,NULL,ir.ePBC,box);
784 }
785 dump_disre_matrix(opt2fn_null("-x",NFILE,fnm),&dr,fcd.disres.nres,
786 j,&top->idef,&mtop,max_dr,nlevels,bThird);
787 ffclose(out);
788 ffclose(aver);
789 ffclose(numv);
790 ffclose(maxxv);
791 if (isize > 0) {
792 ffclose(xvg);
793 do_view(oenv,opt2fn("-dr",NFILE,fnm),"-nxy");
794 }
795 do_view(oenv,opt2fn("-dn",NFILE,fnm),"-nxy");
796 do_view(oenv,opt2fn("-da",NFILE,fnm),"-nxy");
797 do_view(oenv,opt2fn("-ds",NFILE,fnm),"-nxy");
798 do_view(oenv,opt2fn("-dm",NFILE,fnm),"-nxy");
799 }
800 thanx(stderr);
801
802 gmx_finalize_par();
803
804 gmx_log_close(fplog);
805
806 return 0;
807 }
The ultimate molecular dynamics simulation package