?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / gmx_disre.c
blob47d9335cf69fa4eb79d4bf928739621cd20d0655
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 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 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 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,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 bool is_core(int i,int isize,atom_id index[])
311 {
312 int kk;
313 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,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 "g_disre computes violations of distance restraints.",
547 "If necessary all protons can be added to a protein molecule ",
548 "using the protonate 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 pdb 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 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 status,ntopatoms,natoms,i,j,kkk;
592 real t;
593 rvec *x,*f,*xav=NULL;
594 matrix box;
595 bool bPDB;
596 int isize;
597 atom_id *index=NULL,*ind_fit=NULL;
598 char *grpname;
599 t_cluster_ndx *clust=NULL;
600 t_dr_result dr,*dr_clust=NULL;
601 char **leg;
602 real *vvindex=NULL,*w_rls=NULL;
603 t_mdatoms *mdatoms;
604 t_pbc pbc,*pbc_null;
605 int my_clust;
606 FILE *fplog;
607 output_env_t oenv;
608
609 t_filenm fnm[] = {
610 { efTPX, NULL, NULL, ffREAD },
611 { efTRX, "-f", NULL, ffREAD },
612 { efXVG, "-ds", "drsum", ffWRITE },
613 { efXVG, "-da", "draver", ffWRITE },
614 { efXVG, "-dn", "drnum", ffWRITE },
615 { efXVG, "-dm", "drmax", ffWRITE },
616 { efXVG, "-dr", "restr", ffWRITE },
617 { efLOG, "-l", "disres", ffWRITE },
618 { efNDX, NULL, "viol", ffOPTRD },
619 { efPDB, "-q", "viol", ffOPTWR },
620 { efNDX, "-c", "clust", ffOPTRD },
621 { efXPM, "-x", "matrix", ffOPTWR }
622 };
623 #define NFILE asize(fnm)
624
625 cr = init_par(&argc,&argv);
626 CopyRight(stderr,argv[0]);
627 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
628 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
629
630 gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,FALSE,0,&fplog);
631
632 if (ntop)
633 init5(ntop);
634
635 read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header,FALSE,NULL,NULL);
636 snew(xtop,header.natoms);
637 read_tpx(ftp2fn(efTPX,NFILE,fnm),&ir,box,&ntopatoms,xtop,NULL,NULL,&mtop);
638 bPDB = opt2bSet("-q",NFILE,fnm);
639 if (bPDB) {
640 snew(xav,ntopatoms);
641 snew(ind_fit,ntopatoms);
642 snew(w_rls,ntopatoms);
643 for(kkk=0; (kkk<ntopatoms); kkk++) {
644 w_rls[kkk] = 1;
645 ind_fit[kkk] = kkk;
646 }
647
648 snew(atoms,1);
649 *atoms = gmx_mtop_global_atoms(&mtop);
650
651 if (atoms->pdbinfo == NULL) {
652 snew(atoms->pdbinfo,atoms->nr);
653 }
654 }
655
656 top = gmx_mtop_generate_local_top(&mtop,&ir);
657
658 g = NULL;
659 pbc_null = NULL;
660 if (ir.ePBC != epbcNONE) {
661 if (ir.bPeriodicMols)
662 pbc_null = &pbc;
663 else
664 g = mk_graph(fplog,&top->idef,0,mtop.natoms,FALSE,FALSE);
665 }
666
667 if (ftp2bSet(efNDX,NFILE,fnm)) {
668 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
669 xvg=xvgropen(opt2fn("-dr",NFILE,fnm),"Inidividual Restraints","Time (ps)",
670 "nm",oenv);
671 snew(vvindex,isize);
672 snew(leg,isize);
673 for(i=0; (i<isize); i++) {
674 index[i]++;
675 snew(leg[i],12);
676 sprintf(leg[i],"index %d",index[i]);
677 }
678 xvgr_legend(xvg,isize,leg,oenv);
679 }
680 else
681 isize=0;
682
683 ir.dr_tau=0.0;
684 init_disres(fplog,&mtop,&ir,NULL,FALSE,&fcd,NULL);
685
686 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
687 snew(f,5*natoms);
688
689 init_dr_res(&dr,fcd.disres.nres);
690 if (opt2bSet("-c",NFILE,fnm)) {
691 clust = cluster_index(fplog,opt2fn("-c",NFILE,fnm));
692 snew(dr_clust,clust->clust->nr+1);
693 for(i=0; (i<=clust->clust->nr); i++)
694 init_dr_res(&dr_clust[i],fcd.disres.nres);
695 }
696 else {
697 out =xvgropen(opt2fn("-ds",NFILE,fnm),
698 "Sum of Violations","Time (ps)","nm",oenv);
699 aver=xvgropen(opt2fn("-da",NFILE,fnm),
700 "Average Violation","Time (ps)","nm",oenv);
701 numv=xvgropen(opt2fn("-dn",NFILE,fnm),
702 "# Violations","Time (ps)","#",oenv);
703 maxxv=xvgropen(opt2fn("-dm",NFILE,fnm),
704 "Largest Violation","Time (ps)","nm",oenv);
705 }
706
707 mdatoms = init_mdatoms(fplog,&mtop,ir.efep!=efepNO);
708 atoms2md(&mtop,&ir,0,NULL,0,mtop.natoms,mdatoms);
709 update_mdatoms(mdatoms,ir.init_lambda);
710 fr = mk_forcerec();
711 fprintf(fplog,"Made forcerec\n");
712 init_forcerec(fplog,oenv,fr,NULL,&ir,&mtop,cr,box,FALSE,NULL,NULL,NULL,
713 FALSE,-1);
714 init_nrnb(&nrnb);
715 j=0;
716 do {
717 if (ir.ePBC != epbcNONE) {
718 if (ir.bPeriodicMols)
719 set_pbc(&pbc,ir.ePBC,box);
720 else
721 rm_pbc(&top->idef,ir.ePBC,natoms,box,x,x);
722 }
723
724 if (clust) {
725 if (j > clust->maxframe)
726 gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file. t = %8f\n",t);
727 my_clust = clust->inv_clust[j];
728 range_check(my_clust,0,clust->clust->nr);
729 check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
730 top->idef.iparams,top->idef.functype,
731 x,f,fr,pbc_null,g,dr_clust,my_clust,isize,index,vvindex,&fcd);
732 }
733 else
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,0,isize,index,vvindex,&fcd);
737 if (bPDB) {
738 reset_x(atoms->nr,ind_fit,atoms->nr,NULL,x,w_rls);
739 do_fit(atoms->nr,w_rls,x,x);
740 if (j == 0) {
741 /* Store the first frame of the trajectory as 'characteristic'
742 * for colouring with violations.
743 */
744 for(kkk=0; (kkk<atoms->nr); kkk++)
745 copy_rvec(x[kkk],xav[kkk]);
746 }
747 }
748 if (!clust) {
749 if (isize > 0) {
750 fprintf(xvg,"%10g",t);
751 for(i=0; (i<isize); i++)
752 fprintf(xvg," %10g",vvindex[i]);
753 fprintf(xvg,"\n");
754 }
755 fprintf(out, "%10g %10g\n",t,dr.sumv);
756 fprintf(aver, "%10g %10g\n",t,dr.averv);
757 fprintf(maxxv,"%10g %10g\n",t,dr.maxv);
758 fprintf(numv, "%10g %10d\n",t,dr.nv);
759 }
760 j++;
761 } while (read_next_x(oenv,status,&t,natoms,x,box));
762 close_trj(status);
763
764 if (clust) {
765 dump_clust_stats(fplog,fcd.disres.nres,&(top->idef.il[F_DISRES]),
766 top->idef.iparams,clust->clust,dr_clust,
767 clust->grpname,isize,index);
768 }
769 else {
770 dump_stats(fplog,j,fcd.disres.nres,&(top->idef.il[F_DISRES]),
771 top->idef.iparams,&dr,isize,index,
772 bPDB ? atoms : NULL);
773 if (bPDB) {
774 write_sto_conf(opt2fn("-q",NFILE,fnm),
775 "Coloured by average violation in Angstrom",
776 atoms,xav,NULL,ir.ePBC,box);
777 }
778 dump_disre_matrix(opt2fn_null("-x",NFILE,fnm),&dr,fcd.disres.nres,
779 j,&top->idef,&mtop,max_dr,nlevels,bThird);
780 ffclose(out);
781 ffclose(aver);
782 ffclose(numv);
783 ffclose(maxxv);
784 if (isize > 0) {
785 ffclose(xvg);
786 do_view(oenv,opt2fn("-dr",NFILE,fnm),"-nxy");
787 }
788 do_view(oenv,opt2fn("-dn",NFILE,fnm),"-nxy");
789 do_view(oenv,opt2fn("-da",NFILE,fnm),"-nxy");
790 do_view(oenv,opt2fn("-ds",NFILE,fnm),"-nxy");
791 do_view(oenv,opt2fn("-dm",NFILE,fnm),"-nxy");
792 }
793 thanx(stderr);
794
795 if (gmx_parallel_env_initialized())
796 gmx_finalize();
797
798 gmx_log_close(fplog);
799
800 return 0;
801 }