mshift.c
1 
/*


2 
* $Id: mshift.c,v 1.62.2.3 2008/12/11 14:24:43 hess Exp $

3 
*

4 
* This source code is part of

5 
*

6 
* G R O M A C S

7 
*

8 
* GROningen MAchine for Chemical Simulations

9 
*

10 
* VERSION 3.2.0

11 
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.

12 
* Copyright (c) 19912000, University of Groningen, The Netherlands.

13 
* Copyright (c) 20012004, The GROMACS development team,

14 
* check out http://www.gromacs.org for more information.

15 

16 
* This program is free software; you can redistribute it and/or

17 
* modify it under the terms of the GNU General Public License

18 
* as published by the Free Software Foundation; either version 2

19 
* of the License, or (at your option) any later version.

20 
*

21 
* If you want to redistribute modifications, please consider that

22 
* scientific software is very special. Version control is crucial 

23 
* bugs must be traceable. We will be happy to consider code for

24 
* inclusion in the official distribution, but derived work must not

25 
* be called official GROMACS. Details are found in the README & COPYING

26 
* files  if they are missing, get the official version at www.gromacs.org.

27 
*

28 
* To help us fund GROMACS development, we humbly ask that you cite

29 
* the papers on the package  you can find them in the top README file.

30 
*

31 
* For more info, check our website at http://www.gromacs.org

32 
*

33 
* And Hey:

34 
* GROningen Mixture of Alchemy and Childrens' Stories

35 
*/

36 
#ifdef HAVE_CONFIG_H

37 
#include <config.h> 
38 
#endif

39  
40 
#include <string.h> 
41 
#include "smalloc.h" 
42 
#include "gmx_fatal.h" 
43 
#include "macros.h" 
44 
#include "vec.h" 
45 
#include "futil.h" 
46 
#include "copyrite.h" 
47 
#include "mshift.h" 
48 
#include "main.h" 
49 
#include "pbc.h" 
50  
51 
/************************************************************

52 
*

53 
* S H I F T U T I L I T I E S

54 
*

55 
************************************************************/

56 

57  
58 
/************************************************************

59 
*

60 
* G R A P H G E N E R A T I O N C O D E

61 
*

62 
************************************************************/

63  
64 
static void add_gbond(t_graph *g,atom_id a0,atom_id a1) 
65 
{ 
66 
int i;

67 
atom_id inda0,inda1; 
68 
bool bFound;

69  
70 
inda0 = a0  g>start; 
71 
inda1 = a1  g>start; 
72 
bFound = FALSE; 
73 
/* Search for a direct edge between a0 and a1.

74 
* All egdes are bidirectional, so we only need to search one way.

75 
*/

76 
for(i=0; (i<g>nedge[inda0] && !bFound); i++) { 
77 
bFound = (g>edge[inda0][i] == a1); 
78 
} 
79  
80 
if (!bFound) {

81 
g>edge[inda0][g>nedge[inda0]++] = a1; 
82 
g>edge[inda1][g>nedge[inda1]++] = a0; 
83 
} 
84 
} 
85  
86 
static void mk_igraph(t_graph *g,int ftype,t_ilist *il, 
87 
int at_start,int at_end, 
88 
int *part)

89 
{ 
90 
t_iatom *ia; 
91 
int i,j,np;

92 
int end;

93  
94 
end=il>nr; 
95 
ia=il>iatoms; 
96  
97 
i = 0;

98 
while (i < end) {

99 
np = interaction_function[ftype].nratoms; 
100 

101 
if (ia[1] >= at_start && ia[1] < at_end) { 
102 
if (ia[np] >= at_end)

103 
gmx_fatal(FARGS, 
104 
"Molecule in topology has atom numbers below and "

105 
"above natoms (%d).\n"

106 
"You are probably trying to use a trajectory which does "

107 
"not match the first %d atoms of the run input file.\n"

108 
"You can make a matching run input file with tpbconv.",

109 
at_end,at_end); 
110 
if (ftype == F_SETTLE) {

111 
/* Bond all the atoms in the settle */

112 
add_gbond(g,ia[1],ia[1]+1); 
113 
add_gbond(g,ia[1],ia[1]+2); 
114 
} else if (part == NULL) { 
115 
/* Simply add this bond */

116 
for(j=1; j<np; j++) { 
117 
add_gbond(g,ia[j],ia[j+1]);

118 
} 
119 
} else {

120 
/* Add this bond when it connects two unlinked parts of the graph */

121 
for(j=1; j<np; j++) { 
122 
if (part[ia[j]] != part[ia[j+1]]) { 
123 
add_gbond(g,ia[j],ia[j+1]);

124 
} 
125 
} 
126 
} 
127 
} 
128 
ia+=np+1;

129 
i+=np+1;

130 
} 
131 
} 
132  
133 
static void g_error(int line,char *file) 
134 
{ 
135 
gmx_fatal(FARGS,"Tring to print non existant graph (file %s, line %d)",

136 
file,line); 
137 
} 
138 
#define GCHECK(g) if (g == NULL) g_error(__LINE__,__FILE__) 
139  
140 
void p_graph(FILE *log,char *title,t_graph *g) 
141 
{ 
142 
int i,j;

143 
char *cc[egcolNR] = { "W", "G", "B" }; 
144 

145 
GCHECK(g); 
146 
fprintf(log,"graph: %s\n",title);

147 
fprintf(log,"nnodes: %d\n",g>nnodes);

148 
fprintf(log,"nbound: %d\n",g>nbound);

149 
fprintf(log,"start: %d\n",g>start);

150 
fprintf(log,"end: %d\n",g>end);

151 
fprintf(log," atom shiftx shifty shiftz C nedg e1 e2 etc.\n");

152 
for(i=0; (i<g>nnodes); i++) 
153 
if (g>nedge[i] > 0) { 
154 
fprintf(log,"%5d%7d%7d%7d %1s%5d",g>start+i+1, 
155 
g>ishift[i][XX],g>ishift[i][YY], 
156 
g>ishift[i][ZZ], 
157 
(g>negc > 0) ? cc[g>egc[i]] : " ", 
158 
g>nedge[i]); 
159 
for(j=0; (j<g>nedge[i]); j++) 
160 
fprintf(log," %5u",g>edge[i][j]+1); 
161 
fprintf(log,"\n");

162 
} 
163 
fflush(log); 
164 
} 
165  
166 
static void calc_1se(t_graph *g,int ftype,t_ilist *il, 
167 
int nbond[],int at_start,int at_end) 
168 
{ 
169 
int k,nratoms,end,j;

170 
t_iatom *ia,iaa; 
171  
172 
end=il>nr; 
173  
174 
ia=il>iatoms; 
175 
for(j=0; (j<end); j+=nratoms+1,ia+=nratoms+1) { 
176 
nratoms = interaction_function[ftype].nratoms; 
177 

178 
if (ftype == F_SETTLE) {

179 
iaa = ia[1];

180 
if (iaa >= at_start && iaa < at_end) {

181 
nbond[iaa] += 2;

182 
nbond[iaa+1] += 1; 
183 
nbond[iaa+2] += 1; 
184 
g>start = min(g>start,iaa); 
185 
g>end = max(g>end,iaa+2);

186 
} 
187 
} else {

188 
for(k=1; (k<=nratoms); k++) { 
189 
iaa=ia[k]; 
190 
if (iaa >= at_start && iaa < at_end) {

191 
g>start=min(g>start,iaa); 
192 
g>end =max(g>end, iaa); 
193 
/*if (interaction_function[tp].flags & IF_CHEMBOND)*/

194 
nbond[iaa]++; 
195 
} 
196 
} 
197 
} 
198 
} 
199 
} 
200  
201 
static int calc_start_end(FILE *fplog,t_graph *g,t_ilist il[], 
202 
int at_start,int at_end, 
203 
int nbond[])

204 
{ 
205 
int i,nnb,nbtot;

206 

207 
g>start=at_end; 
208 
g>end=0;

209  
210 
/* First add all the real bonds: they should determine the molecular

211 
* graph.

212 
*/

213 
for(i=0; (i<F_NRE); i++) 
214 
if (interaction_function[i].flags & IF_CHEMBOND)

215 
calc_1se(g,i,&il[i],nbond,at_start,at_end); 
216 
/* Then add all the other interactions in fixed lists, but first

217 
* check to see what's there already.

218 
*/

219 
for(i=0; (i<F_NRE); i++) { 
220 
if (!(interaction_function[i].flags & IF_CHEMBOND)) {

221 
calc_1se(g,i,&il[i],nbond,at_start,at_end); 
222 
} 
223 
} 
224 

225 
nnb = 0;

226 
nbtot = 0;

227 
for(i=g>start; (i<=g>end); i++) {

228 
nbtot += nbond[i]; 
229 
nnb = max(nnb,nbond[i]); 
230 
} 
231 
if (fplog) {

232 
fprintf(fplog,"Max number of connections per atom is %d\n",nnb);

233 
fprintf(fplog,"Total number of connections is %d\n",nbtot);

234 
} 
235 
return nbtot;

236 
} 
237  
238  
239  
240 
static void compact_graph(FILE *fplog,t_graph *g) 
241 
{ 
242 
int i,j,n,max_nedge;

243 
atom_id *e; 
244  
245 
max_nedge = 0;

246 
n = 0;

247 
for(i=0; i<g>nnodes; i++) { 
248 
srenew(g>edge[i],g>nedge[i]); 
249 
max_nedge = max(max_nedge,g>nedge[i]); 
250 
n+=g>nedge[i]; 
251 
} 
252 
if (fplog) {

253 
fprintf(fplog,"Max number of graph edges per atom is %d\n",

254 
max_nedge); 
255 
fprintf(fplog,"Total number of graph edges is %d\n",n);

256 
} 
257 
} 
258  
259 
static bool determine_graph_parts(t_graph *g,int *part) 
260 
{ 
261 
int i,e;

262 
int nchanged;

263 
atom_id at_i,*at_i2; 
264 
bool bMultiPart;

265  
266 
/* Initialize the part array with all entries different */

267 
for(at_i=g>start; at_i<g>end; at_i++) {

268 
part[at_i] = at_i; 
269 
} 
270  
271 
/* Loop over the graph until the part array is fixed */

272 
do {

273 
bMultiPart = FALSE; 
274 
nchanged = 0;

275 
for(i=0; (i<g>nnodes); i++) { 
276 
at_i = g>start + i; 
277 
at_i2 = g>edge[i]; 
278 
for(e=0; e<g>nedge[i]; e++) { 
279 
/* Set part for both nodes to the minimum */

280 
if (part[at_i2[e]] > part[at_i]) {

281 
part[at_i2[e]] = part[at_i]; 
282 
nchanged++; 
283 
} else if (part[at_i2[e]] < part[at_i]) { 
284 
part[at_i] = part[at_i2[e]]; 
285 
nchanged++; 
286 
} 
287 
} 
288 
if (part[at_i] != part[g>start]) {

289 
bMultiPart = TRUE; 
290 
} 
291 
} 
292 
if (debug) {

293 
fprintf(debug,"graph part[] nchanged=%d, bMultiPart=%d\n",

294 
nchanged,bMultiPart); 
295 
} 
296 
} while (nchanged > 0); 
297  
298 
return bMultiPart;

299 
} 
300  
301 
void mk_graph_ilist(FILE *fplog,

302 
t_ilist *ilist,int at_start,int at_end, 
303 
bool bShakeOnly,bool bSettle, 
304 
t_graph *g) 
305 
{ 
306 
int *nbond;

307 
int i,nbtot;

308 
bool bMultiPart;

309  
310 
snew(nbond,at_end); 
311 
nbtot = calc_start_end(fplog,g,ilist,at_start,at_end,nbond); 
312 

313 
if (g>start >= g>end) {

314 
g>nnodes = 0;

315 
g>nbound = 0;

316 
} 
317 
else {

318 
g>nnodes = g>end  g>start + 1;

319 
snew(g>ishift,g>nnodes); 
320 
snew(g>nedge,g>nnodes); 
321 
snew(g>edge,g>nnodes); 
322 
for(i=0; (i<g>nnodes); i++) 
323 
snew(g>edge[i],nbond[g>start+i]); 
324  
325 
if (!bShakeOnly) {

326 
/* First add all the real bonds: they should determine the molecular

327 
* graph.

328 
*/

329 
for(i=0; (i<F_NRE); i++) 
330 
if (interaction_function[i].flags & IF_CHEMBOND)

331 
mk_igraph(g,i,&(ilist[i]),at_start,at_end,NULL);

332  
333 
/* Determine of which separated parts the IF_CHEMBOND graph consists.

334 
* Store the parts in the nbond array.

335 
*/

336 
bMultiPart = determine_graph_parts(g,nbond); 
337  
338 
if (bMultiPart) {

339 
/* Then add all the other interactions in fixed lists,

340 
* but only when they connect parts of the graph

341 
* that are not connected through IF_CHEMBOND interactions.

342 
*/

343 
for(i=0; (i<F_NRE); i++) { 
344 
if (!(interaction_function[i].flags & IF_CHEMBOND)) {

345 
mk_igraph(g,i,&(ilist[i]),at_start,at_end,nbond); 
346 
} 
347 
} 
348 
} 
349 

350 
/* Removed all the unused space from the edge array */

351 
compact_graph(fplog,g); 
352 
} 
353 
else {

354 
/* This is a special thing used in splitter.c to generate shakeblocks */

355 
mk_igraph(g,F_CONSTR,&(ilist[F_CONSTR]),at_start,at_end,NULL);

356 
if (bSettle)

357 
mk_igraph(g,F_SETTLE,&(ilist[F_SETTLE]),at_start,at_end,NULL);

358 
} 
359 
g>nbound = 0;

360 
for(i=0; (i<g>nnodes); i++) 
361 
if (g>nedge[i] > 0) 
362 
g>nbound++; 
363 
} 
364  
365 
g>negc = 0;

366 
g>egc = NULL;

367  
368 
sfree(nbond); 
369  
370 
if (gmx_debug_at)

371 
p_graph(debug,"graph",g);

372 
} 
373  
374 
t_graph *mk_graph(FILE *fplog, 
375 
t_idef *idef,int at_start,int at_end, 
376 
bool bShakeOnly,bool bSettle) 
377 
{ 
378 
t_graph *g; 
379  
380 
snew(g,1);

381  
382 
mk_graph_ilist(fplog,idef>il,at_start,at_end,bShakeOnly,bSettle,g); 
383  
384 
return g;

385 
} 
386  
387 
void done_graph(t_graph *g)

388 
{ 
389 
int i;

390 

391 
GCHECK(g); 
392 
if (g>nnodes > 0) { 
393 
sfree(g>ishift); 
394 
sfree(g>nedge); 
395 
for(i=0; (i<g>nnodes); i++) 
396 
sfree(g>edge[i]); 
397 
sfree(g>edge); 
398 
sfree(g>egc); 
399 
} 
400 
} 
401  
402 
/************************************************************

403 
*

404 
* S H I F T C A L C U L A T I O N C O D E

405 
*

406 
************************************************************/

407  
408 
static void mk_1shift_tric(int npbcdim,matrix box,rvec hbox, 
409 
rvec xi,rvec xj,int *mi,int *mj) 
410 
{ 
411 
/* Calculate periodicity for triclinic box... */

412 
int m,d;

413 
rvec dx; 
414 

415 
rvec_sub(xi,xj,dx); 
416  
417 
mj[ZZ] = 0;

418 
for(m=npbcdim1; (m>=0); m) { 
419 
/* If dx < hbox, then xj will be reduced by box, so that

420 
* xi  xj will be bigger

421 
*/

422 
if (dx[m] < hbox[m]) {

423 
mj[m]=mi[m]1;

424 
for(d=m1; d>=0; d) 
425 
dx[d]+=box[m][d]; 
426 
} else if (dx[m] >= hbox[m]) { 
427 
mj[m]=mi[m]+1;

428 
for(d=m1; d>=0; d) 
429 
dx[d]=box[m][d]; 
430 
} else

431 
mj[m]=mi[m]; 
432 
} 
433 
} 
434  
435 
static void mk_1shift(int npbcdim,rvec hbox,rvec xi,rvec xj,int *mi,int *mj) 
436 
{ 
437 
/* Calculate periodicity for rectangular box... */

438 
int m;

439 
rvec dx; 
440 

441 
rvec_sub(xi,xj,dx); 
442  
443 
mj[ZZ] = 0;

444 
for(m=0; (m<npbcdim); m++) { 
445 
/* If dx < hbox, then xj will be reduced by box, so that

446 
* xi  xj will be bigger

447 
*/

448 
if (dx[m] < hbox[m])

449 
mj[m]=mi[m]1;

450 
else if (dx[m] >= hbox[m]) 
451 
mj[m]=mi[m]+1;

452 
else

453 
mj[m]=mi[m]; 
454 
} 
455 
} 
456  
457 
static void mk_1shift_screw(matrix box,rvec hbox, 
458 
rvec xi,rvec xj,int *mi,int *mj) 
459 
{ 
460 
/* Calculate periodicity for rectangular box... */

461 
int signi,m;

462 
rvec dx; 
463  
464 
if ((mi[XX] > 0 && mi[XX] % 2 == 1)  
465 
(mi[XX] < 0 && mi[XX] % 2 == 1)) { 
466 
signi = 1;

467 
} else {

468 
signi = 1;

469 
} 
470  
471 
rvec_sub(xi,xj,dx); 
472  
473 
if (dx[XX] < hbox[XX])

474 
mj[XX] = mi[XX]  1;

475 
else if (dx[XX] >= hbox[XX]) 
476 
mj[XX] = mi[XX] + 1;

477 
else

478 
mj[XX] = mi[XX]; 
479 
if (mj[XX] != mi[XX]) {

480 
/* Rotate */

481 
dx[YY] = xi[YY]  (box[YY][YY] + box[ZZ][YY]  xj[YY]); 
482 
dx[ZZ] = xi[ZZ]  (box[ZZ][ZZ]  xj[ZZ]); 
483 
} 
484 
for(m=1; (m<DIM); m++) { 
485 
/* The signs are taken such that we can first shift x and rotate

486 
* and then shift y and z.

487 
*/

488 
if (dx[m] < hbox[m])

489 
mj[m] = mi[m]  signi; 
490 
else if (dx[m] >= hbox[m]) 
491 
mj[m] = mi[m] + signi; 
492 
else

493 
mj[m] = mi[m]; 
494 
} 
495 
} 
496  
497 
static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI, 
498 
int npbcdim,matrix box,rvec x[],int *nerror) 
499 
{ 
500 
int m,j,ng,ai,aj,g0;

501 
rvec dx,hbox; 
502 
bool bTriclinic;

503 
ivec is_aj; 
504 
t_pbc pbc; 
505 

506 
for(m=0; (m<DIM); m++) 
507 
hbox[m]=box[m][m]*0.5; 
508 
bTriclinic = TRICLINIC(box); 
509 

510 
ng=0;

511 
ai=*AtomI; 
512 

513 
g0=g>start; 
514 
/* Loop over all the bonds */

515 
for(j=0; (j<g>nedge[ai]); j++) { 
516 
aj=g>edge[ai][j]g0; 
517 
/* If there is a white one, make it grey and set pbc */

518 
if (g>bScrewPBC)

519 
mk_1shift_screw(box,hbox,x[g0+ai],x[g0+aj],g>ishift[ai],is_aj); 
520 
else if (bTriclinic) 
521 
mk_1shift_tric(npbcdim,box,hbox,x[g0+ai],x[g0+aj],g>ishift[ai],is_aj); 
522 
else

523 
mk_1shift(npbcdim,hbox,x[g0+ai],x[g0+aj],g>ishift[ai],is_aj); 
524 

525 
if (egc[aj] == egcolWhite) {

526 
if (aj < *AtomI)

527 
*AtomI = aj; 
528 
egc[aj] = egcolGrey; 
529 

530 
copy_ivec(is_aj,g>ishift[aj]); 
531  
532 
ng++; 
533 
} 
534 
else if ((is_aj[XX] != g>ishift[aj][XX])  
535 
(is_aj[YY] != g>ishift[aj][YY])  
536 
(is_aj[ZZ] != g>ishift[aj][ZZ])) { 
537 
if (gmx_debug_at) {

538 
set_pbc(&pbc,1,box);

539 
pbc_dx(&pbc,x[g0+ai],x[g0+aj],dx); 
540 
fprintf(debug,"mk_grey: shifts for atom %d due to atom %d\n"

541 
"are (%d,%d,%d), should be (%d,%d,%d)\n"

542 
"dx = (%g,%g,%g)\n",

543 
aj+g0+1,ai+g0+1,is_aj[XX],is_aj[YY],is_aj[ZZ], 
544 
g>ishift[aj][XX],g>ishift[aj][YY],g>ishift[aj][ZZ], 
545 
dx[XX],dx[YY],dx[ZZ]); 
546 
} 
547 
(*nerror)++; 
548 
} 
549 
} 
550 
return ng;

551 
} 
552  
553 
static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[]) 
554 
/* Return the first node with colour Col starting at fC.

555 
* return 1 if none found.

556 
*/

557 
{ 
558 
int i;

559 

560 
for(i=fC; (i<g>nnodes); i++)

561 
if ((g>nedge[i] > 0) && (egc[i]==Col)) 
562 
return i;

563 

564 
return 1; 
565 
} 
566  
567 
void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[]) 
568 
{ 
569 
static int nerror_tot = 0; 
570 
int npbcdim;

571 
int ng,nnodes,i;

572 
int nW,nG,nB; /* Number of Grey, Black, White */ 
573 
int fW,fG; /* First of each category */ 
574 
int nerror=0; 
575  
576 
g>bScrewPBC = (ePBC == epbcSCREW); 
577  
578 
if (ePBC == epbcXY)

579 
npbcdim = 2;

580 
else

581 
npbcdim = 3;

582  
583 
GCHECK(g); 
584 
/* This puts everything in the central box, that is does not move it

585 
* at all. If we return without doing this for a system without bonds

586 
* (i.e. only settles) all water molecules are moved to the opposite octant

587 
*/

588 
for(i=0; (i<g>nnodes); i++) { 
589 
g>ishift[i][XX]=g>ishift[i][YY]=g>ishift[i][ZZ]=0;

590 
} 
591 

592 
if (!g>nbound)

593 
return;

594  
595 
nnodes=g>nnodes; 
596 
if (nnodes > g>negc) {

597 
g>negc = nnodes; 
598 
srenew(g>egc,g>negc); 
599 
} 
600 
memset(g>egc,0,(size_t)(nnodes*sizeof(g>egc[0]))); 
601  
602 
nW=g>nbound; 
603 
nG=0;

604 
nB=0;

605  
606 
fW=0;

607  
608 
/* We even have a loop invariant:

609 
* nW+nG+nB == g>nbound

610 
*/

611 
#ifdef DEBUG2

612 
fprintf(log,"Starting W loop\n");

613 
#endif

614 
while (nW > 0) { 
615 
/* Find the first white, this will allways be a larger

616 
* number than before, because no nodes are made white

617 
* in the loop

618 
*/

619 
if ((fW=first_colour(fW,egcolWhite,g,g>egc)) == 1) 
620 
gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW);

621 

622 
/* Make the first white node grey */

623 
g>egc[fW]=egcolGrey; 
624 
nG++; 
625 
nW; 
626  
627 
/* Initial value for the first grey */

628 
fG=fW; 
629 
#ifdef DEBUG2

630 
fprintf(log,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",

631 
nW,nG,nB,nW+nG+nB); 
632 
#endif

633 
while (nG > 0) { 
634 
if ((fG=first_colour(fG,egcolGrey,g,g>egc)) == 1) 
635 
gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG);

636 

637 
/* Make the first grey node black */

638 
g>egc[fG]=egcolBlack; 
639 
nB++; 
640 
nG; 
641  
642 
/* Make all the neighbours of this black node grey

643 
* and set their periodicity

644 
*/

645 
ng=mk_grey(log,nnodes,g>egc,g,&fG,npbcdim,box,x,&nerror); 
646 
/* ng is the number of white nodes made grey */

647 
nG+=ng; 
648 
nW=ng; 
649 
} 
650 
} 
651 
if (nerror > 0) { 
652 
nerror_tot++; 
653 
if (nerror_tot <= 100) { 
654 
fprintf(stderr,"There were %d inconsistent shifts. Check your topology\n",

655 
nerror); 
656 
if (log) {

657 
fprintf(log,"There were %d inconsistent shifts. Check your topology\n",

658 
nerror); 
659 
} 
660 
} 
661 
if (nerror_tot == 100) { 
662 
fprintf(stderr,"Will stop reporting inconsistent shifts\n");

663 
if (log) {

664 
fprintf(log,"Will stop reporting inconsistent shifts\n");

665 
} 
666 
} 
667 
} 
668 
} 
669  
670 
/************************************************************

671 
*

672 
* A C T U A L S H I F T C O D E

673 
*

674 
************************************************************/

675  
676 
void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])

677 
{ 
678 
ivec *is; 
679 
int g0,gn;

680 
int i,j,tx,ty,tz;

681  
682 
GCHECK(g); 
683 
g0=g>start; 
684 
gn=g>nnodes; 
685 
is=g>ishift; 
686 

687 
if (g>bScrewPBC) {

688 
for(i=0,j=g0; (i<gn); i++,j++) { 
689 
tx=is[i][XX]; 
690 
ty=is[i][YY]; 
691 
tz=is[i][ZZ]; 
692 

693 
if ((tx > 0 && tx % 2 == 1)  
694 
(tx < 0 && tx %2 == 1)) { 
695 
x_s[j][XX] = x[j][XX] + tx*box[XX][XX]; 
696 
x_s[j][YY] = box[YY][YY] + box[ZZ][YY]  x[j][YY]; 
697 
x_s[j][ZZ] = box[ZZ][ZZ]  x[j][ZZ]; 
698 
} else {

699 
x_s[j][XX] = x[j][XX]; 
700 
} 
701 
x_s[j][YY] = x[j][YY] + ty*box[YY][YY] + tz*box[ZZ][YY]; 
702 
x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ]; 
703 
} 
704 
} else if (TRICLINIC(box)) { 
705 
for(i=0,j=g0; (i<gn); i++,j++) { 
706 
tx=is[i][XX]; 
707 
ty=is[i][YY]; 
708 
tz=is[i][ZZ]; 
709 

710 
x_s[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX]; 
711 
x_s[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY]; 
712 
x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ]; 
713 
} 
714 
} else {

715 
for(i=0,j=g0; (i<gn); i++,j++) { 
716 
tx=is[i][XX]; 
717 
ty=is[i][YY]; 
718 
tz=is[i][ZZ]; 
719 

720 
x_s[j][XX]=x[j][XX]+tx*box[XX][XX]; 
721 
x_s[j][YY]=x[j][YY]+ty*box[YY][YY]; 
722 
x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ]; 
723 
} 
724 
} 
725 

726 
} 
727  
728 
void shift_self(t_graph *g,matrix box,rvec x[])

729 
{ 
730 
ivec *is; 
731 
int g0,gn;

732 
int i,j,tx,ty,tz;

733  
734 
if (g>bScrewPBC)

735 
gmx_incons("screw pbc not implemented for shift_self");

736  
737 
g0=g>start; 
738 
gn=g>nnodes; 
739 
is=g>ishift; 
740  
741 
#ifdef DEBUG

742 
fprintf(stderr,"Shifting atoms %d to %d\n",g0,g0+gn);

743 
#endif

744 
if(TRICLINIC(box)) {

745 
for(i=0,j=g0; (i<gn); i++,j++) { 
746 
tx=is[i][XX]; 
747 
ty=is[i][YY]; 
748 
tz=is[i][ZZ]; 
749 

750 
x[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX]; 
751 
x[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY]; 
752 
x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ]; 
753 
} 
754 
} else {

755 
for(i=0,j=g0; (i<gn); i++,j++) { 
756 
tx=is[i][XX]; 
757 
ty=is[i][YY]; 
758 
tz=is[i][ZZ]; 
759 

760 
x[j][XX]=x[j][XX]+tx*box[XX][XX]; 
761 
x[j][YY]=x[j][YY]+ty*box[YY][YY]; 
762 
x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ]; 
763 
} 
764 
} 
765 

766 
} 
767  
768 
void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])

769 
{ 
770 
ivec *is; 
771 
int g0,gn;

772 
int i,j,tx,ty,tz;

773  
774 
if (g>bScrewPBC)

775 
gmx_incons("screw pbc not implemented for unshift_x");

776  
777 
g0=g>start; 
778 
gn=g>nnodes; 
779 
is=g>ishift; 
780 
if(TRICLINIC(box)) {

781 
for(i=0,j=g0; (i<gn); i++,j++) { 
782 
tx=is[i][XX]; 
783 
ty=is[i][YY]; 
784 
tz=is[i][ZZ]; 
785 

786 
x[j][XX]=x_s[j][XX]tx*box[XX][XX]ty*box[YY][XX]tz*box[ZZ][XX]; 
787 
x[j][YY]=x_s[j][YY]ty*box[YY][YY]tz*box[ZZ][YY]; 
788 
x[j][ZZ]=x_s[j][ZZ]tz*box[ZZ][ZZ]; 
789 
} 
790 
} else {

791 
for(i=0,j=g0; (i<gn); i++,j++) { 
792 
tx=is[i][XX]; 
793 
ty=is[i][YY]; 
794 
tz=is[i][ZZ]; 
795 

796 
x[j][XX]=x_s[j][XX]tx*box[XX][XX]; 
797 
x[j][YY]=x_s[j][YY]ty*box[YY][YY]; 
798 
x[j][ZZ]=x_s[j][ZZ]tz*box[ZZ][ZZ]; 
799 
} 
800 
} 
801 
} 
802  
803 
void unshift_self(t_graph *g,matrix box,rvec x[])

804 
{ 
805 
ivec *is; 
806 
int g0,gn;

807 
int i,j,tx,ty,tz;

808  
809 
if (g>bScrewPBC)

810 
gmx_incons("screw pbc not implemented for unshift_self");

811  
812 
g0=g>start; 
813 
gn=g>nnodes; 
814 
is=g>ishift; 
815 
if(TRICLINIC(box)) {

816 
for(i=0,j=g0; (i<gn); i++,j++) { 
817 
tx=is[i][XX]; 
818 
ty=is[i][YY]; 
819 
tz=is[i][ZZ]; 
820 

821 
x[j][XX]=x[j][XX]tx*box[XX][XX]ty*box[YY][XX]tz*box[ZZ][XX]; 
822 
x[j][YY]=x[j][YY]ty*box[YY][YY]tz*box[ZZ][YY]; 
823 
x[j][ZZ]=x[j][ZZ]tz*box[ZZ][ZZ]; 
824 
} 
825 
} else {

826 
for(i=0,j=g0; (i<gn); i++,j++) { 
827 
tx=is[i][XX]; 
828 
ty=is[i][YY]; 
829 
tz=is[i][ZZ]; 
830 

831 
x[j][XX]=x[j][XX]tx*box[XX][XX]; 
832 
x[j][YY]=x[j][YY]ty*box[YY][YY]; 
833 
x[j][ZZ]=x[j][ZZ]tz*box[ZZ][ZZ]; 
834 
} 
835 
} 
836 
} 
837 
#undef GCHECK

838  
839 
#ifdef DEBUGMSHIFT

840 
void main(int argc,char *argv[]) 
841 
{ 
842 
FILE *out; 
843 
t_args targ; 
844 
t_topology top; 
845 
t_statheader sh; 
846 
rvec *x; 
847 
ivec *mshift; 
848 
matrix box; 
849  
850 
t_graph *g; 
851 
int i,idum,pid;

852 
real rdum; 
853  
854 
CopyRight(stderr,argv[0]);

855 
parse_common_args(&argc,argv,&targ,PCA_NEED_INOUT,NULL);

856 
if (argc > 1) 
857 
pid=atoi(argv[1]);

858 
else

859 
pid=0;

860 

861 
read_status_header(targ.infile,&sh); 
862 
snew(x,sh.natoms); 
863 
snew(mshift,sh.natoms); 
864  
865 
fprintf(stderr,"Reading Status %s\n",targ.infile);

866 
read_status(targ.infile,&idum,&rdum,&rdum,NULL,

867 
box,NULL,NULL,&idum,x,NULL,NULL,&idum,NULL,&top); 
868  
869 
fprintf(stderr,"Making Graph Structure...\n");

870 
g=mk_graph(&(top.idef),top.atoms.nr,FALSE,FALSE); 
871  
872 
out=gmx_fio_fopen(targ.outfile,"w");

873  
874 
fprintf(stderr,"Making Shift...\n");

875 
mk_mshift(out,g,box,x,mshift); 
876  
877 
p_graph(out,"In Den Haag daar woont een graaf...",g);

878 
gmx_fio_fclose(out); 
879 
} 
880 
#endif

881 