3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
41 #include "gmx_fatal.h"
50 /************************************************************
52 * S H I F T U T I L I T I E S
54 ************************************************************/
57 /************************************************************
59 * G R A P H G E N E R A T I O N C O D E
61 ************************************************************/
63 static void add_gbond(t_graph *g,atom_id a0,atom_id a1)
69 inda0 = a0 - g->at_start;
70 inda1 = a1 - g->at_start;
72 /* Search for a direct edge between a0 and a1.
73 * All egdes are bidirectional, so we only need to search one way.
75 for(i=0; (i<g->nedge[inda0] && !bFound); i++) {
76 bFound = (g->edge[inda0][i] == a1);
80 g->edge[inda0][g->nedge[inda0]++] = a1;
81 g->edge[inda1][g->nedge[inda1]++] = a0;
85 static void mk_igraph(t_graph *g,int ftype,t_ilist *il,
86 int at_start,int at_end,
98 np = interaction_function[ftype].nratoms;
100 if (ia[1] >= at_start && ia[1] < at_end) {
101 if (ia[np] >= at_end)
103 "Molecule in topology has atom numbers below and "
104 "above natoms (%d).\n"
105 "You are probably trying to use a trajectory which does "
106 "not match the first %d atoms of the run input file.\n"
107 "You can make a matching run input file with tpbconv.",
109 if (ftype == F_SETTLE) {
110 /* Bond all the atoms in the settle */
111 add_gbond(g,ia[1],ia[2]);
112 add_gbond(g,ia[1],ia[3]);
113 } else if (part == NULL) {
114 /* Simply add this bond */
115 for(j=1; j<np; j++) {
116 add_gbond(g,ia[j],ia[j+1]);
119 /* Add this bond when it connects two unlinked parts of the graph */
120 for(j=1; j<np; j++) {
121 if (part[ia[j]] != part[ia[j+1]]) {
122 add_gbond(g,ia[j],ia[j+1]);
132 GMX_ATTRIBUTE_NORETURN static void g_error(int line,const char *file)
134 gmx_fatal(FARGS,"Tring to print non existant graph (file %s, line %d)",
137 #define GCHECK(g) if (g == NULL) g_error(__LINE__,__FILE__)
139 void p_graph(FILE *log,const char *title,t_graph *g)
142 const char *cc[egcolNR] = { "W", "G", "B" };
145 fprintf(log,"graph: %s\n",title);
146 fprintf(log,"nnodes: %d\n",g->nnodes);
147 fprintf(log,"nbound: %d\n",g->nbound);
148 fprintf(log,"start: %d\n",g->at_start);
149 fprintf(log,"end: %d\n",g->at_end);
150 fprintf(log," atom shiftx shifty shiftz C nedg e1 e2 etc.\n");
151 for(i=0; (i<g->nnodes); i++)
152 if (g->nedge[i] > 0) {
153 fprintf(log,"%5d%7d%7d%7d %1s%5d",g->at_start+i+1,
154 g->ishift[g->at_start+i][XX],
155 g->ishift[g->at_start+i][YY],
156 g->ishift[g->at_start+i][ZZ],
157 (g->negc > 0) ? cc[g->egc[i]] : " ",
159 for(j=0; (j<g->nedge[i]); j++)
160 fprintf(log," %5u",g->edge[i][j]+1);
166 static void calc_1se(t_graph *g,int ftype,t_ilist *il,
167 int nbond[],int at_start,int at_end)
175 for(j=0; (j<end); j+=nratoms+1,ia+=nratoms+1) {
176 nratoms = interaction_function[ftype].nratoms;
178 if (ftype == F_SETTLE) {
180 if (iaa >= at_start && iaa < at_end) {
184 g->at_start = min(g->at_start,iaa);
185 g->at_end = max(g->at_end,iaa+2+1);
188 for(k=1; (k<=nratoms); k++) {
190 if (iaa >= at_start && iaa < at_end) {
191 g->at_start = min(g->at_start,iaa);
192 g->at_end = max(g->at_end, iaa+1);
193 /* When making the graph we (might) link all atoms in an interaction
194 * sequentially. Therefore the end atoms add 1 to the count,
195 * the middle atoms 2.
197 if (k == 1 || k == nratoms) {
208 static int calc_start_end(FILE *fplog,t_graph *g,t_ilist il[],
209 int at_start,int at_end,
214 g->at_start = at_end;
217 /* First add all the real bonds: they should determine the molecular
220 for(i=0; (i<F_NRE); i++)
221 if (interaction_function[i].flags & IF_CHEMBOND)
222 calc_1se(g,i,&il[i],nbond,at_start,at_end);
223 /* Then add all the other interactions in fixed lists, but first
224 * check to see what's there already.
226 for(i=0; (i<F_NRE); i++) {
227 if (!(interaction_function[i].flags & IF_CHEMBOND)) {
228 calc_1se(g,i,&il[i],nbond,at_start,at_end);
234 for(i=g->at_start; (i<g->at_end); i++) {
236 nnb = max(nnb,nbond[i]);
239 fprintf(fplog,"Max number of connections per atom is %d\n",nnb);
240 fprintf(fplog,"Total number of connections is %d\n",nbtot);
247 static void compact_graph(FILE *fplog,t_graph *g)
254 for(i=0; i<g->nnodes; i++) {
255 for(j=0; j<g->nedge[i]; j++) {
256 g->edge[0][n++] = g->edge[i][j];
258 max_nedge = max(max_nedge,g->nedge[i]);
260 srenew(g->edge[0],n);
261 /* set pointers after srenew because edge[0] might move */
262 for(i=1; i<g->nnodes; i++) {
263 g->edge[i] = g->edge[i-1] + g->nedge[i-1];
267 fprintf(fplog,"Max number of graph edges per atom is %d\n",
269 fprintf(fplog,"Total number of graph edges is %d\n",n);
273 static gmx_bool determine_graph_parts(t_graph *g,int *part)
280 /* Initialize the part array with all entries different */
281 for(at_i=g->at_start; at_i<g->at_end; at_i++) {
285 /* Loop over the graph until the part array is fixed */
289 for(i=0; (i<g->nnodes); i++) {
290 at_i = g->at_start + i;
292 for(e=0; e<g->nedge[i]; e++) {
293 /* Set part for both nodes to the minimum */
294 if (part[at_i2[e]] > part[at_i]) {
295 part[at_i2[e]] = part[at_i];
297 } else if (part[at_i2[e]] < part[at_i]) {
298 part[at_i] = part[at_i2[e]];
302 if (part[at_i] != part[g->at_start]) {
307 fprintf(debug,"graph part[] nchanged=%d, bMultiPart=%d\n",
308 nchanged,bMultiPart);
310 } while (nchanged > 0);
315 void mk_graph_ilist(FILE *fplog,
316 t_ilist *ilist,int at_start,int at_end,
317 gmx_bool bShakeOnly,gmx_bool bSettle,
325 gmx_incons("In mk_graph_ilist at_start can not be != 0");
330 nbtot = calc_start_end(fplog,g,ilist,at_start,at_end,nbond);
332 if (g->at_start >= g->at_end) {
333 g->at_start = at_start;
339 g->nnodes = g->at_end - g->at_start;
340 snew(g->nedge,g->nnodes);
341 snew(g->edge,g->nnodes);
342 /* Allocate a single array and set pointers into it */
343 snew(g->edge[0],nbtot);
344 for(i=1; (i<g->nnodes); i++) {
345 g->edge[i] = g->edge[i-1] + nbond[g->at_start+i-1];
349 /* First add all the real bonds: they should determine the molecular
352 for(i=0; (i<F_NRE); i++)
353 if (interaction_function[i].flags & IF_CHEMBOND)
354 mk_igraph(g,i,&(ilist[i]),at_start,at_end,NULL);
356 /* Determine of which separated parts the IF_CHEMBOND graph consists.
357 * Store the parts in the nbond array.
359 bMultiPart = determine_graph_parts(g,nbond);
362 /* Then add all the other interactions in fixed lists,
363 * but only when they connect parts of the graph
364 * that are not connected through IF_CHEMBOND interactions.
366 for(i=0; (i<F_NRE); i++) {
367 if (!(interaction_function[i].flags & IF_CHEMBOND)) {
368 mk_igraph(g,i,&(ilist[i]),at_start,at_end,nbond);
373 /* Removed all the unused space from the edge array */
374 compact_graph(fplog,g);
377 /* This is a special thing used in splitter.c to generate shake-blocks */
378 mk_igraph(g,F_CONSTR,&(ilist[F_CONSTR]),at_start,at_end,NULL);
380 mk_igraph(g,F_SETTLE,&(ilist[F_SETTLE]),at_start,at_end,NULL);
383 for(i=0; (i<g->nnodes); i++)
393 snew(g->ishift,g->natoms);
396 p_graph(debug,"graph",g);
399 t_graph *mk_graph(FILE *fplog,
400 t_idef *idef,int at_start,int at_end,
401 gmx_bool bShakeOnly,gmx_bool bSettle)
407 mk_graph_ilist(fplog,idef->il,at_start,at_end,bShakeOnly,bSettle,g);
412 void done_graph(t_graph *g)
426 /************************************************************
428 * S H I F T C A L C U L A T I O N C O D E
430 ************************************************************/
432 static void mk_1shift_tric(int npbcdim,matrix box,rvec hbox,
433 rvec xi,rvec xj,int *mi,int *mj)
435 /* Calculate periodicity for triclinic box... */
442 for(m=npbcdim-1; (m>=0); m--) {
443 /* If dx < hbox, then xj will be reduced by box, so that
444 * xi - xj will be bigger
446 if (dx[m] < -hbox[m]) {
448 for(d=m-1; d>=0; d--)
450 } else if (dx[m] >= hbox[m]) {
452 for(d=m-1; d>=0; d--)
459 static void mk_1shift(int npbcdim,rvec hbox,rvec xi,rvec xj,int *mi,int *mj)
461 /* Calculate periodicity for rectangular box... */
468 for(m=0; (m<npbcdim); m++) {
469 /* If dx < hbox, then xj will be reduced by box, so that
470 * xi - xj will be bigger
472 if (dx[m] < -hbox[m])
474 else if (dx[m] >= hbox[m])
481 static void mk_1shift_screw(matrix box,rvec hbox,
482 rvec xi,rvec xj,int *mi,int *mj)
484 /* Calculate periodicity for rectangular box... */
488 if ((mi[XX] > 0 && mi[XX] % 2 == 1) ||
489 (mi[XX] < 0 && -mi[XX] % 2 == 1)) {
497 if (dx[XX] < -hbox[XX])
499 else if (dx[XX] >= hbox[XX])
503 if (mj[XX] != mi[XX]) {
505 dx[YY] = xi[YY] - (box[YY][YY] + box[ZZ][YY] - xj[YY]);
506 dx[ZZ] = xi[ZZ] - (box[ZZ][ZZ] - xj[ZZ]);
508 for(m=1; (m<DIM); m++) {
509 /* The signs are taken such that we can first shift x and rotate
510 * and then shift y and z.
512 if (dx[m] < -hbox[m])
513 mj[m] = mi[m] - signi;
514 else if (dx[m] >= hbox[m])
515 mj[m] = mi[m] + signi;
521 static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI,
522 int npbcdim,matrix box,rvec x[],int *nerror)
530 for(m=0; (m<DIM); m++)
531 hbox[m]=box[m][m]*0.5;
532 bTriclinic = TRICLINIC(box);
538 /* Loop over all the bonds */
539 for(j=0; (j<g->nedge[ai-g0]); j++) {
540 aj = g->edge[ai-g0][j];
541 /* If there is a white one, make it grey and set pbc */
543 mk_1shift_screw(box,hbox,x[ai],x[aj],g->ishift[ai],is_aj);
545 mk_1shift_tric(npbcdim,box,hbox,x[ai],x[aj],g->ishift[ai],is_aj);
547 mk_1shift(npbcdim,hbox,x[ai],x[aj],g->ishift[ai],is_aj);
549 if (egc[aj-g0] == egcolWhite) {
550 if (aj - g0 < *AtomI)
552 egc[aj-g0] = egcolGrey;
554 copy_ivec(is_aj,g->ishift[aj]);
558 else if ((is_aj[XX] != g->ishift[aj][XX]) ||
559 (is_aj[YY] != g->ishift[aj][YY]) ||
560 (is_aj[ZZ] != g->ishift[aj][ZZ])) {
562 set_pbc(&pbc,-1,box);
563 pbc_dx(&pbc,x[ai],x[aj],dx);
564 fprintf(debug,"mk_grey: shifts for atom %d due to atom %d\n"
565 "are (%d,%d,%d), should be (%d,%d,%d)\n"
567 aj+1,ai+1,is_aj[XX],is_aj[YY],is_aj[ZZ],
568 g->ishift[aj][XX],g->ishift[aj][YY],g->ishift[aj][ZZ],
569 dx[XX],dx[YY],dx[ZZ]);
577 static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[])
578 /* Return the first node with colour Col starting at fC.
579 * return -1 if none found.
584 for(i=fC; (i<g->nnodes); i++)
585 if ((g->nedge[i] > 0) && (egc[i]==Col))
591 void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[])
593 static int nerror_tot = 0;
596 int nW,nG,nB; /* Number of Grey, Black, White */
597 int fW,fG; /* First of each category */
600 g->bScrewPBC = (ePBC == epbcSCREW);
608 /* This puts everything in the central box, that is does not move it
609 * at all. If we return without doing this for a system without bonds
610 * (i.e. only settles) all water molecules are moved to the opposite octant
612 for(i=0; (i<g->natoms); i++) {
613 g->ishift[i][XX]=g->ishift[i][YY]=g->ishift[i][ZZ]=0;
620 if (nnodes > g->negc) {
622 srenew(g->egc,g->negc);
624 memset(g->egc,0,(size_t)(nnodes*sizeof(g->egc[0])));
632 /* We even have a loop invariant:
633 * nW+nG+nB == g->nbound
636 fprintf(log,"Starting W loop\n");
639 /* Find the first white, this will allways be a larger
640 * number than before, because no nodes are made white
643 if ((fW=first_colour(fW,egcolWhite,g,g->egc)) == -1)
644 gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW);
646 /* Make the first white node grey */
647 g->egc[fW]=egcolGrey;
651 /* Initial value for the first grey */
654 fprintf(log,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
658 if ((fG=first_colour(fG,egcolGrey,g,g->egc)) == -1)
659 gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG);
661 /* Make the first grey node black */
662 g->egc[fG]=egcolBlack;
666 /* Make all the neighbours of this black node grey
667 * and set their periodicity
669 ng=mk_grey(log,nnodes,g->egc,g,&fG,npbcdim,box,x,&nerror);
670 /* ng is the number of white nodes made grey */
677 if (nerror_tot <= 100) {
678 fprintf(stderr,"There were %d inconsistent shifts. Check your topology\n",
681 fprintf(log,"There were %d inconsistent shifts. Check your topology\n",
685 if (nerror_tot == 100) {
686 fprintf(stderr,"Will stop reporting inconsistent shifts\n");
688 fprintf(log,"Will stop reporting inconsistent shifts\n");
694 /************************************************************
696 * A C T U A L S H I F T C O D E
698 ************************************************************/
700 void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
711 for(j=0; j<g0; j++) {
712 copy_rvec(x[j],x_s[j]);
716 for(j=g0; (j<g1); j++) {
721 if ((tx > 0 && tx % 2 == 1) ||
722 (tx < 0 && -tx %2 == 1)) {
723 x_s[j][XX] = x[j][XX] + tx*box[XX][XX];
724 x_s[j][YY] = box[YY][YY] + box[ZZ][YY] - x[j][YY];
725 x_s[j][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
727 x_s[j][XX] = x[j][XX];
729 x_s[j][YY] = x[j][YY] + ty*box[YY][YY] + tz*box[ZZ][YY];
730 x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ];
732 } else if (TRICLINIC(box)) {
733 for(j=g0; (j<g1); j++) {
738 x_s[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
739 x_s[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
740 x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
743 for(j=g0; (j<g1); j++) {
748 x_s[j][XX]=x[j][XX]+tx*box[XX][XX];
749 x_s[j][YY]=x[j][YY]+ty*box[YY][YY];
750 x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
754 for(j=g1; j<g->natoms; j++) {
755 copy_rvec(x[j],x_s[j]);
759 void shift_self(t_graph *g,matrix box,rvec x[])
766 gmx_incons("screw pbc not implemented for shift_self");
773 fprintf(stderr,"Shifting atoms %d to %d\n",g0,g0+gn);
776 for(j=g0; (j<g1); j++) {
781 x[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
782 x[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
783 x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
786 for(j=g0; (j<g1); j++) {
791 x[j][XX]=x[j][XX]+tx*box[XX][XX];
792 x[j][YY]=x[j][YY]+ty*box[YY][YY];
793 x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
798 void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
805 gmx_incons("screw pbc not implemented (yet) for unshift_x");
811 for(j=0; j<g0; j++) {
812 copy_rvec(x_s[j],x[j]);
816 for(j=g0; (j<g1); j++) {
821 x[j][XX]=x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
822 x[j][YY]=x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
823 x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
826 for(j=g0; (j<g1); j++) {
831 x[j][XX]=x_s[j][XX]-tx*box[XX][XX];
832 x[j][YY]=x_s[j][YY]-ty*box[YY][YY];
833 x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
837 for(j=g1; j<g->natoms; j++) {
838 copy_rvec(x_s[j],x[j]);
842 void unshift_self(t_graph *g,matrix box,rvec x[])
849 gmx_incons("screw pbc not implemented for unshift_self");
856 for(j=g0; (j<g1); j++) {
861 x[j][XX]=x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
862 x[j][YY]=x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
863 x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];
866 for(j=g0; (j<g1); j++) {
871 x[j][XX]=x[j][XX]-tx*box[XX][XX];
872 x[j][YY]=x[j][YY]-ty*box[YY][YY];
873 x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];