?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
merge with upstream gromacs
[gromacs/adressmacs.git] / src / kernel / tpbcmp.c
blob15353089a1fdf222acbb9d6bf44b8ecdc012f8a9
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
56 static void cmp_int(FILE *fp,const char *s,int index,int i1,int i2)
57 {
58 if (i1 != i2) {
59 if (index != -1)
60 fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
61 else
62 fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
63 }
64 }
65
66 static void cmp_gmx_large_int(FILE *fp,const char *s,gmx_large_int_t i1,gmx_large_int_t i2)
67 {
68 if (i1 != i2) {
69 fprintf(fp,"%s (",s);
70 fprintf(fp,gmx_large_int_pfmt,i1);
71 fprintf(fp," - ");
72 fprintf(fp,gmx_large_int_pfmt,i2);
73 fprintf(fp,")\n");
74 }
75 }
76
77 static void cmp_us(FILE *fp,const char *s,int index,unsigned short i1,unsigned short i2)
78 {
79 if (i1 != i2) {
80 if (index != -1)
81 fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
82 else
83 fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
84 }
85 }
86
87 static void cmp_uc(FILE *fp,const char *s,int index,unsigned char i1,unsigned char i2)
88 {
89 if (i1 != i2) {
90 if (index != -1)
91 fprintf(fp,"%s[%d] (%d - %d)\n",s,index,i1,i2);
92 else
93 fprintf(fp,"%s (%d - %d)\n",s,i1,i2);
94 }
95 }
96
97 static bool cmp_bool(FILE *fp, const char *s, int index, bool b1, bool b2)
98 {
99 if (b1) {
100 b1 = 1;
101 } else {
102 b1 = 0;
103 }
104 if (b2) {
105 b2 = 1;
106 } else {
107 b2 = 0;
108 }
109 if (b1 != b2) {
110 if (index != -1)
111 fprintf(fp,"%s[%d] (%s - %s)\n",s,index,
112 bool_names[b1],bool_names[b2]);
113 else
114 fprintf(fp,"%s (%s - %s)\n",s,
115 bool_names[b1],bool_names[b2]);
116 }
117 return b1 && b2;
118 }
119
120 static void cmp_str(FILE *fp, const char *s, int index,
121 const char *s1, const char *s2)
122 {
123 if (strcmp(s1,s2) != 0) {
124 if (index != -1)
125 fprintf(fp,"%s[%d] (%s - %s)\n",s,index,s1,s2);
126 else
127 fprintf(fp,"%s (%s - %s)\n",s,s1,s2);
128 }
129 }
130
131 static bool equal_real(real i1,real i2,real ftol,real abstol)
132 {
133 return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
134 }
135
136 static bool equal_double(double i1,double i2,real ftol,real abstol)
137 {
138 return ( ( 2*fabs(i1 - i2) <= (fabs(i1) + fabs(i2))*ftol ) || fabs(i1-i2)<=abstol );
139 }
140
141 static void
142 cmp_real(FILE *fp,const char *s,int index,real i1,real i2,real ftol,real abstol)
143 {
144 if (!equal_real(i1,i2,ftol,abstol)) {
145 if (index != -1)
146 fprintf(fp,"%s[%2d] (%e - %e)\n",s,index,i1,i2);
147 else
148 fprintf(fp,"%s (%e - %e)\n",s,i1,i2);
149 }
150 }
151
152 static void
153 cmp_double(FILE *fp,const char *s,int index,double i1,double i2,double ftol,double abstol)
154 {
155 if (!equal_double(i1,i2,ftol,abstol)) {
156 if (index != -1)
157 fprintf(fp,"%s[%2d] (%16.9e - %16.9e)\n",s,index,i1,i2);
158 else
159 fprintf(fp,"%s (%16.9e - %16.9e)\n",s,i1,i2);
160 }
161 }
162
163 static void cmp_rvec(FILE *fp,const char *s,int index,rvec i1,rvec i2,real ftol,real abstol)
164 {
165 if(!equal_real(i1[XX],i2[XX],ftol,abstol) ||
166 !equal_real(i1[YY],i2[YY],ftol,abstol) ||
167 !equal_real(i1[ZZ],i2[ZZ],ftol,abstol))
168 {
169 if (index != -1)
170 fprintf(fp,"%s[%5d] (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
171 s,index,i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
172 else
173 fprintf(fp,"%s (%12.5e %12.5e %12.5e) - (%12.5e %12.5e %12.5e)\n",
174 s,i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
175 }
176 }
177
178 static void cmp_ivec(FILE *fp,const char *s,int index,ivec i1,ivec i2)
179 {
180 if ((i1[XX] != i2[XX]) || (i1[YY] != i2[YY]) || (i1[ZZ] != i2[ZZ])) {
181 if (index != -1)
182 fprintf(fp,"%s[%5d] (%8d,%8d,%8d - %8d,%8d,%8d)\n",s,index,
183 i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
184 else
185 fprintf(fp,"%s (%8d,%8d,%8d - %8d,%8d,%8d)\n",s,
186 i1[XX],i1[YY],i1[ZZ],i2[XX],i2[YY],i2[ZZ]);
187 }
188 }
189
190 static void cmp_ilist(FILE *fp,int ftype,t_ilist *il1,t_ilist *il2)
191 {
192 int i;
193 char buf[256];
194
195 fprintf(fp,"comparing ilist %s\n",interaction_function[ftype].name);
196 sprintf(buf,"%s->nr",interaction_function[ftype].name);
197 cmp_int(fp,buf,-1,il1->nr,il2->nr);
198 sprintf(buf,"%s->iatoms",interaction_function[ftype].name);
199 if (((il1->nr > 0) && (!il1->iatoms)) ||
200 ((il2->nr > 0) && (!il2->iatoms)) ||
201 ((il1->nr != il2->nr)))
202 fprintf(fp,"Comparing radically different topologies - %s is different\n",
203 buf);
204 else
205 for(i=0; (i<il1->nr); i++)
206 cmp_int(fp,buf,i,il1->iatoms[i],il2->iatoms[i]);
207 }
208
209 void cmp_iparm(FILE *fp,const char *s,t_functype ft,
210 t_iparams ip1,t_iparams ip2,real ftol,real abstol)
211 {
212 int i;
213 bool bDiff;
214
215 bDiff=FALSE;
216 for(i=0; i<MAXFORCEPARAM && !bDiff; i++)
217 bDiff = !equal_real(ip1.generic.buf[i],ip2.generic.buf[i],ftol,abstol);
218 if (bDiff) {
219 fprintf(fp,"%s1: ",s);
220 pr_iparams(fp,ft,&ip1);
221 fprintf(fp,"%s2: ",s);
222 pr_iparams(fp,ft,&ip2);
223 }
224 }
225
226 void cmp_iparm_AB(FILE *fp,const char *s,t_functype ft,t_iparams ip1,real ftol,real abstol)
227 {
228 int nrfpA,nrfpB,p0,i;
229 bool bDiff;
230
231 /* Normally the first parameter is perturbable */
232 p0 = 0;
233 nrfpA = interaction_function[ft].nrfpA;
234 nrfpB = interaction_function[ft].nrfpB;
235 if (ft == F_PDIHS) {
236 nrfpB = 2;
237 } else if (interaction_function[ft].flags & IF_TABULATED) {
238 /* For tabulated interactions only the second parameter is perturbable */
239 p0 = 1;
240 nrfpB = 1;
241 }
242 bDiff=FALSE;
243 for(i=0; i<nrfpB && !bDiff; i++) {
244 bDiff = !equal_real(ip1.generic.buf[p0+i],ip1.generic.buf[nrfpA+i],ftol,abstol);
245 }
246 if (bDiff) {
247 fprintf(fp,"%s: ",s);
248 pr_iparams(fp,ft,&ip1);
249 }
250 }
251
252 static void cmp_idef(FILE *fp,t_idef *id1,t_idef *id2,real ftol,real abstol)
253 {
254 int i;
255 char buf1[64],buf2[64];
256
257 fprintf(fp,"comparing idef\n");
258 if (id2) {
259 cmp_int(fp,"idef->ntypes",-1,id1->ntypes,id2->ntypes);
260 cmp_int(fp,"idef->atnr", -1,id1->atnr,id2->atnr);
261 for(i=0; (i<id1->ntypes); i++) {
262 sprintf(buf1,"idef->functype[%d]",i);
263 sprintf(buf2,"idef->iparam[%d]",i);
264 cmp_int(fp,buf1,i,(int)id1->functype[i],(int)id2->functype[i]);
265 cmp_iparm(fp,buf2,id1->functype[i],
266 id1->iparams[i],id2->iparams[i],ftol,abstol);
267 }
268 cmp_real(fp,"fudgeQQ",-1,id1->fudgeQQ,id2->fudgeQQ,ftol,abstol);
269 for(i=0; (i<F_NRE); i++)
270 cmp_ilist(fp,i,&(id1->il[i]),&(id2->il[i]));
271 } else {
272 for(i=0; (i<id1->ntypes); i++)
273 cmp_iparm_AB(fp,"idef->iparam",id1->functype[i],id1->iparams[i],ftol,abstol);
274 }
275 }
276
277 static void cmp_block(FILE *fp,t_block *b1,t_block *b2,const char *s)
278 {
279 int i,j,k;
280 char buf[32];
281
282 fprintf(fp,"comparing block %s\n",s);
283 sprintf(buf,"%s.nr",s);
284 cmp_int(fp,buf,-1,b1->nr,b2->nr);
285 }
286
287 static void cmp_blocka(FILE *fp,t_blocka *b1,t_blocka *b2,const char *s)
288 {
289 int i,j,k;
290 char buf[32];
291
292 fprintf(fp,"comparing blocka %s\n",s);
293 sprintf(buf,"%s.nr",s);
294 cmp_int(fp,buf,-1,b1->nr,b2->nr);
295 sprintf(buf,"%s.nra",s);
296 cmp_int(fp,buf,-1,b1->nra,b2->nra);
297 }
298
299 static void cmp_atom(FILE *fp,int index,t_atom *a1,t_atom *a2,real ftol,real abstol)
300 {
301 int i;
302 char buf[256];
303
304 if (a2) {
305 cmp_us(fp,"atom.type",index,a1->type,a2->type);
306 cmp_us(fp,"atom.ptype",index,a1->ptype,a2->ptype);
307 cmp_int(fp,"atom.resind",index,a1->resind,a2->resind);
308 cmp_int(fp,"atom.atomnumber",index,a1->atomnumber,a2->atomnumber);
309 cmp_real(fp,"atom.m",index,a1->m,a2->m,ftol,abstol);
310 cmp_real(fp,"atom.q",index,a1->q,a2->q,ftol,abstol);
311 cmp_us(fp,"atom.typeB",index,a1->typeB,a2->typeB);
312 cmp_real(fp,"atom.mB",index,a1->mB,a2->mB,ftol,abstol);
313 cmp_real(fp,"atom.qB",index,a1->qB,a2->qB,ftol,abstol);
314 } else {
315 cmp_us(fp,"atom.type",index,a1->type,a1->typeB);
316 cmp_real(fp,"atom.m",index,a1->m,a1->mB,ftol,abstol);
317 cmp_real(fp,"atom.q",index,a1->q,a1->qB,ftol,abstol);
318 }
319 }
320
321 static void cmp_atoms(FILE *fp,t_atoms *a1,t_atoms *a2,real ftol, real abstol)
322 {
323 int i;
324
325 fprintf(fp,"comparing atoms\n");
326
327 if (a2) {
328 cmp_int(fp,"atoms->nr",-1,a1->nr,a2->nr);
329 for(i=0; (i<a1->nr); i++)
330 cmp_atom(fp,i,&(a1->atom[i]),&(a2->atom[i]),ftol,abstol);
331 } else {
332 for(i=0; (i<a1->nr); i++)
333 cmp_atom(fp,i,&(a1->atom[i]),NULL,ftol,abstol);
334 }
335 }
336
337 static void cmp_top(FILE *fp,t_topology *t1,t_topology *t2,real ftol, real abstol)
338 {
339 int i;
340
341 fprintf(fp,"comparing top\n");
342 if (t2) {
343 cmp_idef(fp,&(t1->idef),&(t2->idef),ftol,abstol);
344 cmp_atoms(fp,&(t1->atoms),&(t2->atoms),ftol,abstol);
345 cmp_block(fp,&t1->cgs,&t2->cgs,"cgs");
346 cmp_block(fp,&t1->mols,&t2->mols,"mols");
347 cmp_blocka(fp,&t1->excls,&t2->excls,"excls");
348 } else {
349 cmp_idef(fp,&(t1->idef),NULL,ftol,abstol);
350 cmp_atoms(fp,&(t1->atoms),NULL,ftol,abstol);
351 }
352 }
353
354 static void cmp_rvecs(FILE *fp,const char *title,int n,rvec x1[],rvec x2[],
355 bool bRMSD,real ftol,real abstol)
356 {
357 int i,m;
358 double d,ssd;
359
360 if (bRMSD) {
361 ssd = 0;
362 for(i=0; (i<n); i++) {
363 for(m=0; m<DIM; m++) {
364 d = x1[i][m] - x2[i][m];
365 ssd += d*d;
366 }
367 }
368 fprintf(fp,"%s RMSD %f\n",title,sqrt(ssd/n));
369 } else {
370 for(i=0; (i<n); i++) {
371 cmp_rvec(fp,title,i,x1[i],x2[i],ftol,abstol);
372 }
373 }
374 }
375
376 static void cmp_grpopts(FILE *fp,t_grpopts *opt1,t_grpopts *opt2,real ftol, real abstol)
377 {
378 int i,j;
379 char buf1[256],buf2[256];
380
381 cmp_int(fp,"inputrec->grpopts.ngtc",-1, opt1->ngtc,opt2->ngtc);
382 cmp_int(fp,"inputrec->grpopts.ngacc",-1, opt1->ngacc,opt2->ngacc);
383 cmp_int(fp,"inputrec->grpopts.ngfrz",-1, opt1->ngfrz,opt2->ngfrz);
384 cmp_int(fp,"inputrec->grpopts.ngener",-1,opt1->ngener,opt2->ngener);
385 for(i=0; (i<min(opt1->ngtc,opt2->ngtc)); i++) {
386 cmp_real(fp,"inputrec->grpopts.nrdf",i,opt1->nrdf[i],opt2->nrdf[i],ftol,abstol);
387 cmp_real(fp,"inputrec->grpopts.ref_t",i,opt1->ref_t[i],opt2->ref_t[i],ftol,abstol);
388 cmp_real(fp,"inputrec->grpopts.tau_t",i,opt1->tau_t[i],opt2->tau_t[i],ftol,abstol);
389 cmp_int(fp,"inputrec->grpopts.annealing",i,opt1->annealing[i],opt2->annealing[i]);
390 cmp_int(fp,"inputrec->grpopts.anneal_npoints",i,
391 opt1->anneal_npoints[i],opt2->anneal_npoints[i]);
392 if(opt1->anneal_npoints[i]==opt2->anneal_npoints[i]) {
393 sprintf(buf1,"inputrec->grpopts.anneal_time[%d]",i);
394 sprintf(buf2,"inputrec->grpopts.anneal_temp[%d]",i);
395 for(j=0;j<opt1->anneal_npoints[i];j++) {
396 cmp_real(fp,buf1,j,opt1->anneal_time[i][j],opt2->anneal_time[i][j],ftol,abstol);
397 cmp_real(fp,buf2,j,opt1->anneal_temp[i][j],opt2->anneal_temp[i][j],ftol,abstol);
398 }
399 }
400 }
401 if (opt1->ngener == opt2->ngener) {
402 for(i=0; i<opt1->ngener; i++)
403 for(j=i; j<opt1->ngener; j++) {
404 sprintf(buf1,"inputrec->grpopts.egp_flags[%d]",i);
405 cmp_int(fp,buf1,j,
406 opt1->egp_flags[opt1->ngener*i+j],
407 opt2->egp_flags[opt1->ngener*i+j]);
408 }
409 }
410 for(i=0; (i<min(opt1->ngacc,opt2->ngacc)); i++)
411 cmp_rvec(fp,"inputrec->grpopts.acc",i,opt1->acc[i],opt2->acc[i],ftol,abstol);
412 for(i=0; (i<min(opt1->ngfrz,opt2->ngfrz)); i++)
413 cmp_ivec(fp,"inputrec->grpopts.nFreeze",i,opt1->nFreeze[i],opt2->nFreeze[i]);
414 }
415
416 static void cmp_cosines(FILE *fp,const char *s,t_cosines c1[DIM],t_cosines c2[DIM],real ftol, real abstol)
417 {
418 int i,m;
419 char buf[256];
420
421 for(m=0; (m<DIM); m++) {
422 sprintf(buf,"inputrec->%s[%d]",s,m);
423 cmp_int(fp,buf,0,c1->n,c2->n);
424 for(i=0; (i<min(c1->n,c2->n)); i++) {
425 cmp_real(fp,buf,i,c1->a[i],c2->a[i],ftol,abstol);
426 cmp_real(fp,buf,i,c1->phi[i],c2->phi[i],ftol,abstol);
427 }
428 }
429 }
430
431 static void cmp_pull(FILE *fp,t_pull *pull1,t_pull *pull2,real ftol, real abstol)
432 {
433 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");
434 }
435
436 static void cmp_inputrec(FILE *fp,t_inputrec *ir1,t_inputrec *ir2,real ftol, real abstol)
437 {
438 int i;
439
440 fprintf(fp,"comparing inputrec\n");
441
442 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
443 * of warnings. Maybe it will change in future versions, but for the
444 * moment I've spelled them out instead. /EL 000820
445 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
446 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
447 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
448 */
449 cmp_int(fp,"inputrec->eI",-1,ir1->eI,ir2->eI);
450 cmp_gmx_large_int(fp,"inputrec->nsteps",ir1->nsteps,ir2->nsteps);
451 cmp_gmx_large_int(fp,"inputrec->init_step",ir1->init_step,ir2->init_step);
452 cmp_int(fp,"inputrec->simulation_part",-1,ir1->simulation_part,ir2->simulation_part);
453 cmp_int(fp,"inputrec->nstcalcenergy",-1,ir1->nstcalcenergy,ir2->nstcalcenergy);
454 cmp_int(fp,"inputrec->ePBC",-1,ir1->ePBC,ir2->ePBC);
455 cmp_int(fp,"inputrec->bPeriodicMols",-1,ir1->bPeriodicMols,ir2->bPeriodicMols);
456 cmp_int(fp,"inputrec->ns_type",-1,ir1->ns_type,ir2->ns_type);
457 cmp_int(fp,"inputrec->nstlist",-1,ir1->nstlist,ir2->nstlist);
458 cmp_int(fp,"inputrec->ndelta",-1,ir1->ndelta,ir2->ndelta);
459 cmp_int(fp,"inputrec->nstcomm",-1,ir1->nstcomm,ir2->nstcomm);
460 cmp_int(fp,"inputrec->comm_mode",-1,ir1->comm_mode,ir2->comm_mode);
461 cmp_int(fp,"inputrec->nstcheckpoint",-1,ir1->nstcheckpoint,ir2->nstcheckpoint);
462 cmp_int(fp,"inputrec->nstlog",-1,ir1->nstlog,ir2->nstlog);
463 cmp_int(fp,"inputrec->nstxout",-1,ir1->nstxout,ir2->nstxout);
464 cmp_int(fp,"inputrec->nstvout",-1,ir1->nstvout,ir2->nstvout);
465 cmp_int(fp,"inputrec->nstfout",-1,ir1->nstfout,ir2->nstfout);
466 cmp_int(fp,"inputrec->nstenergy",-1,ir1->nstenergy,ir2->nstenergy);
467 cmp_int(fp,"inputrec->nstxtcout",-1,ir1->nstxtcout,ir2->nstxtcout);
468 cmp_double(fp,"inputrec->init_t",-1,ir1->init_t,ir2->init_t,ftol,abstol);
469 cmp_double(fp,"inputrec->delta_t",-1,ir1->delta_t,ir2->delta_t,ftol,abstol);
470 cmp_real(fp,"inputrec->xtcprec",-1,ir1->xtcprec,ir2->xtcprec,ftol,abstol);
471 cmp_int(fp,"inputrec->nkx",-1,ir1->nkx,ir2->nkx);
472 cmp_int(fp,"inputrec->nky",-1,ir1->nky,ir2->nky);
473 cmp_int(fp,"inputrec->nkz",-1,ir1->nkz,ir2->nkz);
474 cmp_int(fp,"inputrec->pme_order",-1,ir1->pme_order,ir2->pme_order);
475 cmp_real(fp,"inputrec->ewald_rtol",-1,ir1->ewald_rtol,ir2->ewald_rtol,ftol,abstol);
476 cmp_int(fp,"inputrec->ewald_geometry",-1,ir1->ewald_geometry,ir2->ewald_geometry);
477 cmp_real(fp,"inputrec->epsilon_surface",-1,ir1->epsilon_surface,ir2->epsilon_surface,ftol,abstol);
478 cmp_int(fp,"inputrec->bOptFFT",-1,ir1->bOptFFT,ir2->bOptFFT);
479 cmp_int(fp,"inputrec->bContinuation",-1,ir1->bContinuation,ir2->bContinuation);
480 cmp_int(fp,"inputrec->bShakeSOR",-1,ir1->bShakeSOR,ir2->bShakeSOR);
481 cmp_int(fp,"inputrec->etc",-1,ir1->etc,ir2->etc);
482 cmp_int(fp,"inputrec->epc",-1,ir1->epc,ir2->epc);
483 cmp_int(fp,"inputrec->epct",-1,ir1->epct,ir2->epct);
484 cmp_real(fp,"inputrec->tau_p",-1,ir1->tau_p,ir2->tau_p,ftol,abstol);
485 cmp_rvec(fp,"inputrec->ref_p(x)",-1,ir1->ref_p[XX],ir2->ref_p[XX],ftol,abstol);
486 cmp_rvec(fp,"inputrec->ref_p(y)",-1,ir1->ref_p[YY],ir2->ref_p[YY],ftol,abstol);
487 cmp_rvec(fp,"inputrec->ref_p(z)",-1,ir1->ref_p[ZZ],ir2->ref_p[ZZ],ftol,abstol);
488 cmp_rvec(fp,"inputrec->compress(x)",-1,ir1->compress[XX],ir2->compress[XX],ftol,abstol);
489 cmp_rvec(fp,"inputrec->compress(y)",-1,ir1->compress[YY],ir2->compress[YY],ftol,abstol);
490 cmp_rvec(fp,"inputrec->compress(z)",-1,ir1->compress[ZZ],ir2->compress[ZZ],ftol,abstol);
491 cmp_int(fp,"refcoord_scaling",-1,ir1->refcoord_scaling,ir2->refcoord_scaling);
492 cmp_rvec(fp,"inputrec->posres_com",-1,ir1->posres_com,ir2->posres_com,ftol,abstol);
493 cmp_rvec(fp,"inputrec->posres_comB",-1,ir1->posres_comB,ir2->posres_comB,ftol,abstol);
494 cmp_int(fp,"inputrec->andersen_seed",-1,ir1->andersen_seed,ir2->andersen_seed);
495 cmp_real(fp,"inputrec->rlist",-1,ir1->rlist,ir2->rlist,ftol,abstol);
496 cmp_real(fp,"inputrec->rlistlong",-1,ir1->rlistlong,ir2->rlistlong,ftol,abstol);
497 cmp_real(fp,"inputrec->rtpi",-1,ir1->rtpi,ir2->rtpi,ftol,abstol);
498 cmp_int(fp,"inputrec->coulombtype",-1,ir1->coulombtype,ir2->coulombtype);
499 cmp_real(fp,"inputrec->rcoulomb_switch",-1,ir1->rcoulomb_switch,ir2->rcoulomb_switch,ftol,abstol);
500 cmp_real(fp,"inputrec->rcoulomb",-1,ir1->rcoulomb,ir2->rcoulomb,ftol,abstol);
501 cmp_int(fp,"inputrec->vdwtype",-1,ir1->vdwtype,ir2->vdwtype);
502 cmp_real(fp,"inputrec->rvdw_switch",-1,ir1->rvdw_switch,ir2->rvdw_switch,ftol,abstol);
503 cmp_real(fp,"inputrec->rvdw",-1,ir1->rvdw,ir2->rvdw,ftol,abstol);
504 cmp_real(fp,"inputrec->epsilon_r",-1,ir1->epsilon_r,ir2->epsilon_r,ftol,abstol);
505 cmp_real(fp,"inputrec->epsilon_rf",-1,ir1->epsilon_rf,ir2->epsilon_rf,ftol,abstol);
506 cmp_real(fp,"inputrec->tabext",-1,ir1->tabext,ir2->tabext,ftol,abstol);
507 cmp_int(fp,"inputrec->implicit_solvent",-1,ir1->implicit_solvent,ir2->implicit_solvent);
508 cmp_int(fp,"inputrec->gb_algorithm",-1,ir1->gb_algorithm,ir2->gb_algorithm);
509 cmp_int(fp,"inputrec->nstgbradii",-1,ir1->nstgbradii,ir2->nstgbradii);
510 cmp_real(fp,"inputrec->rgbradii",-1,ir1->rgbradii,ir2->rgbradii,ftol,abstol);
511 cmp_real(fp,"inputrec->gb_saltconc",-1,ir1->gb_saltconc,ir2->gb_saltconc,ftol,abstol);
512 cmp_real(fp,"inputrec->gb_epsilon_solvent",-1,ir1->gb_epsilon_solvent,ir2->gb_epsilon_solvent,ftol,abstol);
513 cmp_real(fp,"inputrec->gb_obc_alpha",-1,ir1->gb_obc_alpha,ir2->gb_obc_alpha,ftol,abstol);
514 cmp_real(fp,"inputrec->gb_obc_beta",-1,ir1->gb_obc_beta,ir2->gb_obc_beta,ftol,abstol);
515 cmp_real(fp,"inputrec->gb_obc_gamma",-1,ir1->gb_obc_gamma,ir2->gb_obc_gamma,ftol,abstol);
516 cmp_real(fp,"inputrec->gb_dielectric_offset",-1,ir1->gb_dielectric_offset,ir2->gb_dielectric_offset,ftol,abstol);
517 cmp_int(fp,"inputrec->sa_algorithm",-1,ir1->sa_algorithm,ir2->sa_algorithm);
518 cmp_real(fp,"inputrec->sa_surface_tension",-1,ir1->sa_surface_tension,ir2->sa_surface_tension,ftol,abstol);
519
520 cmp_int(fp,"inputrec->eDispCorr",-1,ir1->eDispCorr,ir2->eDispCorr);
521 cmp_real(fp,"inputrec->shake_tol",-1,ir1->shake_tol,ir2->shake_tol,ftol,abstol);
522 cmp_int(fp,"inputrec->efep",-1,ir1->efep,ir2->efep);
523 cmp_double(fp,"inputrec->init_lambda",-1,ir1->init_lambda,ir2->init_lambda,ftol,abstol);
524 cmp_double(fp,"inputrec->delta_lambda",-1,ir1->delta_lambda,ir2->delta_lambda,ftol,abstol);
525 cmp_int(fp,"inputrec->n_foreign_lambda",-1,ir1->n_flambda,ir2->n_flambda);
526 for(i=0; i<min(ir1->n_flambda,ir2->n_flambda); i++) {
527 cmp_double(fp,"inputrec->foreign_lambda",-1,ir1->flambda[i],ir2->flambda[i],ftol,abstol);
528 }
529 cmp_real(fp,"inputrec->sc_alpha",-1,ir1->sc_alpha,ir2->sc_alpha,ftol,abstol);
530 cmp_int(fp,"inputrec->sc_power",-1,ir1->sc_power,ir2->sc_power);
531 cmp_real(fp,"inputrec->sc_sigma",-1,ir1->sc_sigma,ir2->sc_sigma,ftol,abstol);
532 cmp_int(fp,"inputrec->nstdhdl",-1,ir1->nstdhdl,ir2->nstdhdl);
533
534 cmp_int(fp,"inputrec->nwall",-1,ir1->nwall,ir2->nwall);
535 cmp_int(fp,"inputrec->wall_type",-1,ir1->wall_type,ir2->wall_type);
536 cmp_int(fp,"inputrec->wall_atomtype[0]",-1,ir1->wall_atomtype[0],ir2->wall_atomtype[0]);
537 cmp_int(fp,"inputrec->wall_atomtype[1]",-1,ir1->wall_atomtype[1],ir2->wall_atomtype[1]);
538 cmp_real(fp,"inputrec->wall_density[0]",-1,ir1->wall_density[0],ir2->wall_density[0],ftol,abstol);
539 cmp_real(fp,"inputrec->wall_density[1]",-1,ir1->wall_density[1],ir2->wall_density[1],ftol,abstol);
540 cmp_real(fp,"inputrec->wall_ewald_zfac",-1,ir1->wall_ewald_zfac,ir2->wall_ewald_zfac,ftol,abstol);
541
542 cmp_int(fp,"inputrec->ePull",-1,ir1->ePull,ir2->ePull);
543 if (ir1->ePull == ir2->ePull && ir1->ePull != epullNO)
544 cmp_pull(fp,ir1->pull,ir2->pull,ftol,abstol);
545
546 cmp_int(fp,"inputrec->eDisre",-1,ir1->eDisre,ir2->eDisre);
547 cmp_real(fp,"inputrec->dr_fc",-1,ir1->dr_fc,ir2->dr_fc,ftol,abstol);
548 cmp_int(fp,"inputrec->eDisreWeighting",-1,ir1->eDisreWeighting,ir2->eDisreWeighting);
549 cmp_int(fp,"inputrec->bDisreMixed",-1,ir1->bDisreMixed,ir2->bDisreMixed);
550 cmp_int(fp,"inputrec->nstdisreout",-1,ir1->nstdisreout,ir2->nstdisreout);
551 cmp_real(fp,"inputrec->dr_tau",-1,ir1->dr_tau,ir2->dr_tau,ftol,abstol);
552 cmp_real(fp,"inputrec->orires_fc",-1,ir1->orires_fc,ir2->orires_fc,ftol,abstol);
553 cmp_real(fp,"inputrec->orires_tau",-1,ir1->orires_tau,ir2->orires_tau,ftol,abstol);
554 cmp_int(fp,"inputrec->nstorireout",-1,ir1->nstorireout,ir2->nstorireout);
555 cmp_real(fp,"inputrec->dihre_fc",-1,ir1->dihre_fc,ir2->dihre_fc,ftol,abstol);
556 cmp_real(fp,"inputrec->em_stepsize",-1,ir1->em_stepsize,ir2->em_stepsize,ftol,abstol);
557 cmp_real(fp,"inputrec->em_tol",-1,ir1->em_tol,ir2->em_tol,ftol,abstol);
558 cmp_int(fp,"inputrec->niter",-1,ir1->niter,ir2->niter);
559 cmp_real(fp,"inputrec->fc_stepsize",-1,ir1->fc_stepsize,ir2->fc_stepsize,ftol,abstol);
560 cmp_int(fp,"inputrec->nstcgsteep",-1,ir1->nstcgsteep,ir2->nstcgsteep);
561 cmp_int(fp,"inputrec->nbfgscorr",0,ir1->nbfgscorr,ir2->nbfgscorr);
562 cmp_int(fp,"inputrec->eConstrAlg",-1,ir1->eConstrAlg,ir2->eConstrAlg);
563 cmp_int(fp,"inputrec->nProjOrder",-1,ir1->nProjOrder,ir2->nProjOrder);
564 cmp_real(fp,"inputrec->LincsWarnAngle",-1,ir1->LincsWarnAngle,ir2->LincsWarnAngle,ftol,abstol);
565 cmp_int(fp,"inputrec->nLincsIter",-1,ir1->nLincsIter,ir2->nLincsIter);
566 cmp_real(fp,"inputrec->bd_fric",-1,ir1->bd_fric,ir2->bd_fric,ftol,abstol);
567 cmp_int(fp,"inputrec->ld_seed",-1,ir1->ld_seed,ir2->ld_seed);
568 cmp_real(fp,"inputrec->cos_accel",-1,ir1->cos_accel,ir2->cos_accel,ftol,abstol);
569 cmp_rvec(fp,"inputrec->deform(a)",-1,ir1->deform[XX],ir2->deform[XX],ftol,abstol);
570 cmp_rvec(fp,"inputrec->deform(b)",-1,ir1->deform[YY],ir2->deform[YY],ftol,abstol);
571 cmp_rvec(fp,"inputrec->deform(c)",-1,ir1->deform[ZZ],ir2->deform[ZZ],ftol,abstol);
572
573 cmp_int(fp,"ir->adress_type" ,-1,ir1->adress_type,ir2->adress_type);
574 cmp_int(fp,"ir->badress_new_wf" ,-1,ir1->badress_new_wf,ir2->badress_new_wf);
575 cmp_int(fp,"ir->badress_chempot_dx" ,-1,ir1->badress_chempot_dx,ir2->badress_chempot_dx);
576 cmp_real(fp,"ir->adress_const_wf" ,-1,ir1->adress_const_wf,ir2->adress_const_wf,ftol,abstol);
577 cmp_real(fp,"ir->adress_ex_width" ,-1,ir1->adress_ex_width,ir2->adress_ex_width,ftol,abstol);
578 cmp_real(fp,"ir->adress_hy_width" ,-1,ir1->adress_hy_width,ir2->adress_hy_width,ftol,abstol);
579 cmp_int(fp,"ir->adress_icor" ,-1,ir1->adress_icor,ir2->adress_icor);
580 cmp_int(fp,"ir->adress_ivdw" ,-1,ir1->adress_ivdw,ir2->adress_ivdw);
581 cmp_int(fp,"ir->adress_site" ,-1,ir1->adress_site,ir2->adress_site);
582 cmp_rvec(fp,"ir->adress_refs" ,-1,ir1->adress_refs,ir2->adress_refs,ftol,abstol);
583
584 cmp_int(fp,"inputrec->userint1",-1,ir1->userint1,ir2->userint1);
585 cmp_int(fp,"inputrec->userint2",-1,ir1->userint2,ir2->userint2);
586 cmp_int(fp,"inputrec->userint3",-1,ir1->userint3,ir2->userint3);
587 cmp_int(fp,"inputrec->userint4",-1,ir1->userint4,ir2->userint4);
588 cmp_real(fp,"inputrec->userreal1",-1,ir1->userreal1,ir2->userreal1,ftol,abstol);
589 cmp_real(fp,"inputrec->userreal2",-1,ir1->userreal2,ir2->userreal2,ftol,abstol);
590 cmp_real(fp,"inputrec->userreal3",-1,ir1->userreal3,ir2->userreal3,ftol,abstol);
591 cmp_real(fp,"inputrec->userreal4",-1,ir1->userreal4,ir2->userreal4,ftol,abstol);
592 cmp_grpopts(fp,&(ir1->opts),&(ir2->opts),ftol,abstol);
593 cmp_cosines(fp,"ex",ir1->ex,ir2->ex,ftol,abstol);
594 cmp_cosines(fp,"et",ir1->et,ir2->et,ftol,abstol);
595 }
596
597 static void comp_pull_AB(FILE *fp,t_pull *pull,real ftol,real abstol)
598 {
599 int i;
600
601 for(i=0; i<pull->ngrp+1; i++) {
602 fprintf(fp,"comparing pull group %d\n",i);
603 cmp_real(fp,"pullgrp->k",-1,pull->grp[i].k,pull->grp[i].kB,ftol,abstol);
604 }
605 }
606
607 static void comp_state(t_state *st1, t_state *st2,
608 bool bRMSD,real ftol,real abstol)
609 {
610 int i,j,nc;
611
612 fprintf(stdout,"comparing flags\n");
613 cmp_int(stdout,"flags",-1,st1->flags,st2->flags);
614 fprintf(stdout,"comparing box\n");
615 cmp_rvecs(stdout,"box",DIM,st1->box,st2->box,FALSE,ftol,abstol);
616 fprintf(stdout,"comparing box_rel\n");
617 cmp_rvecs(stdout,"box_rel",DIM,st1->box_rel,st2->box_rel,FALSE,ftol,abstol);
618 fprintf(stdout,"comparing boxv\n");
619 cmp_rvecs(stdout,"boxv",DIM,st1->boxv,st2->boxv,FALSE,ftol,abstol);
620 if (st1->flags & (1<<estSVIR_PREV)) {
621 fprintf(stdout,"comparing shake vir_prev\n");
622 cmp_rvecs(stdout,"svir_prev",DIM,st1->svir_prev,st2->svir_prev,FALSE,ftol,abstol);
623 }
624 if (st1->flags & (1<<estFVIR_PREV)) {
625 fprintf(stdout,"comparing force vir_prev\n");
626 cmp_rvecs(stdout,"fvir_prev",DIM,st1->fvir_prev,st2->fvir_prev,FALSE,ftol,abstol);
627 }
628 if (st1->flags & (1<<estPRES_PREV)) {
629 fprintf(stdout,"comparing prev_pres\n");
630 cmp_rvecs(stdout,"pres_prev",DIM,st1->pres_prev,st2->pres_prev,FALSE,ftol,abstol);
631 }
632 cmp_int(stdout,"ngtc",-1,st1->ngtc,st2->ngtc);
633 cmp_int(stdout,"nhchainlength",-1,st1->nhchainlength,st2->nhchainlength);
634 if (st1->ngtc == st2->ngtc && st1->nhchainlength == st2->nhchainlength){
635 for(i=0; i<st1->ngtc; i++) {
636 nc = i*st1->nhchainlength;
637 for(j=0; j<nc; j++) {
638 cmp_real(stdout,"nosehoover_xi",
639 i,st1->nosehoover_xi[nc+j],st2->nosehoover_xi[nc+j],ftol,abstol);
640 }
641 }
642 }
643 cmp_int(stdout,"nnhpres",-1,st1->nnhpres,st2->nnhpres);
644 if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength) {
645 for(i=0; i<st1->nnhpres; i++) {
646 nc = i*st1->nhchainlength;
647 for(j=0; j<nc; j++) {
648 cmp_real(stdout,"nosehoover_xi",
649 i,st1->nhpres_xi[nc+j],st2->nhpres_xi[nc+j],ftol,abstol);
650 }
651 }
652 }
653
654 cmp_int(stdout,"natoms",-1,st1->natoms,st2->natoms);
655 if (st1->natoms == st2->natoms) {
656 if ((st1->flags & (1<<estX)) && (st2->flags & (1<<estX))) {
657 fprintf(stdout,"comparing x\n");
658 cmp_rvecs(stdout,"x",st1->natoms,st1->x,st2->x,bRMSD,ftol,abstol);
659 }
660 if ((st1->flags & (1<<estV)) && (st2->flags & (1<<estV))) {
661 fprintf(stdout,"comparing v\n");
662 cmp_rvecs(stdout,"v",st1->natoms,st1->v,st2->v,bRMSD,ftol,abstol);
663 }
664 }
665 }
666
667 void comp_tpx(const char *fn1,const char *fn2,
668 bool bRMSD,real ftol,real abstol)
669 {
670 const char *ff[2];
671 t_tpxheader sh[2];
672 t_inputrec ir[2];
673 t_state state[2];
674 gmx_mtop_t mtop[2];
675 t_topology top[2];
676 int i;
677
678 ff[0]=fn1;
679 ff[1]=fn2;
680 for(i=0; i<(fn2 ? 2 : 1); i++) {
681 read_tpx_state(ff[i],&(ir[i]),&state[i],NULL,&(mtop[i]));
682 }
683 if (fn2) {
684 cmp_inputrec(stdout,&ir[0],&ir[1],ftol,abstol);
685 /* Convert gmx_mtop_t to t_topology.
686 * We should implement direct mtop comparison,
687 * but it might be useful to keep t_topology comparison as an option.
688 */
689 top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
690 top[1] = gmx_mtop_t_to_t_topology(&mtop[1]);
691 cmp_top(stdout,&top[0],&top[1],ftol,abstol);
692 comp_state(&state[0],&state[1],bRMSD,ftol,abstol);
693 } else {
694 if (ir[0].efep == efepNO) {
695 fprintf(stdout,"inputrec->efep = %s\n",efep_names[ir[0].efep]);
696 } else {
697 if (ir[0].ePull != epullNO) {
698 comp_pull_AB(stdout,ir->pull,ftol,abstol);
699 }
700 /* Convert gmx_mtop_t to t_topology.
701 * We should implement direct mtop comparison,
702 * but it might be useful to keep t_topology comparison as an option.
703 */
704 top[0] = gmx_mtop_t_to_t_topology(&mtop[0]);
705 cmp_top(stdout,&top[0],NULL,ftol,abstol);
706 }
707 }
708 }
709
710 void comp_frame(FILE *fp, t_trxframe *fr1, t_trxframe *fr2,
711 bool bRMSD, real ftol,real abstol)
712 {
713 fprintf(fp,"\n");
714 cmp_int(fp,"flags",-1,fr1->flags,fr2->flags);
715 cmp_int(fp,"not_ok",-1,fr1->not_ok,fr2->not_ok);
716 cmp_int(fp,"natoms",-1,fr1->natoms,fr2->natoms);
717 cmp_real(fp,"t0",-1,fr1->t0,fr2->t0,ftol,abstol);
718 if (cmp_bool(fp,"bTitle",-1,fr1->bTitle,fr2->bTitle))
719 cmp_str(fp,"title", -1, fr1->title, fr2->title);
720 if (cmp_bool(fp,"bStep",-1,fr1->bStep,fr2->bStep))
721 cmp_int(fp,"step",-1,fr1->step,fr2->step);
722 cmp_int(fp,"step",-1,fr1->step,fr2->step);
723 if (cmp_bool(fp,"bTime",-1,fr1->bTime,fr2->bTime))
724 cmp_real(fp,"time",-1,fr1->time,fr2->time,ftol,abstol);
725 if (cmp_bool(fp,"bLambda",-1,fr1->bLambda,fr2->bLambda))
726 cmp_real(fp,"lambda",-1,fr1->lambda,fr2->lambda,ftol,abstol);
727 if (cmp_bool(fp,"bAtoms",-1,fr1->bAtoms,fr2->bAtoms))
728 cmp_atoms(fp,fr1->atoms,fr2->atoms,ftol,abstol);
729 if (cmp_bool(fp,"bPrec",-1,fr1->bPrec,fr2->bPrec))
730 cmp_real(fp,"prec",-1,fr1->prec,fr2->prec,ftol,abstol);
731 if (cmp_bool(fp,"bX",-1,fr1->bX,fr2->bX))
732 cmp_rvecs(fp,"x",min(fr1->natoms,fr2->natoms),fr1->x,fr2->x,bRMSD,ftol,abstol);
733 if (cmp_bool(fp,"bV",-1,fr1->bV,fr2->bV))
734 cmp_rvecs(fp,"v",min(fr1->natoms,fr2->natoms),fr1->v,fr2->v,bRMSD,ftol,abstol);
735 if (cmp_bool(fp,"bF",-1,fr1->bF,fr2->bF))
736 cmp_rvecs(fp,"f",min(fr1->natoms,fr2->natoms),fr1->f,fr2->f,bRMSD,ftol,abstol);
737 if (cmp_bool(fp,"bBox",-1,fr1->bBox,fr2->bBox))
738 cmp_rvecs(fp,"box",3,fr1->box,fr2->box,FALSE,ftol,abstol);
739 }
740
741 void comp_trx(const output_env_t oenv,const char *fn1, const char *fn2,
742 bool bRMSD,real ftol,real abstol)
743 {
744 int i;
745 const char *fn[2];
746 t_trxframe fr[2];
747 int status[2];
748 bool b[2];
749
750 fn[0]=fn1;
751 fn[1]=fn2;
752 fprintf(stderr,"Comparing trajectory files %s and %s\n",fn1,fn2);
753 for (i=0; i<2; i++)
754 b[i] = read_first_frame(oenv,&status[i],fn[i],&fr[i],TRX_READ_X|TRX_READ_V|TRX_READ_F);
755
756 if (b[0] && b[1]) {
757 do {
758 comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
759
760 for (i=0; i<2; i++)
761 b[i] = read_next_frame(oenv,status[i],&fr[i]);
762 } while (b[0] && b[1]);
763
764 for (i=0; i<2; i++) {
765 if (b[i] && !b[1-i])
766 fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn[i],fn[1-i]);
767 close_trj(status[i]);
768 }
769 }
770 if (!b[0] && !b[1])
771 fprintf(stdout,"\nBoth files read correctly\n");
772 }
773
774 static void cmp_energies(FILE *fp,int step1,int step2,int nre,
775 t_energy e1[],t_energy e2[],
776 gmx_enxnm_t *enm1,gmx_enxnm_t *enm2,
777 real ftol,real abstol,
778 int maxener)
779 {
780 int i;
781
782 for(i=0; (i<maxener); i++) {
783 if (!equal_real(e1[i].e,e2[i].e,ftol,abstol))
784 fprintf(fp,"%-15s step %3d: %12g, %s step %3d: %12g\n",
785 enm1[i].name,step1,e1[i].e,
786 strcmp(enm1[i].name,enm2[i].name)!=0 ? enm2[i].name:"",
787 step2,e2[i].e);
788 }
789 }
790
791 static void cmp_disres(t_enxframe *fr1,t_enxframe *fr2,real ftol, real abstol)
792 {
793 int i;
794 char bav[64],bt[64],bs[22];
795
796 cmp_int(stdout,"ndisre",-1,fr1->ndisre,fr2->ndisre);
797 if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0)) {
798 sprintf(bav,"step %s: disre rav",gmx_step_str(fr1->step,bs));
799 sprintf(bt, "step %s: disre rt",gmx_step_str(fr1->step,bs));
800 for(i=0; (i<fr1->ndisre); i++) {
801 cmp_real(stdout,bav,i,fr1->disre_rm3tav[i],fr2->disre_rm3tav[i],ftol,abstol);
802 cmp_real(stdout,bt ,i,fr1->disre_rt[i] ,fr2->disre_rt[i] ,ftol,abstol);
803 }
804 }
805 }
806
807 static void cmp_eblocks(t_enxframe *fr1,t_enxframe *fr2,real ftol, real abstol)
808 {
809 int i,j;
810 char buf[64],bs[22];
811
812 cmp_int(stdout,"nblock",-1,fr1->nblock,fr2->nblock);
813 if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0)) {
814 for(j=0; (j<fr1->nblock); j++) {
815 sprintf(buf,"step %s: block[%d]",gmx_step_str(fr1->step,bs),j);
816 cmp_int(stdout,buf,-1,fr1->nr[j],fr2->nr[j]);
817 if ((fr1->nr[j] == fr2->nr[j]) && (fr1->nr[j] > 0)) {
818 for(i=0; (i<fr1->nr[j]); i++) {
819 cmp_real(stdout,buf,i,fr1->block[i][j],fr2->block[i][j],ftol,abstol);
820 }
821 }
822 }
823 }
824 }
825
826 void comp_enx(const char *fn1,const char *fn2,real ftol,real abstol,const char *lastener)
827 {
828 int nre,nre1,nre2,block;
829 ener_file_t in1, in2;
830 int i,maxener;
831 char buf[256];
832 gmx_enxnm_t *enm1=NULL,*enm2=NULL;
833 t_enxframe *fr1,*fr2;
834 bool b1,b2;
835
836 fprintf(stdout,"comparing energy file %s and %s\n\n",fn1,fn2);
837
838 in1 = open_enx(fn1,"r");
839 in2 = open_enx(fn2,"r");
840 do_enxnms(in1,&nre1,&enm1);
841 do_enxnms(in2,&nre2,&enm2);
842 if (nre1 != nre2) {
843 fprintf(stdout,"%s: nre=%d, %s: nre=%d\n",fn1,nre1,fn2,nre2);
844 return;
845 }
846 nre = nre1;
847 fprintf(stdout,"There are %d terms in the energy files\n\n",nre);
848
849 maxener = nre;
850 for(i=0; (i<nre); i++) {
851 cmp_str(stdout,"enm",i,enm1[i].name,enm2[i].name);
852 cmp_str(stdout,"unit",i,enm1[i].unit,enm2[i].unit);
853 if ((lastener != NULL) && (strstr(enm1[i].name,lastener) != NULL)) {
854 maxener=i+1;
855 break;
856 }
857 }
858 fprintf(stdout,"There are %d terms to compare in the energy files\n\n",
859 maxener);
860
861 snew(fr1,1);
862 snew(fr2,1);
863 do {
864 b1 = do_enx(in1,fr1);
865 b2 = do_enx(in2,fr2);
866 if (b1 && !b2)
867 fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn2,fn1);
868 else if (!b1 && b2)
869 fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn1,fn2);
870 else if (!b1 && !b2)
871 fprintf(stdout,"\nFiles read succesfully\n");
872 else {
873 cmp_real(stdout,"t",-1,fr1->t,fr2->t,ftol,abstol);
874 cmp_int(stdout,"step",-1,fr1->step,fr2->step);
875 cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre);
876 if ((fr1->nre == nre) && (fr2->nre == nre))
877 cmp_energies(stdout,fr1->step,fr1->step,nre,fr1->ener,fr2->ener,
878 enm1,enm2,ftol,abstol,maxener);
879 cmp_disres(fr1,fr2,ftol,abstol);
880 cmp_eblocks(fr1,fr2,ftol,abstol);
881 }
882 } while (b1 && b2);
883
884 close_enx(in1);
885 close_enx(in2);
886
887 free_enxframe(fr2);
888 sfree(fr2);
889 free_enxframe(fr1);
890 sfree(fr1);
891 }
got merge a long time ago