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
63 #include "mtop_util.h"
78 real
*aver1
,*aver2
,*aver_3
,*aver_6
;
81 static void init5(int n
)
87 static void reset5(void)
91 for(i
=0; (i
<ntop
); i
++) {
97 int tpcomp(const void *a
,const void *b
)
105 return (1e7
*(tpb
->v
-tpa
->v
));
108 static void add5(int ndr
,real viol
)
113 for(i
=1; (i
<ntop
); i
++)
114 if (top
[i
].v
< top
[mini
].v
)
116 if (viol
> top
[mini
].v
) {
122 static void print5(FILE *fp
)
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
);
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
[],
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
;
156 forceatoms
=disres
->iatoms
;
157 for(j
=0; (j
<isize
); j
++) {
160 nat
= interaction_function
[F_DISRES
].nratoms
+1;
161 for(i
=0; (i
<disres
->nr
); ) {
162 type
= forceatoms
[i
];
164 label
= forceparams
[type
].disres
.label
;
166 fprintf(debug
,"DISRE: ndr = %d, label = %d i=%d, n =%d\n",
169 gmx_fatal(FARGS
,"tpr inconsistency. ndr = %d, label = %d\n",ndr
,label
);
172 } while (((i
+n
) < disres
->nr
) &&
173 (forceparams
[forceatoms
[i
+n
]].disres
.label
== label
));
175 calc_disres_R_6(cr
->ms
,n
,&forceatoms
[i
],forceparams
,
176 (const rvec
*)x
,pbc
,fcd
,NULL
);
178 if (fcd
->disres
.Rt_6
[0] <= 0)
179 gmx_fatal(FARGS
,"ndr = %d, rt_6 = %f",ndr
,fcd
->disres
.Rt_6
[0]);
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
);
185 dr
[clust_id
].aver_3
[ndr
] += drt
;
186 dr
[clust_id
].aver_6
[ndr
] += fcd
->disres
.Rt_6
[0];
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
;
196 add5(forceparams
[type
].disres
.label
,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);
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
++;
215 fprintf(stderr
,"\nThere are %d restraints and %d pairs\n",
226 real up1
,r
,rT3
,rT6
,viol
,violT3
,violT6
;
229 static int drs_comp(const void *a
,const void *b
)
233 da
= (t_dr_stats
*)a
;
234 db
= (t_dr_stats
*)b
;
236 if (da
->viol
> db
->viol
)
238 else if (da
->viol
< db
->viol
)
244 static void dump_dump(FILE *log
,int ndr
,t_dr_stats drs
[])
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;
253 for(bCore
= FALSE
; (bCore
<= TRUE
); bCore
++) {
254 for(kkk
=0; (kkk
<3); kkk
++) {
259 for(i
=0; (i
<ndr
); i
++) {
260 if (!bCore
|| (bCore
&& drs
[i
].bCore
)) {
266 viol
= drs
[i
].violT3
;
269 viol
= drs
[i
].violT6
;
272 gmx_incons("Dumping violations");
274 viol_max
= max(viol_max
,viol
);
281 if ((nrestr
> 0) || (bCore
&& (nrestr
< ndr
))) {
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
);
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
);
295 static void dump_viol(FILE *log
,int ndr
,t_dr_stats
*drs
,gmx_bool bLinear
)
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))
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
);
310 static gmx_bool
is_core(int i
,int isize
,atom_id index
[])
313 gmx_bool bIC
= FALSE
;
315 for(kk
=0; !bIC
&& (kk
<isize
); kk
++)
316 bIC
= (index
[kk
] == i
);
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
)
329 fprintf(log
,"++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
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
)
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
);
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;
355 dump_viol(log
,ndr
,drs
,FALSE
);
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
);
361 dump_dump(log
,ndr
,drs
);
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
[])
371 double sumV
,maxV
,sumVT3
,sumVT6
,maxVT3
,maxVT6
;
375 fprintf(fp
,"++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
376 fprintf(fp
,"Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
380 for(k
=0; (k
<clust
->nr
); k
++) {
381 if (dr
[k
].nframes
== 0)
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
]);
389 gmx_fatal(FARGS
,"Inconsistency with cluster %d. Invalid name",k
);
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
)
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
);
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
);
417 if (strcmp(clust_name
[k
],"1000") == 0) {
420 fprintf(fp
,"%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
422 dr
[k
].nframes
,sumV
,maxV
,sumVT3
,maxVT3
,sumVT6
,maxVT6
);
429 static void init_dr_res(t_dr_result
*dr
,int ndr
)
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);
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
)
448 int n_res
,a_offset
,mb
,mol
,a
;
450 int iii
,i
,j
,nra
,nratoms
,tp
,ri
,rj
,index
,nlabel
,label
;
452 real
**matrix
,*t_res
,hi
,*w_dr
,rav
,rviol
;
453 t_rgb rlo
= { 1, 1, 1 };
454 t_rgb rhi
= { 0, 0, 0 };
458 snew(resnr
,mtop
->natoms
);
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
;
467 n_res
+= atoms
->nres
;
468 a_offset
+= atoms
->nr
;
473 for(i
=0; (i
<n_res
); i
++)
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));
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
;
489 if (label
!= index
) {
490 /* Set index pointer */
493 gmx_fatal(FARGS
,"nlabel is %d, label = %d",nlabel
,label
);
495 gmx_fatal(FARGS
,"ndr = %d, index = %d",ndr
,index
);
496 /* Update the weight */
497 w_dr
[index
] = 1.0/nlabel
;
505 printf("nlabel = %d, index = %d, ndr = %d\n",nlabel
,index
,ndr
);
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];
516 rav
= pow(dr
->aver_3
[i
]/nsteps
,-1.0/3.0);
518 rav
= dr
->aver1
[i
]/nsteps
;
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
]);
533 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n",max_dr
,hi
);
536 printf("Highest level in the matrix will be %g\n",hi
);
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
);
543 int gmx_disre(int argc
,char *argv
[])
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",
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."
565 static int nlevels
= 20;
566 static real max_dr
= 0;
567 static gmx_bool bThird
= TRUE
;
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" }
579 FILE *out
=NULL
,*aver
=NULL
,*numv
=NULL
,*maxxv
=NULL
,*xvg
=NULL
;
591 int ntopatoms
,natoms
,i
,j
,kkk
;
594 rvec
*x
,*f
,*xav
=NULL
;
598 atom_id
*index
=NULL
,*ind_fit
=NULL
;
600 t_cluster_ndx
*clust
=NULL
;
601 t_dr_result dr
,*dr_clust
=NULL
;
603 real
*vvindex
=NULL
,*w_rls
=NULL
;
609 gmx_rmpbc_t gpbc
=NULL
;
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
}
625 #define NFILE asize(fnm)
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
);
632 gmx_log_open(ftp2fn(efLOG
,NFILE
,fnm
),cr
,FALSE
,0,&fplog
);
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
);
643 snew(ind_fit
,ntopatoms
);
644 snew(w_rls
,ntopatoms
);
645 for(kkk
=0; (kkk
<ntopatoms
); kkk
++) {
651 *atoms
= gmx_mtop_global_atoms(&mtop
);
653 if (atoms
->pdbinfo
== NULL
) {
654 snew(atoms
->pdbinfo
,atoms
->nr
);
658 top
= gmx_mtop_generate_local_top(&mtop
,&ir
);
662 if (ir
.ePBC
!= epbcNONE
) {
663 if (ir
.bPeriodicMols
)
666 g
= mk_graph(fplog
,&top
->idef
,0,mtop
.natoms
,FALSE
,FALSE
);
669 if (ftp2bSet(efNDX
,NFILE
,fnm
)) {
670 rd_index(ftp2fn(efNDX
,NFILE
,fnm
),1,&isize
,&index
,&grpname
);
671 xvg
=xvgropen(opt2fn("-dr",NFILE
,fnm
),"Inidividual Restraints","Time (ps)",
675 for(i
=0; (i
<isize
); i
++) {
678 sprintf(leg
[i
],"index %d",index
[i
]);
680 xvgr_legend(xvg
,isize
,(const char**)leg
,oenv
);
686 init_disres(fplog
,&mtop
,&ir
,NULL
,FALSE
,&fcd
,NULL
);
688 natoms
=read_first_x(oenv
,&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
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
);
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
);
709 mdatoms
= init_mdatoms(fplog
,&mtop
,ir
.efep
!=efepNO
);
710 atoms2md(&mtop
,&ir
,0,NULL
,0,mtop
.natoms
,mdatoms
);
711 update_mdatoms(mdatoms
,ir
.init_lambda
);
713 fprintf(fplog
,"Made forcerec\n");
714 init_forcerec(fplog
,oenv
,fr
,NULL
,&ir
,&mtop
,cr
,box
,FALSE
,NULL
,NULL
,NULL
,
717 if (ir
.ePBC
!= epbcNONE
)
718 gpbc
= gmx_rmpbc_init(&top
->idef
,ir
.ePBC
,natoms
,box
);
722 if (ir
.ePBC
!= epbcNONE
) {
723 if (ir
.bPeriodicMols
)
724 set_pbc(&pbc
,ir
.ePBC
,box
);
726 gmx_rmpbc(gpbc
,natoms
,box
,x
);
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
);
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
);
743 reset_x(atoms
->nr
,ind_fit
,atoms
->nr
,NULL
,x
,w_rls
);
744 do_fit(atoms
->nr
,w_rls
,x
,x
);
746 /* Store the first frame of the trajectory as 'characteristic'
747 * for colouring with violations.
749 for(kkk
=0; (kkk
<atoms
->nr
); kkk
++)
750 copy_rvec(x
[kkk
],xav
[kkk
]);
755 fprintf(xvg
,"%10g",t
);
756 for(i
=0; (i
<isize
); i
++)
757 fprintf(xvg
," %10g",vvindex
[i
]);
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
);
766 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
768 if (ir
.ePBC
!= epbcNONE
)
769 gmx_rmpbc_done(gpbc
);
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
);
777 dump_stats(fplog
,j
,fcd
.disres
.nres
,&(top
->idef
.il
[F_DISRES
]),
778 top
->idef
.iparams
,&dr
,isize
,index
,
779 bPDB
? atoms
: NULL
);
781 write_sto_conf(opt2fn("-q",NFILE
,fnm
),
782 "Coloured by average violation in Angstrom",
783 atoms
,xav
,NULL
,ir
.ePBC
,box
);
785 dump_disre_matrix(opt2fn_null("-x",NFILE
,fnm
),&dr
,fcd
.disres
.nres
,
786 j
,&top
->idef
,&mtop
,max_dr
,nlevels
,bThird
);
793 do_view(oenv
,opt2fn("-dr",NFILE
,fnm
),"-nxy");
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");
802 if (gmx_parallel_env_initialized())
805 gmx_log_close(fplog
);