?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / kernel / tpbcmp.c
blob4b975d5ec33b79841fb3ec0855ea0219b408e7bd
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 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
34 */
35 /* This file is completely threadsafe - keep it that way! */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include <stdio.h>
42 #include <string.h>
43 #include "main.h"
44 #include "macros.h"
45 #include "smalloc.h"
46 #include "futil.h"
47 #include "statutil.h"
48 #include "sysstuff.h"
49 #include "txtdump.h"
50 #include "gmx_fatal.h"
51 #include "names.h"
52 #include "tpxio.h"
53 #include "enxio.h"
54 #include "mtop_util.h"
55 #include "string2.h"
56
57 static void cmp_int(FILE *fp,const char *s,int index,int i1,int i2)
58 {
59 if (i1 != i2) {
60 if (index != -1)
61 fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
62 else
63 fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
64 }
65 }
66
67 static void cmp_gmx_large_int(FILE *fp,const char *s,gmx_large_int_t i1,gmx_large_int_t i2)
68 {
69 if (i1 != i2) {
70 fprintf(fp,"%s (",s);
71 fprintf(fp,gmx_large_int_pfmt,i1);
72 fprintf(fp," - ");
73 fprintf(fp,gmx_large_int_pfmt,i2);
74 fprintf(fp,")\n");
75 }
76 }
77
78 static void cmp_us(FILE *fp,const char *s,int index,unsigned short i1,unsigned short i2)
79 {
80 if (i1 != i2) {
81 if (index != -1)
82 fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
83 else
84 fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
85 }
86 }
87
88 static void cmp_uc(FILE *fp,const char *s,int index,unsigned char i1,unsigned char i2)
89 {
90 if (i1 != i2) {
91 if (index != -1)
92 fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
93 else
94 fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
95 }
96 }
97
98 static gmx_bool cmp_bool(FILE *fp, const char *s, int index, gmx_bool b1, gmx_bool b2)
99 {
100 if (b1) {
101 b1 = 1;
102 } else {
103 b1 = 0;
104 }
105 if (b2) {
106 b2 = 1;
107 } else {
108 b2 = 0;
109 }
110 if (b1 != b2) {
111 if (index != -1)
112 fprintf(fp,"%s[%d] (%s - %s)\n",s,index,
113 bool_names[b1],bool_names[b2]);
114 else
115 fprintf(fp,"%s (%s - %s)\n",s,
116 bool_names[b1],bool_names[b2]);
117 }
118 return b1 && b2;
119 }
120
121 static void cmp_str(FILE *fp, const char *s, int index,
122 const char *s1, const char *s2)
123 {
124 if (strcmp(s1,s2) != 0) {
125 if (index != -1)
126 fprintf(fp,"%s[%d] (%s - %s)\n",s,index,s1,s2);
127 else
128 fprintf(fp,"%s (%s - %s)\n",s,s1,s2);
129 }
130 }
131
132 static gmx_bool equal_real(real i1,real i2,real ftol,real abstol)
133 {
134 return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
135 }
136
137 static gmx_bool equal_float(float i1,float i2,float ftol,float abstol)
138 {
139 return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
140 }
141
142 static gmx_bool equal_double(double i1,double i2,real ftol,real abstol)
143 {
144 return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
145 }
146
147 static void
148 cmp_real(FILE *fp,const char *s,int index,real i1,real i2,real ftol,real abstol)
149 {
150 if (!equal_real(i1,i2,ftol,abstol)) {
151 if (index != -1)
152 fprintf(fp,"%s[%2d] (%e - %e)\n",s,index,i1,i2);
153 else
154 fprintf(fp,"%s (%e - %e)\n",s,i1,i2);
155 }
156 }
157
158 static void
159 cmp_float(FILE *fp,const char *s,int index,float i1,float i2,float ftol,float abstol)
160 {
161 if (!equal_float(i1,i2,ftol,abstol)) {
162 if (index != -1)
163 fprintf(fp,"%s[%2d] (%e - %e)\n",s,index,i1,i2);
164 else
165 fprintf(fp,"%s (%e - %e)\n",s,i1,i2);
166 }
167 }
168
169
170
171 static void
172 cmp_double(FILE *fp,const char *s,int index,double i1,double i2,double ftol,double abstol)
173 {
174 if (!equal_double(i1,i2,ftol,abstol)) {
175 if (index != -1)
176 fprintf(fp,"%s[%2d] (%16.9e - %16.9e)\n",s,index,i1,i2);
177 else
178 fprintf(fp,"%s (%16.9e - %16.9e)\n",s,i1,i2);
179 }
180 }
181
182 static void cmp_rvec(FILE *fp,const char *s,int index,rvec i1,rvec i2,real ftol,real abstol)
183 {
184 if(!equal_real(i1[XX],i2[XX],ftol,abstol) ||
185 !equal_real(i1[YY],i2[YY],ftol,abstol) ||
186 !equal_real(i1[ZZ],i2[ZZ],ftol,abstol))
187 {
188 if (index != -1)
189 fprintf(fp,"%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
190 s,index,i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
191 else
192 fprintf(fp,"%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
193 s,i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
194 }
195 }
196
197 static void cmp_ivec(FILE *fp,const char *s,int index,ivec i1,ivec i2)
198 {
199 if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ])) {
200 if (index != -1)
201 fprintf(fp,"%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n",s,index,
202 i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
203 else
204 fprintf(fp,"%s (%8d,%8d,%8d - %8d,%8d,%8d)\n",s,
205 i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
206 }
207 }
208
209 static void cmp_ilist(FILE *fp,int ftype,t_ilist *il1,t_ilist *il2)
210 {
211 int i;
212 char buf[256];
213
214 fprintf(fp,"comparing ilist %s\n",interaction_function[ftype].name);
215 sprintf(buf,"%s->nr",interaction_function[ftype].name);
216 cmp_int(fp,buf,-1,il1->nr,il2->nr);
217 sprintf(buf,"%s->iatoms",interaction_function[ftype].name);
218 if (((il1->nr > 0) && (!il1->iatoms)) ||
219 ((il2->nr > 0) && (!il2->iatoms)) ||
220 ((il1->nr != il2->nr)))
221 fprintf(fp,"Comparing radically different topologies - %s is different\n",
222 buf);
223 else
224 for(i=0; (i<il1->nr); i++)
225 cmp_int(fp,buf,i,il1->iatoms[i],il2->iatoms[i]);
226 }
227
228 void cmp_iparm(FILE *fp,const char *s,t_functype ft,
229 t_iparams ip1,t_iparams ip2,real ftol,real abstol)
230 {
231 int i;
232 gmx_bool bDiff;
233
234 bDiff=FALSE;
235 for(i=0; i<MAXFORCEPARAM && !bDiff; i++)
236 bDiff = !equal_real(ip1.generic.buf[i],ip2.generic.buf[i],ftol,abstol);
237 if (bDiff) {
238 fprintf(fp,"%s1: ",s);
239 pr_iparams(fp,ft,&ip1);
240 fprintf(fp,"%s2: ",s);
241 pr_iparams(fp,ft,&ip2);
242 }
243 }
244
245 void cmp_iparm_AB(FILE *fp,const char *s,t_functype ft,t_iparams ip1,real ftol,real abstol)
246 {
247 int nrfpA,nrfpB,p0,i;
248 gmx_bool bDiff;
249
250 /* Normally the first parameter is perturbable */
251 p0 = 0;
252 nrfpA = interaction_function[ft].nrfpA;
253 nrfpB = interaction_function[ft].nrfpB;
254 if (ft == F_PDIHS) {
255 nrfpB = 2;
256 } else if (interaction_function[ft].flags & IF_TABULATED) {
257 /* For tabulated interactions only the second parameter is perturbable */
258 p0 = 1;
259 nrfpB = 1;
260 }
261 bDiff=FALSE;
262 for(i=0; i<nrfpB && !bDiff; i++) {
263 bDiff = !equal_real(ip1.generic.buf[p0+i],ip1.generic.buf[nrfpA+i],ftol,abstol);
264 }
265 if (bDiff) {
266 fprintf(fp,"%s: ",s);
267 pr_iparams(fp,ft,&ip1);
268 }
269 }
270
271 static void cmp_idef(FILE *fp,t_idef *id1,t_idef *id2,real ftol,real abstol)
272 {
273 int i;
274 char buf1[64],buf2[64];
275
276 fprintf(fp,"comparing idef\n");
277 if (id2) {
278 cmp_int(fp,"idef->ntypes",-1,id1->ntypes,id2->ntypes);
279 cmp_int(fp,"idef->atnr", -1,id1->atnr,id2->atnr);
280 for(i=0; (i<id1->ntypes); i++) {
281 sprintf(buf1,"idef->functype[%d]",i);
282 sprintf(buf2,"idef->iparam[%d]",i);
283 cmp_int(fp,buf1,i,(int)id1->functype[i],(int)id2->functype[i]);
284 cmp_iparm(fp,buf2,id1->functype[i],
285 id1->iparams[i],id2->iparams[i],ftol,abstol);
286 }
287 cmp_real(fp,"fudgeQQ",-1,id1->fudgeQQ,id2->fudgeQQ,ftol,abstol);
288 for(i=0; (i<F_NRE); i++)
289 cmp_ilist(fp,i,&(id1->il[i]),&(id2->il[i]));
290 } else {
291 for(i=0; (i<id1->ntypes); i++)
292 cmp_iparm_AB(fp,"idef->iparam",id1->functype[i],id1->iparams[i],ftol,abstol);
293 }
294 }
295
296 static void cmp_block(FILE *fp,t_block *b1,t_block *b2,const char *s)
297 {
298 int i,j,k;
299 char buf[32];
300
301 fprintf(fp,"comparing block %s\n",s);
302 sprintf(buf,"%s.nr",s);
303 cmp_int(fp,buf,-1,b1->nr,b2->nr);
304 }
305
306 static void cmp_blocka(FILE *fp,t_blocka *b1,t_blocka *b2,const char *s)
307 {
308 int i,j,k;
309 char buf[32];
310
311 fprintf(fp,"comparing blocka %s\n",s);
312 sprintf(buf,"%s.nr",s);
313 cmp_int(fp,buf,-1,b1->nr,b2->nr);
314 sprintf(buf,"%s.nra",s);
315 cmp_int(fp,buf,-1,b1->nra,b2->nra);
316 }
317
318 static void cmp_atom(FILE *fp,int index,t_atom *a1,t_atom *a2,real ftol,real abstol)
319 {
320 int i;
321 char buf[256];
322
323 if (a2) {
324 cmp_us(fp,"atom.type",index,a1->type,a2->type);
325 cmp_us(fp,"atom.ptype",index,a1->ptype,a2->ptype);
326 cmp_int(fp,"atom.resind",index,a1->resind,a2->resind);
327 cmp_int(fp,"atom.atomnumber",index,a1->atomnumber,a2->atomnumber);
328 cmp_real(fp,"atom.m",index,a1->m,a2->m,ftol,abstol);
329 cmp_real(fp,"atom.q",index,a1->q,a2->q,ftol,abstol);
330 cmp_us(fp,"atom.typeB",index,a1->typeB,a2->typeB);
331 cmp_real(fp,"atom.mB",index,a1->mB,a2->mB,ftol,abstol);
332 cmp_real(fp,"atom.qB",index,a1->qB,a2->qB,ftol,abstol);
333 } else {
334 cmp_us(fp,"atom.type",index,a1->type,a1->typeB);
335 cmp_real(fp,"atom.m",index,a1->m,a1->mB,ftol,abstol);
336 cmp_real(fp,"atom.q",index,a1->q,a1->qB,ftol,abstol);
337 }
338 }
339
340 static void cmp_atoms(FILE *fp,t_atoms *a1,t_atoms *a2,real ftol, real abstol)
341 {
342 int i;
343
344 fprintf(fp,"comparing atoms\n");
345
346 if (a2) {
347 cmp_int(fp,"atoms->nr",-1,a1->nr,a2->nr);
348 for(i=0; (i<a1->nr); i++)
349 cmp_atom(fp,i,&(a1->atom[i]),&(a2->atom[i]),ftol,abstol);
350 } else {
351 for(i=0; (i<a1->nr); i++)
352 cmp_atom(fp,i,&(a1->atom[i]),NULL,ftol,abstol);
353 }
354 }
355
356 static void cmp_top(FILE *fp,t_topology *t1,t_topology *t2,real ftol, real abstol)
357 {
358 int i;
359
360 fprintf(fp,"comparing top\n");
361 if (t2) {
362 cmp_idef(fp,&(t1->idef),&(t2->idef),ftol,abstol);
363 cmp_atoms(fp,&(t1->atoms),&(t2->atoms),ftol,abstol);
364 cmp_block(fp,&t1->cgs,&t2->cgs,"cgs");
365 cmp_block(fp,&t1->mols,&t2->mols,"mols");
366 cmp_blocka(fp,&t1->excls,&t2->excls,"excls");
367 } else {
368 cmp_idef(fp,&(t1->idef),NULL,ftol,abstol);
369 cmp_atoms(fp,&(t1->atoms),NULL,ftol,abstol);
370 }
371 }
372
373 static void cmp_groups(FILE *fp,gmx_groups_t *g0,gmx_groups_t *g1,
374 int natoms0,int natoms1)
375 {
376 int i,j,ndiff;
377 char buf[32];
378
379 fprintf(fp,"comparing groups\n");
380
381 for(i=0; i<egcNR; i++) {
382 sprintf(buf,"grps[%d].nr",i);
383 cmp_int(fp,buf,-1,g0->grps[i].nr,g1->grps[i].nr);
384 if (g0->grps[i].nr == g1->grps[i].nr) {
385 for(j=0; j<g0->grps[i].nr; j++) {
386 sprintf(buf,"grps[%d].name[%d]",i,j);
387 cmp_str(fp,buf,-1,
388 *g0->grpname[g0->grps[i].nm_ind[j]],
389 *g1->grpname[g1->grps[i].nm_ind[j]]);
390 }
391 }
392 cmp_int(fp,"ngrpnr",i,g0->ngrpnr[i],g1->ngrpnr[i]);
393 if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
394 (g0->grpnr[i] != NULL || g1->grpnr[i] != NULL)) {
395 for(j=0; j<natoms0; j++) {
396 cmp_int(fp,gtypes[i],j,ggrpnr(g0,i,j),ggrpnr(g1,i,j));
397 }
398 }
399 }
400 /* We have compared the names in the groups lists,
401 * so we can skip the grpname list comparison.
402 */
403 }
404
405 static void cmp_rvecs(FILE *fp,const char *title,int n,rvec x1[],rvec x2[],
406 gmx_bool bRMSD,real ftol,real abstol)
407 {
408 int i,m;
409 double d,ssd;
410
411 if (bRMSD) {
412 ssd = 0;
413 for(i=0; (i<n); i++) {
414 for(m=0; m<DIM; m++) {
415 d = x1[i][m] - x2[i][m];
416 ssd += d*d;
417 }
418 }
419 fprintf(fp,"%s RMSD %g\n",title,sqrt(ssd/n));
420 } else {
421 for(i=0; (i<n); i++) {
422 cmp_rvec(fp,title,i,x1[i],x2[i],ftol,abstol);
423 }
424 }
425 }
426
427 static void cmp_grpopts(FILE *fp,t_grpopts *opt1,t_grpopts *opt2,real ftol, real abstol)
428 {
429 int i,j;
430 char buf1[256],buf2[256];
431
432 cmp_int(fp,"inputrec->grpopts.ngtc",-1, opt1->ngtc,opt2->ngtc);
433 cmp_int(fp,"inputrec->grpopts.ngacc",-1, opt1->ngacc,opt2->ngacc);
434 cmp_int(fp,"inputrec->grpopts.ngfrz",-1, opt1->ngfrz,opt2->ngfrz);
435 cmp_int(fp,"inputrec->grpopts.ngener",-1,opt1->ngener,opt2->ngener);
436 for(i=0; (i<min(opt1->ngtc,opt2->ngtc)); i++) {
437 cmp_real(fp,"inputrec->grpopts.nrdf",i,opt1->nrdf[i],opt2->nrdf[i],ftol,abstol);
438 cmp_real(fp,"inputrec->grpopts.ref_t",i,opt1->ref_t[i],opt2->ref_t[i],ftol,abstol);
439 cmp_real(fp,"inputrec->grpopts.tau_t",i,opt1->tau_t[i],opt2->tau_t[i],ftol,abstol);
440 cmp_int(fp,"inputrec->grpopts.annealing",i,opt1->annealing[i],opt2->annealing[i]);
441 cmp_int(fp,"inputrec->grpopts.anneal_npoints",i,
442 opt1->anneal_npoints[i],opt2->anneal_npoints[i]);
443 if(opt1->anneal_npoints[i]==opt2->anneal_npoints[i]) {
444 sprintf(buf1,"inputrec->grpopts.anneal_time[%d]",i);
445 sprintf(buf2,"inputrec->grpopts.anneal_temp[%d]",i);
446 for(j=0;j<opt1->anneal_npoints[i];j++) {
447 cmp_real(fp,buf1,j,opt1->anneal_time[i][j],opt2->anneal_time[i][j],ftol,abstol);
448 cmp_real(fp,buf2,j,opt1->anneal_temp[i][j],opt2->anneal_temp[i][j],ftol,abstol);
449 }
450 }
451 }
452 if (opt1->ngener == opt2->ngener) {
453 for(i=0; i<opt1->ngener; i++)
454 for(j=i; j<opt1->ngener; j++) {
455 sprintf(buf1,"inputrec->grpopts.egp_flags[%d]",i);
456 cmp_int(fp,buf1,j,
457 opt1->egp_flags[opt1->ngener*i+j],
458 opt2->egp_flags[opt1->ngener*i+j]);
459 }
460 }
461 for(i=0; (i<min(opt1->ngacc,opt2->ngacc)); i++)
462 cmp_rvec(fp,"inputrec->grpopts.acc",i,opt1->acc[i],opt2->acc[i],ftol,abstol);
463 for(i=0; (i<min(opt1->ngfrz,opt2->ngfrz)); i++)
464 cmp_ivec(fp,"inputrec->grpopts.nFreeze",i,opt1->nFreeze[i],opt2->nFreeze[i]);
465 }
466
467 static void cmp_cosines(FILE *fp,const char *s,t_cosines c1[DIM],t_cosines c2[DIM],real ftol, real abstol)
468 {
469 int i,m;
470 char buf[256];
471
472 for(m=0; (m<DIM); m++) {
473 sprintf(buf,"inputrec->%s[%d]",s,m);
474 cmp_int(fp,buf,0,c1->n,c2->n);
475 for(i=0; (i<min(c1->n,c2->n)); i++) {
476 cmp_real(fp,buf,i,c1->a[i],c2->a[i],ftol,abstol);
477 cmp_real(fp,buf,i,c1->phi[i],c2->phi[i],ftol,abstol);
478 }
479 }
480 }
481
482 static void cmp_pull(FILE *fp,t_pull *pull1,t_pull *pull2,real ftol, real abstol)
483 {
484 fprintf(fp,"WARNING: Both files use COM pulling, but comparing of the pull struct is not implemented (yet). The pull parameters could be the same or different.\n");
485 }
486
487 static void cmp_inputrec(FILE *fp,t_inputrec *ir1,t_inputrec *ir2,real ftol, real abstol)
488 {
489 int i;
490
491 fprintf(fp,"comparing inputrec\n");
492
493 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
494 * of warnings. Maybe it will change in future versions, but for the
495 * moment I've spelled them out instead. /EL 000820
496 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
497 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
498 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
499 */
500 cmp_int(fp,"inputrec->eI",-1,ir1->eI,ir2->eI);
501 cmp_gmx_large_int(fp,"inputrec->nsteps",ir1->nsteps,ir2->nsteps);
502 cmp_gmx_large_int(fp,"inputrec->init_step",ir1->init_step,ir2->init_step);
503 cmp_int(fp,"inputrec->simulation_part",-1,ir1->simulation_part,ir2->simulation_part);
504 cmp_int(fp,"inputrec->ePBC",-1,ir1->ePBC,ir2->ePBC);
505 cmp_int(fp,"inputrec->bPeriodicMols",-1,ir1->bPeriodicMols,ir2->bPeriodicMols);
506 cmp_int(fp,"inputrec->ns_type",-1,ir1->ns_type,ir2->ns_type);
507 cmp_int(fp,"inputrec->nstlist",-1,ir1->nstlist,ir2->nstlist);
508 cmp_int(fp,"inputrec->ndelta",-1,ir1->ndelta,ir2->ndelta);
509 cmp_int(fp,"inputrec->nstcomm",-1,ir1->nstcomm,ir2->nstcomm);
510 cmp_int(fp,"inputrec->comm_mode",-1,ir1->comm_mode,ir2->comm_mode);
511 cmp_int(fp,"inputrec->nstcheckpoint",-1,ir1->nstcheckpoint,ir2->nstcheckpoint);
512 cmp_int(fp,"inputrec->nstlog",-1,ir1->nstlog,ir2->nstlog);
513 cmp_int(fp,"inputrec->nstxout",-1,ir1->nstxout,ir2->nstxout);
514 cmp_int(fp,"inputrec->nstvout",-1,ir1->nstvout,ir2->nstvout);
515 cmp_int(fp,"inputrec->nstfout",-1,ir1->nstfout,ir2->nstfout);
516 cmp_int(fp,"inputrec->nstcalcenergy",-1,ir1->nstcalcenergy,ir2->nstcalcenergy);
517 cmp_int(fp,"inputrec->nstenergy",-1,ir1->nstenergy,ir2->nstenergy);
518 cmp_int(fp,"inputrec->nstxtcout",-1,ir1->nstxtcout,ir2->nstxtcout);
519 cmp_double(fp,"inputrec->init_t",-1,ir1->init_t,ir2->init_t,ftol,abstol);
520 cmp_double(fp,"inputrec->delta_t",-1,ir1->delta_t,ir2->delta_t,ftol,abstol);
521 cmp_real(fp,"inputrec->xtcprec",-1,ir1->xtcprec,ir2->xtcprec,ftol,abstol);
522 cmp_int(fp,"inputrec->nkx",-1,ir1->nkx,ir2->nkx);
523 cmp_int(fp,"inputrec->nky",-1,ir1->nky,ir2->nky);
524 cmp_int(fp,"inputrec->nkz",-1,ir1->nkz,ir2->nkz);
525 cmp_int(fp,"inputrec->pme_order",-1,ir1->pme_order,ir2->pme_order);
526 cmp_real(fp,"inputrec->ewald_rtol",-1,ir1->ewald_rtol,ir2->ewald_rtol,ftol,abstol);
527 cmp_int(fp,"inputrec->ewald_geometry",-1,ir1->ewald_geometry,ir2->ewald_geometry);
528 cmp_real(fp,"inputrec->epsilon_surface",-1,ir1->epsilon_surface,ir2->epsilon_surface,ftol,abstol);
529 cmp_int(fp,"inputrec->bOptFFT",-1,ir1->bOptFFT,ir2->bOptFFT);
530 cmp_int(fp,"inputrec->bContinuation",-1,ir1->bContinuation,ir2->bContinuation);
531 cmp_int(fp,"inputrec->bShakeSOR",-1,ir1->bShakeSOR,ir2->bShakeSOR);
532 cmp_int(fp,"inputrec->etc",-1,ir1->etc,ir2->etc);
533 cmp_int(fp,"inputrec->epc",-1,ir1->epc,ir2->epc);
534 cmp_int(fp,"inputrec->epct",-1,ir1->epct,ir2->epct);
535 cmp_real(fp,"inputrec->tau_p",-1,ir1->tau_p,ir2->tau_p,ftol,abstol);
536 cmp_rvec(fp,"inputrec->ref_p(x)",-1,ir1->ref_p[XX],ir2->ref_p[XX],ftol,abstol);
537 cmp_rvec(fp,"inputrec->ref_p(y)",-1,ir1->ref_p[YY],ir2->ref_p[YY],ftol,abstol);
538 cmp_rvec(fp,"inputrec->ref_p(z)",-1,ir1->ref_p[ZZ],ir2->ref_p[ZZ],ftol,abstol);
539 cmp_rvec(fp,"inputrec->compress(x)",-1,ir1->compress[XX],ir2->compress[XX],ftol,abstol);
540 cmp_rvec(fp,"inputrec->compress(y)",-1,ir1->compress[YY],ir2->compress[YY],ftol,abstol);
541 cmp_rvec(fp,"inputrec->compress(z)",-1,ir1->compress[ZZ],ir2->compress[ZZ],ftol,abstol);
542 cmp_int(fp,"refcoord_scaling",-1,ir1->refcoord_scaling,ir2->refcoord_scaling);
543 cmp_rvec(fp,"inputrec->posres_com",-1,ir1->posres_com,ir2->posres_com,ftol,abstol);
544 cmp_rvec(fp,"inputrec->posres_comB",-1,ir1->posres_comB,ir2->posres_comB,ftol,abstol);
545 cmp_int(fp,"inputrec->andersen_seed",-1,ir1->andersen_seed,ir2->andersen_seed);
546 cmp_real(fp,"inputrec->rlist",-1,ir1->rlist,ir2->rlist,ftol,abstol);
547 cmp_real(fp,"inputrec->rlistlong",-1,ir1->rlistlong,ir2->rlistlong,ftol,abstol);
548 cmp_real(fp,"inputrec->rtpi",-1,ir1->rtpi,ir2->rtpi,ftol,abstol);
549 cmp_int(fp,"inputrec->coulombtype",-1,ir1->coulombtype,ir2->coulombtype);
550 cmp_real(fp,"inputrec->rcoulomb_switch",-1,ir1->rcoulomb_switch,ir2->rcoulomb_switch,ftol,abstol);
551 cmp_real(fp,"inputrec->rcoulomb",-1,ir1->rcoulomb,ir2->rcoulomb,ftol,abstol);
552 cmp_int(fp,"inputrec->vdwtype",-1,ir1->vdwtype,ir2->vdwtype);
553 cmp_real(fp,"inputrec->rvdw_switch",-1,ir1->rvdw_switch,ir2->rvdw_switch,ftol,abstol);
554 cmp_real(fp,"inputrec->rvdw",-1,ir1->rvdw,ir2->rvdw,ftol,abstol);
555 cmp_real(fp,"inputrec->epsilon_r",-1,ir1->epsilon_r,ir2->epsilon_r,ftol,abstol);
556 cmp_real(fp,"inputrec->epsilon_rf",-1,ir1->epsilon_rf,ir2->epsilon_rf,ftol,abstol);
557 cmp_real(fp,"inputrec->tabext",-1,ir1->tabext,ir2->tabext,ftol,abstol);
558 cmp_int(fp,"inputrec->implicit_solvent",-1,ir1->implicit_solvent,ir2->implicit_solvent);
559 cmp_int(fp,"inputrec->gb_algorithm",-1,ir1->gb_algorithm,ir2->gb_algorithm);
560 cmp_int(fp,"inputrec->nstgbradii",-1,ir1->nstgbradii,ir2->nstgbradii);
561 cmp_real(fp,"inputrec->rgbradii",-1,ir1->rgbradii,ir2->rgbradii,ftol,abstol);
562 cmp_real(fp,"inputrec->gb_saltconc",-1,ir1->gb_saltconc,ir2->gb_saltconc,ftol,abstol);
563 cmp_real(fp,"inputrec->gb_epsilon_solvent",-1,ir1->gb_epsilon_solvent,ir2->gb_epsilon_solvent,ftol,abstol);
564 cmp_real(fp,"inputrec->gb_obc_alpha",-1,ir1->gb_obc_alpha,ir2->gb_obc_alpha,ftol,abstol);
565 cmp_real(fp,"inputrec->gb_obc_beta",-1,ir1->gb_obc_beta,ir2->gb_obc_beta,ftol,abstol);
566 cmp_real(fp,"inputrec->gb_obc_gamma",-1,ir1->gb_obc_gamma,ir2->gb_obc_gamma,ftol,abstol);
567 cmp_real(fp,"inputrec->gb_dielectric_offset",-1,ir1->gb_dielectric_offset,ir2->gb_dielectric_offset,ftol,abstol);
568 cmp_int(fp,"inputrec->sa_algorithm",-1,ir1->sa_algorithm,ir2->sa_algorithm);
569 cmp_real(fp,"inputrec->sa_surface_tension",-1,ir1->sa_surface_tension,ir2->sa_surface_tension,ftol,abstol);
570
571 cmp_int(fp,"inputrec->eDispCorr",-1,ir1->eDispCorr,ir2->eDispCorr);
572 cmp_real(fp,"inputrec->shake_tol",-1,ir1->shake_tol,ir2->shake_tol,ftol,abstol);
573 cmp_int(fp,"inputrec->efep",-1,ir1->efep,ir2->efep);
574 cmp_double(fp,"inputrec->init_lambda",-1,ir1->init_lambda,ir2->init_lambda,ftol,abstol);
575 cmp_double(fp,"inputrec->delta_lambda",-1,ir1->delta_lambda,ir2->delta_lambda,ftol,abstol);
576 cmp_int(fp,"inputrec->n_foreign_lambda",-1,ir1->n_flambda,ir2->n_flambda);
577 for(i=0; i<min(ir1->n_flambda,ir2->n_flambda); i++) {
578 cmp_double(fp,"inputrec->foreign_lambda",-1,ir1->flambda[i],ir2->flambda[i],ftol,abstol);
579 }
580 cmp_real(fp,"inputrec->sc_alpha",-1,ir1->sc_alpha,ir2->sc_alpha,ftol,abstol);
581 cmp_int(fp,"inputrec->sc_power",-1,ir1->sc_power,ir2->sc_power);
582 cmp_real(fp,"inputrec->sc_sigma",-1,ir1->sc_sigma,ir2->sc_sigma,ftol,abstol);
583 cmp_real(fp,"inputrec->sc_sigma_min",-1,ir1->sc_sigma_min,ir2->sc_sigma_min,ftol,abstol);
584 cmp_int(fp,"inputrec->nstdhdl",-1,ir1->nstdhdl,ir2->nstdhdl);
585 cmp_int(fp,"inputrec->separate_dhdl_file",-1,ir1->separate_dhdl_file,ir2->nstdhdl);
586 cmp_int(fp,"inputrec->dhdl_derivatives",-1,ir1->dhdl_derivatives,ir2->dhdl_derivatives);
587 cmp_int(fp,"inputrec->dh_hist_size",-1,ir1->dh_hist_size,ir2->dh_hist_size);
588 cmp_double(fp,"inputrec->dh_hist_spacing",-1,ir1->dh_hist_spacing,ir2->dh_hist_spacing,ftol,abstol);
589
590 cmp_int(fp,"inputrec->nwall",-1,ir1->nwall,ir2->nwall);
591 cmp_int(fp,"inputrec->wall_type",-1,ir1->wall_type,ir2->wall_type);
592 cmp_int(fp,"inputrec->wall_atomtype[0]",-1,ir1->wall_atomtype[0],ir2->wall_atomtype[0]);
593 cmp_int(fp,"inputrec->wall_atomtype[1]",-1,ir1->wall_atomtype[1],ir2->wall_atomtype[1]);
594 cmp_real(fp,"inputrec->wall_density[0]",-1,ir1->wall_density[0],ir2->wall_density[0],ftol,abstol);
595 cmp_real(fp,"inputrec->wall_density[1]",-1,ir1->wall_density[1],ir2->wall_density[1],ftol,abstol);
596 cmp_real(fp,"inputrec->wall_ewald_zfac",-1,ir1->wall_ewald_zfac,ir2->wall_ewald_zfac,ftol,abstol);
597
598 cmp_int(fp,"inputrec->ePull",-1,ir1->ePull,ir2->ePull);
599 if (ir1->ePull == ir2->ePull && ir1->ePull != epullNO)
600 cmp_pull(fp,ir1->pull,ir2->pull,ftol,abstol);
601
602 cmp_int(fp,"inputrec->eDisre",-1,ir1->eDisre,ir2->eDisre);
603 cmp_real(fp,"inputrec->dr_fc",-1,ir1->dr_fc,ir2->dr_fc,ftol,abstol);
604 cmp_int(fp,"inputrec->eDisreWeighting",-1,ir1->eDisreWeighting,ir2->eDisreWeighting);
605 cmp_int(fp,"inputrec->bDisreMixed",-1,ir1->bDisreMixed,ir2->bDisreMixed);
606 cmp_int(fp,"inputrec->nstdisreout",-1,ir1->nstdisreout,ir2->nstdisreout);
607 cmp_real(fp,"inputrec->dr_tau",-1,ir1->dr_tau,ir2->dr_tau,ftol,abstol);
608 cmp_real(fp,"inputrec->orires_fc",-1,ir1->orires_fc,ir2->orires_fc,ftol,abstol);
609 cmp_real(fp,"inputrec->orires_tau",-1,ir1->orires_tau,ir2->orires_tau,ftol,abstol);
610 cmp_int(fp,"inputrec->nstorireout",-1,ir1->nstorireout,ir2->nstorireout);
611 cmp_real(fp,"inputrec->dihre_fc",-1,ir1->dihre_fc,ir2->dihre_fc,ftol,abstol);
612 cmp_real(fp,"inputrec->em_stepsize",-1,ir1->em_stepsize,ir2->em_stepsize,ftol,abstol);
613 cmp_real(fp,"inputrec->em_tol",-1,ir1->em_tol,ir2->em_tol,ftol,abstol);
614 cmp_int(fp,"inputrec->niter",-1,ir1->niter,ir2->niter);
615 cmp_real(fp,"inputrec->fc_stepsize",-1,ir1->fc_stepsize,ir2->fc_stepsize,ftol,abstol);
616 cmp_int(fp,"inputrec->nstcgsteep",-1,ir1->nstcgsteep,ir2->nstcgsteep);
617 cmp_int(fp,"inputrec->nbfgscorr",0,ir1->nbfgscorr,ir2->nbfgscorr);
618 cmp_int(fp,"inputrec->eConstrAlg",-1,ir1->eConstrAlg,ir2->eConstrAlg);
619 cmp_int(fp,"inputrec->nProjOrder",-1,ir1->nProjOrder,ir2->nProjOrder);
620 cmp_real(fp,"inputrec->LincsWarnAngle",-1,ir1->LincsWarnAngle,ir2->LincsWarnAngle,ftol,abstol);
621 cmp_int(fp,"inputrec->nLincsIter",-1,ir1->nLincsIter,ir2->nLincsIter);
622 cmp_real(fp,"inputrec->bd_fric",-1,ir1->bd_fric,ir2->bd_fric,ftol,abstol);
623 cmp_int(fp,"inputrec->ld_seed",-1,ir1->ld_seed,ir2->ld_seed);
624 cmp_real(fp,"inputrec->cos_accel",-1,ir1->cos_accel,ir2->cos_accel,ftol,abstol);
625 cmp_rvec(fp,"inputrec->deform(a)",-1,ir1->deform[XX],ir2->deform[XX],ftol,abstol);
626 cmp_rvec(fp,"inputrec->deform(b)",-1,ir1->deform[YY],ir2->deform[YY],ftol,abstol);
627 cmp_rvec(fp,"inputrec->deform(c)",-1,ir1->deform[ZZ],ir2->deform[ZZ],ftol,abstol);
628
629 cmp_int(fp,"ir->adress_type" ,-1,ir1->adress_type,ir2->adress_type);
630 cmp_int(fp,"ir->badress_tf_full_box" ,-1,ir1->badress_tf_full_box,ir2->badress_tf_full_box);
631 cmp_real(fp,"ir->adress_const_wf" ,-1,ir1->adress_const_wf,ir2->adress_const_wf,ftol,abstol);
632 cmp_real(fp,"ir->adress_ex_width" ,-1,ir1->adress_ex_width,ir2->adress_ex_width,ftol,abstol);
633 cmp_real(fp,"ir->adress_hy_width" ,-1,ir1->adress_hy_width,ir2->adress_hy_width,ftol,abstol);
634 cmp_int(fp,"ir->adress_icor" ,-1,ir1->adress_icor,ir2->adress_icor);
635 cmp_int(fp,"ir->adress_ivdw" ,-1,ir1->adress_ivdw,ir2->adress_ivdw);
636 cmp_int(fp,"ir->adress_site" ,-1,ir1->adress_site,ir2->adress_site);
637 cmp_rvec(fp,"ir->adress_refs" ,-1,ir1->adress_refs,ir2->adress_refs,ftol,abstol);
638 cmp_real(fp,"ir->adress_ex_forcecap", -1,ir1->adress_ex_forcecap,ir2->adress_ex_forcecap,ftol,abstol);
639
640 cmp_int(fp,"inputrec->userint1",-1,ir1->userint1,ir2->userint1);
641 cmp_int(fp,"inputrec->userint2",-1,ir1->userint2,ir2->userint2);
642 cmp_int(fp,"inputrec->userint3",-1,ir1->userint3,ir2->userint3);
643 cmp_int(fp,"inputrec->userint4",-1,ir1->userint4,ir2->userint4);
644 cmp_real(fp,"inputrec->userreal1",-1,ir1->userreal1,ir2->userreal1,ftol,abstol);
645 cmp_real(fp,"inputrec->userreal2",-1,ir1->userreal2,ir2->userreal2,ftol,abstol);
646 cmp_real(fp,"inputrec->userreal3",-1,ir1->userreal3,ir2->userreal3,ftol,abstol);
647 cmp_real(fp,"inputrec->userreal4",-1,ir1->userreal4,ir2->userreal4,ftol,abstol);
648 cmp_grpopts(fp,&(ir1->opts),&(ir2->opts),ftol,abstol);
649 cmp_cosines(fp,"ex",ir1->ex,ir2->ex,ftol,abstol);
650 cmp_cosines(fp,"et",ir1->et,ir2->et,ftol,abstol);
651 }
652
653 static void comp_pull_AB(FILE *fp,t_pull *pull,real ftol,real abstol)
654 {
655 int i;
656
657 for(i=0; i<pull->ngrp+1; i++) {
658 fprintf(fp,"comparing pull group %d\n",i);
659 cmp_real(fp,"pullgrp->k",-1,pull->grp[i].k,pull->grp[i].kB,ftol,abstol);
660 }
661 }
662
663 static void comp_state(t_state *st1, t_state *st2,
664 gmx_bool bRMSD,real ftol,real abstol)
665 {
666 int i,j,nc;
667
668 fprintf(stdout,"comparing flags\n");
669 cmp_int(stdout,"flags",-1,st1->flags,st2->flags);
670 fprintf(stdout,"comparing box\n");
671 cmp_rvecs(stdout,"box",DIM,st1->box,st2->box,FALSE,ftol,abstol);
672 fprintf(stdout,"comparing box_rel\n");
673 cmp_rvecs(stdout,"box_rel",DIM,st1->box_rel,st2->box_rel,FALSE,ftol,abstol);
674 fprintf(stdout,"comparing boxv\n");
675 cmp_rvecs(stdout,"boxv",DIM,st1->boxv,st2->boxv,FALSE,ftol,abstol);
676 if (st1->flags & (1<<estSVIR_PREV)) {
677 fprintf(stdout,"comparing shake vir_prev\n");
678 cmp_rvecs(stdout,"svir_prev",DIM,st1->svir_prev,st2->svir_prev,FALSE,ftol,abstol);
679 }
680 if (st1->flags & (1<<estFVIR_PREV)) {
681 fprintf(stdout,"comparing force vir_prev\n");
682 cmp_rvecs(stdout,"fvir_prev",DIM,st1->fvir_prev,st2->fvir_prev,FALSE,ftol,abstol);
683 }
684 if (st1->flags & (1<<estPRES_PREV)) {
685 fprintf(stdout,"comparing prev_pres\n");
686 cmp_rvecs(stdout,"pres_prev",DIM,st1->pres_prev,st2->pres_prev,FALSE,ftol,abstol);
687 }
688 cmp_int(stdout,"ngtc",-1,st1->ngtc,st2->ngtc);
689 cmp_int(stdout,"nhchainlength",-1,st1->nhchainlength,st2->nhchainlength);
690 if (st1->ngtc == st2->ngtc && st1->nhchainlength == st2->nhchainlength){
691 for(i=0; i<st1->ngtc; i++) {
692 nc = i*st1->nhchainlength;
693 for(j=0; j<nc; j++) {
694 cmp_real(stdout,"nosehoover_xi",
695 i,st1->nosehoover_xi[nc+j],st2->nosehoover_xi[nc+j],ftol,abstol);
696 }
697 }
698 }
699 cmp_int(stdout,"nnhpres",-1,st1->nnhpres,st2->nnhpres);
700 if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength) {
701 for(i=0; i<st1->nnhpres; i++) {
702 nc = i*st1->nhchainlength;
703 for(j=0; j<nc; j++) {
704 cmp_real(stdout,"nosehoover_xi",
705 i,st1->nhpres_xi[nc+j],st2->nhpres_xi[nc+j],ftol,abstol);
706 }
707 }
708 }
709
710 cmp_int(stdout,"natoms",-1,st1->natoms,st2->natoms);
711 if (st1->natoms == st2->natoms) {
712 if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX))) {
713 fprintf(stdout,"comparing x\n");
714 cmp_rvecs(stdout,"x",st1->natoms,st1->x,st2->x,bRMSD,ftol,abstol);
715 }
716 if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV))) {
717 fprintf(stdout,"comparing v\n");
718 cmp_rvecs(stdout,"v",st1->natoms,st1->v,st2->v,bRMSD,ftol,abstol);
719 }
720 }
721 }
722
723 void comp_tpx(const char *fn1,const char *fn2,
724 gmx_bool bRMSD,real ftol,real abstol)
725 {
726 const char *ff[2];
727 t_tpxheader sh[2];
728 t_inputrec ir[2];
729 t_state state[2];
730 gmx_mtop_t mtop[2];
731 t_topology top[2];
732 int i;
733
734 ff[0]=fn1;
735 ff[1]=fn2;
736 for(i=0; i<(fn2 ? 2 : 1); i++) {
737 read_tpx_state(ff[i],&(ir[i]),&state[i],NULL,&(mtop[i]));
738 }
739 if (fn2) {
740 cmp_inputrec(stdout,&ir[0],&ir[1],ftol,abstol);
741 /* Convert gmx_mtop_t to t_topology.
742 * We should implement direct mtop comparison,
743 * but it might be useful to keep t_topology comparison as an option.
744 */
745 top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
746 top[1] = gmx_mtop_t_to_t_topology(&mtop[1]);
747 cmp_top(stdout,&top[0],&top[1],ftol,abstol);
748 cmp_groups(stdout,&mtop[0].groups,&mtop[1].groups,
749 mtop[0].natoms,mtop[1].natoms);
750 comp_state(&state[0],&state[1],bRMSD,ftol,abstol);
751 } else {
752 if (ir[0].efep == efepNO) {
753 fprintf(stdout,"inputrec->efep = %s\n",efep_names[ir[0].efep]);
754 } else {
755 if (ir[0].ePull != epullNO) {
756 comp_pull_AB(stdout,ir->pull,ftol,abstol);
757 }
758 /* Convert gmx_mtop_t to t_topology.
759 * We should implement direct mtop comparison,
760 * but it might be useful to keep t_topology comparison as an option.
761 */
762 top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
763 cmp_top(stdout,&top[0],NULL,ftol,abstol);
764 }
765 }
766 }
767
768 void comp_frame(FILE *fp, t_trxframe *fr1, t_trxframe *fr2,
769 gmx_bool bRMSD, real ftol,real abstol)
770 {
771 fprintf(fp,"\n");
772 cmp_int(fp,"flags",-1,fr1->flags,fr2->flags);
773 cmp_int(fp,"not_ok",-1,fr1->not_ok,fr2->not_ok);
774 cmp_int(fp,"natoms",-1,fr1->natoms,fr2->natoms);
775 cmp_real(fp,"t0",-1,fr1->t0,fr2->t0,ftol,abstol);
776 if (cmp_bool(fp,"bTitle",-1,fr1->bTitle,fr2->bTitle))
777 cmp_str(fp,"title", -1, fr1->title, fr2->title);
778 if (cmp_bool(fp,"bStep",-1,fr1->bStep,fr2->bStep))
779 cmp_int(fp,"step",-1,fr1->step,fr2->step);
780 cmp_int(fp,"step",-1,fr1->step,fr2->step);
781 if (cmp_bool(fp,"bTime",-1,fr1->bTime,fr2->bTime))
782 cmp_real(fp,"time",-1,fr1->time,fr2->time,ftol,abstol);
783 if (cmp_bool(fp,"bLambda",-1,fr1->bLambda,fr2->bLambda))
784 cmp_real(fp,"lambda",-1,fr1->lambda,fr2->lambda,ftol,abstol);
785 if (cmp_bool(fp,"bAtoms",-1,fr1->bAtoms,fr2->bAtoms))
786 cmp_atoms(fp,fr1->atoms,fr2->atoms,ftol,abstol);
787 if (cmp_bool(fp,"bPrec",-1,fr1->bPrec,fr2->bPrec))
788 cmp_real(fp,"prec",-1,fr1->prec,fr2->prec,ftol,abstol);
789 if (cmp_bool(fp,"bX",-1,fr1->bX,fr2->bX))
790 cmp_rvecs(fp,"x",min(fr1->natoms,fr2->natoms),fr1->x,fr2->x,bRMSD,ftol,abstol);
791 if (cmp_bool(fp,"bV",-1,fr1->bV,fr2->bV))
792 cmp_rvecs(fp,"v",min(fr1->natoms,fr2->natoms),fr1->v,fr2->v,bRMSD,ftol,abstol);
793 if (cmp_bool(fp,"bF",-1,fr1->bF,fr2->bF))
794 cmp_rvecs(fp,"f",min(fr1->natoms,fr2->natoms),fr1->f,fr2->f,bRMSD,ftol,abstol);
795 if (cmp_bool(fp,"bBox",-1,fr1->bBox,fr2->bBox))
796 cmp_rvecs(fp,"box",3,fr1->box,fr2->box,FALSE,ftol,abstol);
797 }
798
799 void comp_trx(const output_env_t oenv,const char *fn1, const char *fn2,
800 gmx_bool bRMSD,real ftol,real abstol)
801 {
802 int i;
803 const char *fn[2];
804 t_trxframe fr[2];
805 t_trxstatus *status[2];
806 gmx_bool b[2];
807
808 fn[0]=fn1;
809 fn[1]=fn2;
810 fprintf(stderr,"Comparing trajectory files %s and %s\n",fn1,fn2);
811 for (i=0; i<2; i++)
812 b[i] = read_first_frame(oenv,&status[i],fn[i],&fr[i],TRX_READ_X|TRX_READ_V|TRX_READ_F);
813
814 if (b[0] && b[1]) {
815 do {
816 comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
817
818 for (i=0; i<2; i++)
819 b[i] = read_next_frame(oenv,status[i],&fr[i]);
820 } while (b[0] && b[1]);
821
822 for (i=0; i<2; i++) {
823 if (b[i] && !b[1-i])
824 fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn[i],fn[1-i]);
825 close_trj(status[i]);
826 }
827 }
828 if (!b[0] && !b[1])
829 fprintf(stdout,"\nBoth files read correctly\n");
830 }
831
832 static real ener_tensor_diag(int n,int *ind1,int *ind2,
833 gmx_enxnm_t *enm1,
834 int *tensi,int i,
835 t_energy e1[],t_energy e2[])
836 {
837 int d1,d2;
838 int len;
839 int j;
840 real prod1,prod2;
841 int nfound;
842
843 d1 = tensi[i]/DIM;
844 d2 = tensi[i] - d1*DIM;
845
846 /* Find the diagonal elements d1 and d2 */
847 len = strlen(enm1[ind1[i]].name);
848 prod1 = 1;
849 prod2 = 1;
850 nfound = 0;
851 for(j=0; j<n; j++) {
852 if (tensi[j] >= 0 &&
853 strlen(enm1[ind1[j]].name) == len &&
854 strncmp(enm1[ind1[i]].name,enm1[ind1[j]].name,len-2) == 0 &&
855 (tensi[j] == d1*DIM+d1 || tensi[j] == d2*DIM+d2)) {
856 prod1 *= fabs(e1[ind1[j]].e);
857 prod2 *= fabs(e2[ind2[j]].e);
858 nfound++;
859 }
860 }
861
862 if (nfound == 2) {
863 return 0.5*(sqrt(prod1) + sqrt(prod2));
864 } else {
865 return 0;
866 }
867 }
868
869 static gmx_bool enernm_equal(const char *nm1,const char *nm2)
870 {
871 int len1,len2;
872
873 len1 = strlen(nm1);
874 len2 = strlen(nm2);
875
876 /* Remove " (bar)" at the end of a name */
877 if (len1 > 6 && strcmp(nm1+len1-6," (bar)") == 0) {
878 len1 -= 6;
879 }
880 if (len2 > 6 && strcmp(nm2+len2-6," (bar)") == 0) {
881 len2 -= 6;
882 }
883
884 return (len1 == len2 && gmx_strncasecmp(nm1,nm2,len1) == 0);
885 }
886
887 static void cmp_energies(FILE *fp,int step1,int step2,
888 t_energy e1[],t_energy e2[],
889 gmx_enxnm_t *enm1,gmx_enxnm_t *enm2,
890 real ftol,real abstol,
891 int nre,int *ind1,int *ind2,int maxener)
892 {
893 int i,ii;
894 int *tensi,len,d1,d2;
895 real ftol_i,abstol_i;
896
897 snew(tensi,maxener);
898 /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
899 for(i=0; (i<maxener); i++) {
900 ii = ind1[i];
901 tensi[i] = -1;
902 len = strlen(enm1[ii].name);
903 if (len > 3 && enm1[ii].name[len-3] == '-') {
904 d1 = enm1[ii].name[len-2] - 'X';
905 d2 = enm1[ii].name[len-1] - 'X';
906 if (d1 >= 0 && d1 < DIM &&
907 d2 >= 0 && d2 < DIM) {
908 tensi[i] = d1*DIM + d2;
909 }
910 }
911 }
912
913 for(i=0; (i<maxener); i++) {
914 /* Check if this is an off-diagonal tensor element */
915 if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8) {
916 /* Turn on the relative tolerance check (4 is maximum relative diff.) */
917 ftol_i = 5;
918 /* Do the relative tolerance through an absolute tolerance times
919 * the size of diagonal components of the tensor.
920 */
921 abstol_i = ftol*ener_tensor_diag(nre,ind1,ind2,enm1,tensi,i,e1,e2);
922 if (debug) {
923 fprintf(debug,"tensor '%s' val %f diag %f\n",
924 enm1[i].name,e1[i].e,abstol_i/ftol);
925 }
926 if (abstol_i > 0) {
927 /* We found a diagonal, we need to check with the minimum tolerance */
928 abstol_i = min(abstol_i,abstol);
929 } else {
930 /* We did not find a diagonal, ignore the relative tolerance check */
931 abstol_i = abstol;
932 }
933 } else {
934 ftol_i = ftol;
935 abstol_i = abstol;
936 }
937 if (!equal_real(e1[ind1[i]].e,e2[ind2[i]].e,ftol_i,abstol_i)) {
938 fprintf(fp,"%-15s step %3d: %12g, step %3d: %12g\n",
939 enm1[ind1[i]].name,
940 step1,e1[ind1[i]].e,
941 step2,e2[ind2[i]].e);
942 }
943 }
944
945 sfree(tensi);
946 }
947
948 #if 0
949 static void cmp_disres(t_enxframe *fr1,t_enxframe *fr2,real ftol, real abstol)
950 {
951 int i;
952 char bav[64],bt[64],bs[22];
953
954 cmp_int(stdout,"ndisre",-1,fr1->ndisre,fr2->ndisre);
955 if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0)) {
956 sprintf(bav,"step %s: disre rav",gmx_step_str(fr1->step,bs));
957 sprintf(bt, "step %s: disre rt",gmx_step_str(fr1->step,bs));
958 for(i=0; (i<fr1->ndisre); i++) {
959 cmp_real(stdout,bav,i,fr1->disre_rm3tav[i],fr2->disre_rm3tav[i],ftol,abstol);
960 cmp_real(stdout,bt ,i,fr1->disre_rt[i] ,fr2->disre_rt[i] ,ftol,abstol);
961 }
962 }
963 }
964 #endif
965
966 static void cmp_eblocks(t_enxframe *fr1,t_enxframe *fr2,real ftol, real abstol)
967 {
968 int i,j,k;
969 char buf[64],bs[22];
970
971 cmp_int(stdout,"nblock",-1,fr1->nblock,fr2->nblock);
972 if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
973 {
974 for(j=0; (j<fr1->nblock); j++)
975 {
976 t_enxblock *b1, *b2; /* convenience vars */
977
978 b1=&(fr1->block[j]);
979 b2=&(fr2->block[j]);
980
981 sprintf(buf,"step %s: block[%d]",gmx_step_str(fr1->step,bs),j);
982 cmp_int(stdout,buf,-1,b1->nsub,b2->nsub);
983 cmp_int(stdout,buf,-1,b1->id,b2->id);
984
985 if ( (b1->nsub==b2->nsub) && (b1->id==b2->id) )
986 {
987 for(i=0;i<b1->nsub;i++)
988 {
989 t_enxsubblock *s1, *s2;
990
991 s1=&(b1->sub[i]);
992 s2=&(b2->sub[i]);
993
994 cmp_int(stdout, buf, -1, (int)s1->type, (int)s2->type);
995 cmp_gmx_large_int(stdout, buf, s1->nr, s2->nr);
996
997 if ((s1->type == s2->type) && (s1->nr == s2->nr))
998 {
999 switch(s1->type)
1000 {
1001 case xdr_datatype_float:
1002 for(k=0;k<s1->nr;k++)
1003 {
1004 cmp_float(stdout, buf, i,
1005 s1->fval[k], s2->fval[k],
1006 ftol, abstol);
1007 }
1008 break;
1009 case xdr_datatype_double:
1010 for(k=0;k<s1->nr;k++)
1011 {
1012 cmp_double(stdout, buf, i,
1013 s1->dval[k], s2->dval[k],
1014 ftol, abstol);
1015 }
1016 break;
1017 case xdr_datatype_int:
1018 for(k=0;k<s1->nr;k++)
1019 {
1020 cmp_int(stdout, buf, i,
1021 s1->ival[k], s2->ival[k]);
1022 }
1023 break;
1024 case xdr_datatype_large_int:
1025 for(k=0;k<s1->nr;k++)
1026 {
1027 cmp_gmx_large_int(stdout, buf,
1028 s1->lval[k], s2->lval[k]);
1029 }
1030 break;
1031 case xdr_datatype_char:
1032 for(k=0;k<s1->nr;k++)
1033 {
1034 cmp_uc(stdout, buf, i,
1035 s1->cval[k], s2->cval[k]);
1036 }
1037 break;
1038 case xdr_datatype_string:
1039 for(k=0;k<s1->nr;k++)
1040 {
1041 cmp_str(stdout, buf, i,
1042 s1->sval[k], s2->sval[k]);
1043 }
1044 break;
1045 default:
1046 gmx_incons("Unknown data type!!");
1047 }
1048 }
1049 }
1050 }
1051 }
1052 }
1053 }
1054
1055 void comp_enx(const char *fn1,const char *fn2,real ftol,real abstol,const char *lastener)
1056 {
1057 int nre,nre1,nre2,block;
1058 ener_file_t in1, in2;
1059 int i,j,maxener,*ind1,*ind2,*have;
1060 char buf[256];
1061 gmx_enxnm_t *enm1=NULL,*enm2=NULL;
1062 t_enxframe *fr1,*fr2;
1063 gmx_bool b1,b2;
1064
1065 fprintf(stdout,"comparing energy file %s and %s\n\n",fn1,fn2);
1066
1067 in1 = open_enx(fn1,"r");
1068 in2 = open_enx(fn2,"r");
1069 do_enxnms(in1,&nre1,&enm1);
1070 do_enxnms(in2,&nre2,&enm2);
1071 if (nre1 != nre2) {
1072 fprintf(stdout,"There are %d and %d terms in the energy files\n\n",
1073 nre1,nre2);
1074 } else {
1075 fprintf(stdout,"There are %d terms in the energy files\n\n",nre1);
1076 }
1077
1078 snew(ind1,nre1);
1079 snew(ind2,nre2);
1080 snew(have,nre2);
1081 nre = 0;
1082 for(i=0; i<nre1; i++) {
1083 for(j=0; j<nre2; j++) {
1084 if (enernm_equal(enm1[i].name,enm2[j].name)) {
1085 ind1[nre] = i;
1086 ind2[nre] = j;
1087 have[j] = 1;
1088 nre++;
1089 break;
1090 }
1091 }
1092 if (nre == 0 || ind1[nre-1] != i) {
1093 cmp_str(stdout,"enm",i,enm1[i].name,"-");
1094 }
1095 }
1096 for(i=0; i<nre2; i++) {
1097 if (have[i] == 0) {
1098 cmp_str(stdout,"enm",i,"-",enm2[i].name);
1099 }
1100 }
1101
1102 maxener = nre;
1103 for(i=0; i<nre; i++) {
1104 if ((lastener != NULL) && (strstr(enm1[i].name,lastener) != NULL)) {
1105 maxener=i+1;
1106 break;
1107 }
1108 }
1109
1110 fprintf(stdout,"There are %d terms to compare in the energy files\n\n",
1111 maxener);
1112
1113 for(i=0; i<maxener; i++) {
1114 cmp_str(stdout,"unit",i,enm1[ind1[i]].unit,enm2[ind2[i]].unit);
1115 }
1116
1117 snew(fr1,1);
1118 snew(fr2,1);
1119 do {
1120 b1 = do_enx(in1,fr1);
1121 b2 = do_enx(in2,fr2);
1122 if (b1 && !b2)
1123 fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn2,fn1);
1124 else if (!b1 && b2)
1125 fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn1,fn2);
1126 else if (!b1 && !b2)
1127 fprintf(stdout,"\nFiles read successfully\n");
1128 else {
1129 cmp_real(stdout,"t",-1,fr1->t,fr2->t,ftol,abstol);
1130 cmp_int(stdout,"step",-1,fr1->step,fr2->step);
1131 /* We don't want to print the nre mismatch for every frame */
1132 /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1133 if ((fr1->nre >= nre) && (fr2->nre >= nre))
1134 cmp_energies(stdout,fr1->step,fr1->step,fr1->ener,fr2->ener,
1135 enm1,enm2,ftol,abstol,nre,ind1,ind2,maxener);
1136 /*cmp_disres(fr1,fr2,ftol,abstol);*/
1137 cmp_eblocks(fr1,fr2,ftol,abstol);
1138 }
1139 } while (b1 && b2);
1140
1141 close_enx(in1);
1142 close_enx(in2);
1143
1144 free_enxframe(fr2);
1145 sfree(fr2);
1146 free_enxframe(fr1);
1147 sfree(fr1);
1148 }
got merge a long time ago