--- /dev/null
- if(NOT GMX_OPENMM)
+set(LIBGROMACS_SOURCES)
+
+add_subdirectory(legacyheaders)
+add_subdirectory(gmxlib)
+add_subdirectory(mdlib)
+add_subdirectory(gmxpreprocess)
+add_subdirectory(analysisdata)
+add_subdirectory(commandline)
+add_subdirectory(options)
+add_subdirectory(selection)
+add_subdirectory(trajectoryanalysis)
+add_subdirectory(utility)
+
+file(GLOB LIBGROMACS_HEADERS *.h)
+install(FILES ${LIBGROMACS_HEADERS} DESTINATION ${INCL_INSTALL_DIR}/gromacs
+ COMPONENT development)
+
+# only fiddle with assembly kernels if we're not doing OpenMM build
- endif(NOT GMX_OPENMM)
++if(GMX_SSEKERNEL_ASM_SRC AND NOT GMX_OPENMM)
+if(GMX_ASM_USEASM_NASM)
+ enable_language(ASM_NASM)
+ # if NASM is used, we need a special build command for windows...
+ FOREACH(SRC ${GMX_SSEKERNEL_ASM_SRC})
+ GET_FILENAME_COMPONENT(FILE_BASE ${SRC} NAME_WE)
+ SET(OBJ ${CMAKE_CURRENT_BINARY_DIR}/${FILE_BASE}${CMAKE_C_OUTPUT_EXTENSION})
+
+ ADD_CUSTOM_COMMAND(OUTPUT ${OBJ}
+ MAIN_DEPENDENCY ${SRC}
+ COMMAND ${CMAKE_ASM_NASM_COMPILER} -f ${CMAKE_ASM_NASM_OBJECT_FORMAT} -o ${OBJ} ${SRC})
+
+ SET(ALL_ASM_OBJS ${ALL_ASM_OBJS} ${OBJ})
+ ENDFOREACH(SRC ${GMX_SSEKERNEL_ASM_SRC})
+ set(GMX_SSEKERNEL_ASM_SRC ${ALL_ASM_OBJS})
+else(GMX_ASM_USEASM_NASM)
+ enable_language(ASM-ATT)
+ SET(CMAKE_ASM-ATT_COMPILER ${CMAKE_C_COMPILER})
+ if(GMX_IA32_ASM)
+ set_source_files_properties(${GMX_SSEKERNEL_ASM_SRC} PROPERTIES COMPILE_FLAGS "-c -m32")
+ else()
+ set_source_files_properties(${GMX_SSEKERNEL_ASM_SRC} PROPERTIES COMPILE_FLAGS "-c -m64")
+ endif()
+endif(GMX_ASM_USEASM_NASM)
++endif(GMX_SSEKERNEL_ASM_SRC AND NOT GMX_OPENMM)
+
+list(APPEND LIBGROMACS_SOURCES ${GMXLIB_SOURCES} ${GMX_SSEKERNEL_ASM_SRC} ${MDLIB_SOURCES})
+
+# add target that generates version.c every time a make is run
+# only do this if we generate the version
+if (USE_VERSION_H)
+ add_custom_target(gmx_version ALL
+ COMMAND ${CMAKE_COMMAND}
+ -D GIT_EXECUTABLE="${GIT_EXECUTABLE}"
+ -D GIT_VERSION="${GIT_VERSION}"
+ -D PROJECT_VERSION="${PROJECT_VERSION}"
+ -D PROJECT_SOURCE_DIR="${PROJECT_SOURCE_DIR}"
+ -D VERSION_C_CMAKEIN="${CMAKE_CURRENT_SOURCE_DIR}/version.c.cmakein"
+ -D VERSION_C_OUT="${CMAKE_CURRENT_BINARY_DIR}/version.c"
+ -P ${CMAKE_SOURCE_DIR}/cmake/gmxGenerateVersionInfo.cmake
+ WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/src/gmxlib
+ DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/version.c.cmakein
+ COMMENT "Generating version information")
+ list(APPEND LIBGROMACS_SOURCES ${CMAKE_CURRENT_BINARY_DIR}/version.c) # auto-generated
+ set_source_files_properties(${CMAKE_CURRENT_BINARY_DIR}/version.c
+ PROPERTIES GENERATED true)
+endif (USE_VERSION_H)
+
+add_library(libgromacs ${LIBGROMACS_SOURCES})
+if (USE_VERSION_H)
+ add_dependencies(libgromacs gmx_version)
+endif (USE_VERSION_H)
+target_link_libraries(libgromacs
+ ${GMX_EXTRA_LIBRARIES} ${FFT_LIBRARIES} ${XML_LIBRARIES}
+ ${THREAD_LIB})
+set_target_properties(libgromacs PROPERTIES
+ OUTPUT_NAME "gromacs${GMX_LIBS_SUFFIX}"
+ SOVERSION ${SOVERSION}
+ INSTALL_NAME_DIR "${LIB_INSTALL_DIR}")
+
+install(TARGETS libgromacs DESTINATION ${LIB_INSTALL_DIR} COMPONENT libraries)
+
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/libgromacs.pc.cmakein
+ ${CMAKE_CURRENT_BINARY_DIR}/libgromacs.pc @ONLY)
+install(FILES ${CMAKE_CURRENT_BINARY_DIR}/libgromacs.pc
+ DESTINATION ${LIB_INSTALL_DIR}/pkgconfig
+ RENAME "libgromacs${GMX_LIBS_SUFFIX}.pc"
+ COMPONENT development)
--- /dev/null
- for(i=0; i<graph->nnodes; i++) {
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * GROningen Mixture of Alchemy and Childrens' Stories
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include "sysstuff.h"
+#include "typedefs.h"
+#include "vec.h"
+#include "maths.h"
+#include "main.h"
+#include "pbc.h"
+#include "smalloc.h"
+#include "txtdump.h"
+#include "gmx_fatal.h"
+#include "names.h"
+#include "macros.h"
+
+/* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
+enum { epbcdxRECTANGULAR=1, epbcdxTRICLINIC,
+ epbcdx2D_RECT, epbcdx2D_TRIC,
+ epbcdx1D_RECT, epbcdx1D_TRIC,
+ epbcdxSCREW_RECT, epbcdxSCREW_TRIC,
+ epbcdxNOPBC, epbcdxUNSUPPORTED };
+
+/* Margin factor for error message and correction if the box is too skewed */
+#define BOX_MARGIN 1.0010
+#define BOX_MARGIN_CORRECT 1.0005
+
+int ePBC2npbcdim(int ePBC)
+{
+ int npbcdim=0;
+
+ switch(ePBC) {
+ case epbcXYZ: npbcdim = 3; break;
+ case epbcXY: npbcdim = 2; break;
+ case epbcSCREW: npbcdim = 3; break;
+ case epbcNONE: npbcdim = 0; break;
+ default: gmx_fatal(FARGS,"Unknown ePBC=%d in ePBC2npbcdim",ePBC);
+ }
+
+ return npbcdim;
+}
+
+int inputrec2nboundeddim(t_inputrec *ir)
+{
+ if (ir->nwall == 2 && ir->ePBC == epbcXY)
+ {
+ return 3;
+ }
+ else
+ {
+ return ePBC2npbcdim(ir->ePBC);
+ }
+}
+
+void dump_pbc(FILE *fp,t_pbc *pbc)
+{
+ rvec sum_box;
+
+ fprintf(fp,"ePBCDX = %d\n",pbc->ePBCDX);
+ pr_rvecs(fp,0,"box",pbc->box,DIM);
+ pr_rvecs(fp,0,"fbox_diag",&pbc->fbox_diag,1);
+ pr_rvecs(fp,0,"hbox_diag",&pbc->hbox_diag,1);
+ pr_rvecs(fp,0,"mhbox_diag",&pbc->mhbox_diag,1);
+ rvec_add(pbc->hbox_diag,pbc->mhbox_diag,sum_box);
+ pr_rvecs(fp,0,"sum of the above two",&sum_box,1);
+ fprintf(fp,"max_cutoff2 = %g\n",pbc->max_cutoff2);
+ fprintf(fp,"bLimitDistance = %s\n",BOOL(pbc->bLimitDistance));
+ fprintf(fp,"limit_distance2 = %g\n",pbc->limit_distance2);
+ fprintf(fp,"ntric_vec = %d\n",pbc->ntric_vec);
+ if (pbc->ntric_vec > 0) {
+ pr_ivecs(fp,0,"tric_shift",pbc->tric_shift,pbc->ntric_vec,FALSE);
+ pr_rvecs(fp,0,"tric_vec",pbc->tric_vec,pbc->ntric_vec);
+ }
+}
+
+const char *check_box(int ePBC,matrix box)
+{
+ const char *ptr;
+
+ if (ePBC == -1)
+ ePBC = guess_ePBC(box);
+
+ if (ePBC == epbcNONE)
+ return NULL;
+
+ if ((box[XX][YY] != 0) || (box[XX][ZZ] != 0) || (box[YY][ZZ] != 0)) {
+ ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
+ } else if (ePBC == epbcSCREW && (box[YY][XX] != 0 || box[ZZ][XX] != 0)) {
+ ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
+ } else if (fabs(box[YY][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
+ (ePBC != epbcXY &&
+ (fabs(box[ZZ][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
+ fabs(box[ZZ][YY]) > BOX_MARGIN*0.5*box[YY][YY]))) {
+ ptr = "Triclinic box is too skewed.";
+ } else {
+ ptr = NULL;
+ }
+
+ return ptr;
+}
+
+real max_cutoff2(int ePBC,matrix box)
+{
+ real min_hv2,min_ss;
+
+ /* Physical limitation of the cut-off
+ * by half the length of the shortest box vector.
+ */
+ min_hv2 = min(0.25*norm2(box[XX]),0.25*norm2(box[YY]));
+ if (ePBC != epbcXY)
+ min_hv2 = min(min_hv2,0.25*norm2(box[ZZ]));
+
+ /* Limitation to the smallest diagonal element due to optimizations:
+ * checking only linear combinations of single box-vectors (2 in x)
+ * in the grid search and pbc_dx is a lot faster
+ * than checking all possible combinations.
+ */
+ if (ePBC == epbcXY) {
+ min_ss = min(box[XX][XX],box[YY][YY]);
+ } else {
+ min_ss = min(box[XX][XX],min(box[YY][YY]-fabs(box[ZZ][YY]),box[ZZ][ZZ]));
+ }
+
+ return min(min_hv2,min_ss*min_ss);
+}
+
+/* this one is mostly harmless... */
+static gmx_bool bWarnedGuess=FALSE;
+
+int guess_ePBC(matrix box)
+{
+ int ePBC;
+
+ if (box[XX][XX]>0 && box[YY][YY]>0 && box[ZZ][ZZ]>0) {
+ ePBC = epbcXYZ;
+ } else if (box[XX][XX]>0 && box[YY][YY]>0 && box[ZZ][ZZ]==0) {
+ ePBC = epbcXY;
+ } else if (box[XX][XX]==0 && box[YY][YY]==0 && box[ZZ][ZZ]==0) {
+ ePBC = epbcNONE;
+ } else {
+ if (!bWarnedGuess) {
+ fprintf(stderr,"WARNING: Unsupported box diagonal %f %f %f, "
+ "will not use periodic boundary conditions\n\n",
+ box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
+ bWarnedGuess = TRUE;
+ }
+ ePBC = epbcNONE;
+ }
+
+ if (debug)
+ fprintf(debug,"Guessed pbc = %s from the box matrix\n",epbc_names[ePBC]);
+
+ return ePBC;
+}
+
+static int correct_box_elem(FILE *fplog,int step,tensor box,int v,int d)
+{
+ int shift,maxshift=10;
+
+ shift = 0;
+
+ /* correct elem d of vector v with vector d */
+ while (box[v][d] > BOX_MARGIN_CORRECT*0.5*box[d][d]) {
+ if (fplog) {
+ fprintf(fplog,"Step %d: correcting invalid box:\n",step);
+ pr_rvecs(fplog,0,"old box",box,DIM);
+ }
+ rvec_dec(box[v],box[d]);
+ shift--;
+ if (fplog) {
+ pr_rvecs(fplog,0,"new box",box,DIM);
+ }
+ if (shift <= -maxshift)
+ gmx_fatal(FARGS,
+ "Box was shifted at least %d times. Please see log-file.",
+ maxshift);
+ }
+ while (box[v][d] < -BOX_MARGIN_CORRECT*0.5*box[d][d]) {
+ if (fplog) {
+ fprintf(fplog,"Step %d: correcting invalid box:\n",step);
+ pr_rvecs(fplog,0,"old box",box,DIM);
+ }
+ rvec_inc(box[v],box[d]);
+ shift++;
+ if (fplog) {
+ pr_rvecs(fplog,0,"new box",box,DIM);
+ }
+ if (shift >= maxshift)
+ gmx_fatal(FARGS,
+ "Box was shifted at least %d times. Please see log-file.",
+ maxshift);
+ }
+
+ return shift;
+}
+
+gmx_bool correct_box(FILE *fplog,int step,tensor box,t_graph *graph)
+{
+ int zy,zx,yx,i;
+ gmx_bool bCorrected;
+
+ /* check if the box still obeys the restrictions, if not, correct it */
+ zy = correct_box_elem(fplog,step,box,ZZ,YY);
+ zx = correct_box_elem(fplog,step,box,ZZ,XX);
+ yx = correct_box_elem(fplog,step,box,YY,XX);
+
+ bCorrected = (zy || zx || yx);
+
+ if (bCorrected && graph) {
+ /* correct the graph */
++ for(i=graph->at_start; i<graph->at_end; i++) {
+ graph->ishift[i][YY] -= graph->ishift[i][ZZ]*zy;
+ graph->ishift[i][XX] -= graph->ishift[i][ZZ]*zx;
+ graph->ishift[i][XX] -= graph->ishift[i][YY]*yx;
+ }
+ }
+
+ return bCorrected;
+}
+
+int ndof_com(t_inputrec *ir)
+{
+ int n=0;
+
+ switch (ir->ePBC) {
+ case epbcXYZ:
+ case epbcNONE:
+ n = 3;
+ break;
+ case epbcXY:
+ n = (ir->nwall == 0 ? 3 : 2);
+ break;
+ case epbcSCREW:
+ n = 1;
+ break;
+ default:
+ gmx_incons("Unknown pbc in calc_nrdf");
+ }
+
+ return n;
+}
+
+static void low_set_pbc(t_pbc *pbc,int ePBC,ivec *dd_nc,matrix box)
+{
+ int order[5]={0,-1,1,-2,2};
+ int ii,jj,kk,i,j,k,d,dd,jc,kc,npbcdim,shift;
+ ivec bPBC;
+ real d2old,d2new,d2new_c;
+ rvec trial,pos;
+ gmx_bool bXY,bUse;
+ const char *ptr;
+
+ pbc->ndim_ePBC = ePBC2npbcdim(ePBC);
+
+ copy_mat(box,pbc->box);
+ pbc->bLimitDistance = FALSE;
+ pbc->max_cutoff2 = 0;
+ pbc->dim = -1;
+
+ for(i=0; (i<DIM); i++) {
+ pbc->fbox_diag[i] = box[i][i];
+ pbc->hbox_diag[i] = pbc->fbox_diag[i]*0.5;
+ pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
+ }
+
+ ptr = check_box(ePBC,box);
+ if (ePBC == epbcNONE) {
+ pbc->ePBCDX = epbcdxNOPBC;
+ } else if (ptr) {
+ fprintf(stderr, "Warning: %s\n",ptr);
+ pr_rvecs(stderr,0," Box",box,DIM);
+ fprintf(stderr, " Can not fix pbc.\n");
+ pbc->ePBCDX = epbcdxUNSUPPORTED;
+ pbc->bLimitDistance = TRUE;
+ pbc->limit_distance2 = 0;
+ } else {
+ if (ePBC == epbcSCREW && dd_nc) {
+ /* This combinated should never appear here */
+ gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
+ }
+
+ npbcdim = 0;
+ for(i=0; i<DIM; i++) {
+ if ((dd_nc && (*dd_nc)[i] > 1) || (ePBC == epbcXY && i == ZZ)) {
+ bPBC[i] = 0;
+ } else {
+ bPBC[i] = 1;
+ npbcdim++;
+ }
+ }
+ switch (npbcdim) {
+ case 1:
+ /* 1D pbc is not an mdp option and it is therefore only used
+ * with single shifts.
+ */
+ pbc->ePBCDX = epbcdx1D_RECT;
+ for(i=0; i<DIM; i++)
+ if (bPBC[i])
+ pbc->dim = i;
+ for(i=0; i<pbc->dim; i++)
+ if (pbc->box[pbc->dim][i] != 0)
+ pbc->ePBCDX = epbcdx1D_TRIC;
+ break;
+ case 2:
+ pbc->ePBCDX = epbcdx2D_RECT;
+ for(i=0; i<DIM; i++)
+ if (!bPBC[i])
+ pbc->dim = i;
+ for(i=0; i<DIM; i++)
+ if (bPBC[i])
+ for(j=0; j<i; j++)
+ if (pbc->box[i][j] != 0)
+ pbc->ePBCDX = epbcdx2D_TRIC;
+ break;
+ case 3:
+ if (ePBC != epbcSCREW) {
+ if (TRICLINIC(box)) {
+ pbc->ePBCDX = epbcdxTRICLINIC;
+ } else {
+ pbc->ePBCDX = epbcdxRECTANGULAR;
+ }
+ } else {
+ pbc->ePBCDX = (box[ZZ][YY]==0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
+ if (pbc->ePBCDX == epbcdxSCREW_TRIC) {
+ fprintf(stderr,
+ "Screw pbc is not yet implemented for triclinic boxes.\n"
+ "Can not fix pbc.\n");
+ pbc->ePBCDX = epbcdxUNSUPPORTED;
+ }
+ }
+ break;
+ default:
+ gmx_fatal(FARGS,"Incorrect number of pbc dimensions with DD: %d",
+ npbcdim);
+ }
+ pbc->max_cutoff2 = max_cutoff2(ePBC,box);
+
+ if (pbc->ePBCDX == epbcdxTRICLINIC ||
+ pbc->ePBCDX == epbcdx2D_TRIC ||
+ pbc->ePBCDX == epbcdxSCREW_TRIC) {
+ if (debug) {
+ pr_rvecs(debug,0,"Box",box,DIM);
+ fprintf(debug,"max cutoff %.3f\n",sqrt(pbc->max_cutoff2));
+ }
+ pbc->ntric_vec = 0;
+ /* We will only use single shifts, but we will check a few
+ * more shifts to see if there is a limiting distance
+ * above which we can not be sure of the correct distance.
+ */
+ for(kk=0; kk<5; kk++) {
+ k = order[kk];
+ if (!bPBC[ZZ] && k != 0)
+ continue;
+ for(jj=0; jj<5; jj++) {
+ j = order[jj];
+ if (!bPBC[YY] && j != 0)
+ continue;
+ for(ii=0; ii<3; ii++) {
+ i = order[ii];
+ if (!bPBC[XX] && i != 0)
+ continue;
+ /* A shift is only useful when it is trilinic */
+ if (j != 0 || k != 0) {
+ d2old = 0;
+ d2new = 0;
+ for(d=0; d<DIM; d++) {
+ trial[d] = i*box[XX][d] + j*box[YY][d] + k*box[ZZ][d];
+ /* Choose the vector within the brick around 0,0,0 that
+ * will become the shortest due to shift try.
+ */
+ if (d == pbc->dim) {
+ trial[d] = 0;
+ pos[d] = 0;
+ } else {
+ if (trial[d] < 0)
+ pos[d] = min( pbc->hbox_diag[d],-trial[d]);
+ else
+ pos[d] = max(-pbc->hbox_diag[d],-trial[d]);
+ }
+ d2old += sqr(pos[d]);
+ d2new += sqr(pos[d] + trial[d]);
+ }
+ if (BOX_MARGIN*d2new < d2old) {
+ if (j < -1 || j > 1 || k < -1 || k > 1) {
+ /* Check if there is a single shift vector
+ * that decreases this distance even more.
+ */
+ jc = 0;
+ kc = 0;
+ if (j < -1 || j > 1)
+ jc = j/2;
+ if (k < -1 || k > 1)
+ kc = k/2;
+ d2new_c = 0;
+ for(d=0; d<DIM; d++)
+ d2new_c += sqr(pos[d] + trial[d]
+ - jc*box[YY][d] - kc*box[ZZ][d]);
+ if (d2new_c > BOX_MARGIN*d2new) {
+ /* Reject this shift vector, as there is no a priori limit
+ * to the number of shifts that decrease distances.
+ */
+ if (!pbc->bLimitDistance || d2new < pbc->limit_distance2)
+ pbc->limit_distance2 = d2new;
+ pbc->bLimitDistance = TRUE;
+ }
+ } else {
+ /* Check if shifts with one box vector less do better */
+ bUse = TRUE;
+ for(dd=0; dd<DIM; dd++) {
+ shift = (dd==0 ? i : (dd==1 ? j : k));
+ if (shift) {
+ d2new_c = 0;
+ for(d=0; d<DIM; d++)
+ d2new_c += sqr(pos[d] + trial[d] - shift*box[dd][d]);
+ if (d2new_c <= BOX_MARGIN*d2new)
+ bUse = FALSE;
+ }
+ }
+ if (bUse) {
+ /* Accept this shift vector. */
+ if (pbc->ntric_vec >= MAX_NTRICVEC) {
+ fprintf(stderr,"\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
+ " There is probably something wrong with your box.\n",MAX_NTRICVEC);
+ pr_rvecs(stderr,0," Box",box,DIM);
+ } else {
+ copy_rvec(trial,pbc->tric_vec[pbc->ntric_vec]);
+ pbc->tric_shift[pbc->ntric_vec][XX] = i;
+ pbc->tric_shift[pbc->ntric_vec][YY] = j;
+ pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
+ pbc->ntric_vec++;
+ }
+ }
+ }
+ if (debug) {
+ fprintf(debug," tricvec %2d = %2d %2d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
+ pbc->ntric_vec,i,j,k,
+ sqrt(d2old),sqrt(d2new),
+ trial[XX],trial[YY],trial[ZZ],
+ pos[XX],pos[YY],pos[ZZ]);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+void set_pbc(t_pbc *pbc,int ePBC,matrix box)
+{
+ if (ePBC == -1)
+ ePBC = guess_ePBC(box);
+
+ low_set_pbc(pbc,ePBC,NULL,box);
+}
+
+t_pbc *set_pbc_dd(t_pbc *pbc,int ePBC,
+ gmx_domdec_t *dd,gmx_bool bSingleDir,matrix box)
+{
+ ivec nc2;
+ int npbcdim,i;
+
+ if (dd == NULL) {
+ npbcdim = DIM;
+ } else {
+ if (ePBC == epbcSCREW && dd->nc[XX] > 1) {
+ /* The rotation has been taken care of during coordinate communication */
+ ePBC = epbcXYZ;
+ }
+ npbcdim = 0;
+ for(i=0; i<DIM; i++) {
+ if (dd->nc[i] <= (bSingleDir ? 1 : 2)) {
+ nc2[i] = 1;
+ if (!(ePBC == epbcXY && i == ZZ))
+ npbcdim++;
+ } else {
+ nc2[i] = dd->nc[i];
+ }
+ }
+ }
+
+ if (npbcdim > 0)
+ low_set_pbc(pbc,ePBC,npbcdim<DIM ? &nc2 : NULL,box);
+
+ return (npbcdim > 0 ? pbc : NULL);
+}
+
+void pbc_dx(const t_pbc *pbc,const rvec x1, const rvec x2, rvec dx)
+{
+ int i,j;
+ rvec dx_start,trial;
+ real d2min,d2trial;
+ gmx_bool bRot;
+
+ rvec_sub(x1,x2,dx);
+
+ switch (pbc->ePBCDX) {
+ case epbcdxRECTANGULAR:
+ for(i=0; i<DIM; i++) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ }
+ }
+ break;
+ case epbcdxTRICLINIC:
+ for(i=DIM-1; i>=0; i--) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] -= pbc->box[i][j];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] += pbc->box[i][j];
+ }
+ }
+ /* dx is the distance in a rectangular box */
+ d2min = norm2(dx);
+ if (d2min > pbc->max_cutoff2) {
+ copy_rvec(dx,dx_start);
+ d2min = norm2(dx);
+ /* Now try all possible shifts, when the distance is within max_cutoff
+ * it must be the shortest possible distance.
+ */
+ i = 0;
+ while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec)) {
+ rvec_add(dx_start,pbc->tric_vec[i],trial);
+ d2trial = norm2(trial);
+ if (d2trial < d2min) {
+ copy_rvec(trial,dx);
+ d2min = d2trial;
+ }
+ i++;
+ }
+ }
+ break;
+ case epbcdx2D_RECT:
+ for(i=0; i<DIM; i++) {
+ if (i != pbc->dim) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ }
+ }
+ }
+ break;
+ case epbcdx2D_TRIC:
+ d2min = 0;
+ for(i=DIM-1; i>=0; i--) {
+ if (i != pbc->dim) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] -= pbc->box[i][j];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] += pbc->box[i][j];
+ }
+ d2min += dx[i]*dx[i];
+ }
+ }
+ if (d2min > pbc->max_cutoff2) {
+ copy_rvec(dx,dx_start);
+ d2min = norm2(dx);
+ /* Now try all possible shifts, when the distance is within max_cutoff
+ * it must be the shortest possible distance.
+ */
+ i = 0;
+ while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec)) {
+ rvec_add(dx_start,pbc->tric_vec[i],trial);
+ d2trial = 0;
+ for(j=0; j<DIM; j++) {
+ if (j != pbc->dim) {
+ d2trial += trial[j]*trial[j];
+ }
+ }
+ if (d2trial < d2min) {
+ copy_rvec(trial,dx);
+ d2min = d2trial;
+ }
+ i++;
+ }
+ }
+ break;
+ case epbcdxSCREW_RECT:
+ /* The shift definition requires x first */
+ bRot = FALSE;
+ while (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ bRot = !bRot;
+ }
+ while (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[YY];
+ bRot = !bRot;
+ }
+ if (bRot) {
+ /* Rotate around the x-axis in the middle of the box */
+ dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
+ dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
+ }
+ /* Normal pbc for y and z */
+ for(i=YY; i<=ZZ; i++) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ }
+ }
+ break;
+ case epbcdxNOPBC:
+ case epbcdxUNSUPPORTED:
+ break;
+ default:
+ gmx_fatal(FARGS,"Internal error in pbc_dx, set_pbc has not been called");
+ break;
+ }
+}
+
+int pbc_dx_aiuc(const t_pbc *pbc,const rvec x1, const rvec x2, rvec dx)
+{
+ int i,j,is;
+ rvec dx_start,trial;
+ real d2min,d2trial;
+ ivec ishift,ishift_start;
+
+ rvec_sub(x1,x2,dx);
+ clear_ivec(ishift);
+
+ switch (pbc->ePBCDX) {
+ case epbcdxRECTANGULAR:
+ for(i=0; i<DIM; i++) {
+ if (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ ishift[i]++;
+ }
+ }
+ break;
+ case epbcdxTRICLINIC:
+ /* For triclinic boxes the performance difference between
+ * if/else and two while loops is negligible.
+ * However, the while version can cause extreme delays
+ * before a simulation crashes due to large forces which
+ * can cause unlimited displacements.
+ * Also allowing multiple shifts would index fshift beyond bounds.
+ */
+ for(i=DIM-1; i>=1; i--) {
+ if (dx[i] > pbc->hbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] -= pbc->box[i][j];
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] += pbc->box[i][j];
+ ishift[i]++;
+ }
+ }
+ /* Allow 2 shifts in x */
+ if (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ ishift[XX]--;
+ if (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ ishift[XX]--;
+ }
+ } else if (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[XX];
+ ishift[XX]++;
+ if (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[XX];
+ ishift[XX]++;
+ }
+ }
+ /* dx is the distance in a rectangular box */
+ d2min = norm2(dx);
+ if (d2min > pbc->max_cutoff2) {
+ copy_rvec(dx,dx_start);
+ copy_ivec(ishift,ishift_start);
+ d2min = norm2(dx);
+ /* Now try all possible shifts, when the distance is within max_cutoff
+ * it must be the shortest possible distance.
+ */
+ i = 0;
+ while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec)) {
+ rvec_add(dx_start,pbc->tric_vec[i],trial);
+ d2trial = norm2(trial);
+ if (d2trial < d2min) {
+ copy_rvec(trial,dx);
+ ivec_add(ishift_start,pbc->tric_shift[i],ishift);
+ d2min = d2trial;
+ }
+ i++;
+ }
+ }
+ break;
+ case epbcdx2D_RECT:
+ for(i=0; i<DIM; i++) {
+ if (i != pbc->dim) {
+ if (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ ishift[i]++;
+ }
+ }
+ }
+ break;
+ case epbcdx2D_TRIC:
+ d2min = 0;
+ for(i=DIM-1; i>=1; i--) {
+ if (i != pbc->dim) {
+ if (dx[i] > pbc->hbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] -= pbc->box[i][j];
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] += pbc->box[i][j];
+ ishift[i]++;
+ }
+ d2min += dx[i]*dx[i];
+ }
+ }
+ if (pbc->dim != XX) {
+ /* Allow 2 shifts in x */
+ if (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ ishift[XX]--;
+ if (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ ishift[XX]--;
+ }
+ } else if (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[XX];
+ ishift[XX]++;
+ if (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[XX];
+ ishift[XX]++;
+ }
+ }
+ d2min += dx[XX]*dx[XX];
+ }
+ if (d2min > pbc->max_cutoff2) {
+ copy_rvec(dx,dx_start);
+ copy_ivec(ishift,ishift_start);
+ /* Now try all possible shifts, when the distance is within max_cutoff
+ * it must be the shortest possible distance.
+ */
+ i = 0;
+ while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec)) {
+ rvec_add(dx_start,pbc->tric_vec[i],trial);
+ d2trial = 0;
+ for(j=0; j<DIM; j++) {
+ if (j != pbc->dim) {
+ d2trial += trial[j]*trial[j];
+ }
+ }
+ if (d2trial < d2min) {
+ copy_rvec(trial,dx);
+ ivec_add(ishift_start,pbc->tric_shift[i],ishift);
+ d2min = d2trial;
+ }
+ i++;
+ }
+ }
+ break;
+ case epbcdx1D_RECT:
+ i = pbc->dim;
+ if (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ ishift[i]++;
+ }
+ break;
+ case epbcdx1D_TRIC:
+ i = pbc->dim;
+ if (dx[i] > pbc->hbox_diag[i]) {
+ rvec_dec(dx,pbc->box[i]);
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ rvec_inc(dx,pbc->box[i]);
+ ishift[i]++;
+ }
+ break;
+ case epbcdxSCREW_RECT:
+ /* The shift definition requires x first */
+ if (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ ishift[XX]--;
+ } else if (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[XX];
+ ishift[XX]++;
+ }
+ if (ishift[XX] == 1 || ishift[XX] == -1) {
+ /* Rotate around the x-axis in the middle of the box */
+ dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
+ dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
+ }
+ /* Normal pbc for y and z */
+ for(i=YY; i<=ZZ; i++) {
+ if (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ ishift[i]--;
+ } else if (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ ishift[i]++;
+ }
+ }
+ break;
+ case epbcdxNOPBC:
+ case epbcdxUNSUPPORTED:
+ break;
+ default:
+ gmx_fatal(FARGS,"Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
+ break;
+ }
+
+ is = IVEC2IS(ishift);
+ if (debug)
+ {
+ range_check_mesg(is,0,SHIFTS,"PBC shift vector index range check.");
+ }
+
+ return is;
+}
+
+void pbc_dx_d(const t_pbc *pbc,const dvec x1, const dvec x2, dvec dx)
+{
+ int i,j;
+ dvec dx_start,trial;
+ double d2min,d2trial;
+ gmx_bool bRot;
+
+ dvec_sub(x1,x2,dx);
+
+ switch (pbc->ePBCDX) {
+ case epbcdxRECTANGULAR:
+ case epbcdx2D_RECT:
+ for(i=0; i<DIM; i++) {
+ if (i != pbc->dim) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ }
+ }
+ }
+ break;
+ case epbcdxTRICLINIC:
+ case epbcdx2D_TRIC:
+ d2min = 0;
+ for(i=DIM-1; i>=0; i--) {
+ if (i != pbc->dim) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] -= pbc->box[i][j];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ for (j=i; j>=0; j--)
+ dx[j] += pbc->box[i][j];
+ }
+ d2min += dx[i]*dx[i];
+ }
+ }
+ if (d2min > pbc->max_cutoff2) {
+ copy_dvec(dx,dx_start);
+ /* Now try all possible shifts, when the distance is within max_cutoff
+ * it must be the shortest possible distance.
+ */
+ i = 0;
+ while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec)) {
+ for(j=0; j<DIM; j++) {
+ trial[j] = dx_start[j] + pbc->tric_vec[i][j];
+ }
+ d2trial = 0;
+ for(j=0; j<DIM; j++) {
+ if (j != pbc->dim) {
+ d2trial += trial[j]*trial[j];
+ }
+ }
+ if (d2trial < d2min) {
+ copy_dvec(trial,dx);
+ d2min = d2trial;
+ }
+ i++;
+ }
+ }
+ break;
+ case epbcdxSCREW_RECT:
+ /* The shift definition requires x first */
+ bRot = FALSE;
+ while (dx[XX] > pbc->hbox_diag[XX]) {
+ dx[XX] -= pbc->fbox_diag[XX];
+ bRot = !bRot;
+ }
+ while (dx[XX] <= pbc->mhbox_diag[XX]) {
+ dx[XX] += pbc->fbox_diag[YY];
+ bRot = !bRot;
+ }
+ if (bRot) {
+ /* Rotate around the x-axis in the middle of the box */
+ dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
+ dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
+ }
+ /* Normal pbc for y and z */
+ for(i=YY; i<=ZZ; i++) {
+ while (dx[i] > pbc->hbox_diag[i]) {
+ dx[i] -= pbc->fbox_diag[i];
+ }
+ while (dx[i] <= pbc->mhbox_diag[i]) {
+ dx[i] += pbc->fbox_diag[i];
+ }
+ }
+ break;
+ case epbcdxNOPBC:
+ case epbcdxUNSUPPORTED:
+ break;
+ default:
+ gmx_fatal(FARGS,"Internal error in pbc_dx, set_pbc has not been called");
+ break;
+ }
+}
+
+gmx_bool image_rect(ivec xi,ivec xj,ivec box_size,real rlong2,int *shift,real *r2)
+{
+ int m,t;
+ int dx,b,b_2;
+ real dxr,rij2;
+
+ rij2=0.0;
+ t=0;
+ for(m=0; (m<DIM); m++) {
+ dx=xi[m]-xj[m];
+ t*=DIM;
+ b=box_size[m];
+ b_2=b/2;
+ if (dx < -b_2) {
+ t+=2;
+ dx+=b;
+ }
+ else if (dx > b_2)
+ dx-=b;
+ else
+ t+=1;
+ dxr=dx;
+ rij2+=dxr*dxr;
+ if (rij2 >= rlong2)
+ return FALSE;
+ }
+
+ *shift = t;
+ *r2 = rij2;
+ return TRUE;
+}
+
+gmx_bool image_cylindric(ivec xi,ivec xj,ivec box_size,real rlong2,
+ int *shift,real *r2)
+{
+ int m,t;
+ int dx,b,b_2;
+ real dxr,rij2;
+
+ rij2=0.0;
+ t=0;
+ for(m=0; (m<DIM); m++) {
+ dx=xi[m]-xj[m];
+ t*=DIM;
+ b=box_size[m];
+ b_2=b/2;
+ if (dx < -b_2) {
+ t+=2;
+ dx+=b;
+ }
+ else if (dx > b_2)
+ dx-=b;
+ else
+ t+=1;
+
+ dxr=dx;
+ rij2+=dxr*dxr;
+ if (m < ZZ) {
+ if (rij2 >= rlong2)
+ return FALSE;
+ }
+ }
+
+ *shift = t;
+ *r2 = rij2;
+ return TRUE;
+}
+
+void calc_shifts(matrix box,rvec shift_vec[])
+{
+ int k,l,m,d,n,test;
+
+ n=0;
+ for(m = -D_BOX_Z; m <= D_BOX_Z; m++)
+ for(l = -D_BOX_Y; l <= D_BOX_Y; l++)
+ for(k = -D_BOX_X; k <= D_BOX_X; k++,n++) {
+ test = XYZ2IS(k,l,m);
+ if (n != test)
+ gmx_incons("inconsistent shift numbering");
+ for(d = 0; d < DIM; d++)
+ shift_vec[n][d] = k*box[XX][d] + l*box[YY][d] + m*box[ZZ][d];
+ }
+}
+
+void calc_box_center(int ecenter,matrix box,rvec box_center)
+{
+ int d,m;
+
+ clear_rvec(box_center);
+ switch (ecenter) {
+ case ecenterTRIC:
+ for(m=0; (m<DIM); m++)
+ for(d=0; d<DIM; d++)
+ box_center[d] += 0.5*box[m][d];
+ break;
+ case ecenterRECT:
+ for(d=0; d<DIM; d++)
+ box_center[d] = 0.5*box[d][d];
+ break;
+ case ecenterZERO:
+ break;
+ default:
+ gmx_fatal(FARGS,"Unsupported value %d for ecenter",ecenter);
+ }
+}
+
+void calc_triclinic_images(matrix box,rvec img[])
+{
+ int i;
+
+ /* Calculate 3 adjacent images in the xy-plane */
+ copy_rvec(box[0],img[0]);
+ copy_rvec(box[1],img[1]);
+ if (img[1][XX] < 0)
+ svmul(-1,img[1],img[1]);
+ rvec_sub(img[1],img[0],img[2]);
+
+ /* Get the next 3 in the xy-plane as mirror images */
+ for(i=0; i<3; i++)
+ svmul(-1,img[i],img[3+i]);
+
+ /* Calculate the first 4 out of xy-plane images */
+ copy_rvec(box[2],img[6]);
+ if (img[6][XX] < 0)
+ svmul(-1,img[6],img[6]);
+ for(i=0; i<3; i++)
+ rvec_add(img[6],img[i+1],img[7+i]);
+
+ /* Mirror the last 4 from the previous in opposite rotation */
+ for(i=0; i<4; i++)
+ svmul(-1,img[6 + (2+i) % 4],img[10+i]);
+}
+
+void calc_compact_unitcell_vertices(int ecenter,matrix box,rvec vert[])
+{
+ rvec img[NTRICIMG],box_center;
+ int n,i,j,tmp[4],d;
+
+ calc_triclinic_images(box,img);
+
+ n=0;
+ for(i=2; i<=5; i+=3) {
+ tmp[0] = i-1;
+ if (i==2)
+ tmp[1] = 8;
+ else
+ tmp[1] = 6;
+ tmp[2] = (i+1) % 6;
+ tmp[3] = tmp[1]+4;
+ for(j=0; j<4; j++) {
+ for(d=0; d<DIM; d++)
+ vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
+ n++;
+ }
+ }
+ for(i=7; i<=13; i+=6) {
+ tmp[0] = (i-7)/2;
+ tmp[1] = tmp[0]+1;
+ if (i==7)
+ tmp[2] = 8;
+ else
+ tmp[2] = 10;
+ tmp[3] = i-1;
+ for(j=0; j<4; j++) {
+ for(d=0; d<DIM; d++)
+ vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
+ n++;
+ }
+ }
+ for(i=9; i<=11; i+=2) {
+ if (i==9)
+ tmp[0] = 3;
+ else
+ tmp[0] = 0;
+ tmp[1] = tmp[0]+1;
+ if (i==9)
+ tmp[2] = 6;
+ else
+ tmp[2] = 12;
+ tmp[3] = i-1;
+ for(j=0; j<4; j++) {
+ for(d=0; d<DIM; d++)
+ vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
+ n++;
+ }
+ }
+
+ calc_box_center(ecenter,box,box_center);
+ for(i=0; i<NCUCVERT; i++)
+ for(d=0; d<DIM; d++)
+ vert[i][d] = vert[i][d]*0.25+box_center[d];
+}
+
+int *compact_unitcell_edges()
+{
+ /* this is an index in vert[] (see calc_box_vertices) */
+ /*static int edge[NCUCEDGE*2];*/
+ int *edge;
+ static const int hexcon[24] = { 0,9, 1,19, 2,15, 3,21,
+ 4,17, 5,11, 6,23, 7,13,
+ 8,20, 10,18, 12,16, 14,22 };
+ int e,i,j;
+ gmx_bool bFirst = TRUE;
+
+ snew(edge,NCUCEDGE*2);
+
+ if (bFirst) {
+ e = 0;
+ for(i=0; i<6; i++)
+ for(j=0; j<4; j++) {
+ edge[e++] = 4*i + j;
+ edge[e++] = 4*i + (j+1) % 4;
+ }
+ for(i=0; i<12*2; i++)
+ edge[e++] = hexcon[i];
+
+ bFirst = FALSE;
+ }
+
+ return edge;
+}
+
+void put_atom_in_box(matrix box,rvec x)
+{
+ int i,m,d;
+
+ for(m=DIM-1; m>=0; m--) {
+ while (x[m] < 0)
+ for(d=0; d<=m; d++)
+ x[d] += box[m][d];
+ while (x[m] >= box[m][m])
+ for(d=0; d<=m; d++)
+ x[d] -= box[m][d];
+ }
+}
+
+void put_atoms_in_box(matrix box,int natoms,rvec x[])
+{
+ int i,m,d;
+
+ for(i=0; (i<natoms); i++)
+ put_atom_in_box(box,x[i]);
+}
+
+void put_atoms_in_triclinic_unitcell(int ecenter,matrix box,
+ int natoms,rvec x[])
+{
+ rvec box_center,shift_center;
+ real shm01,shm02,shm12,shift;
+ int i,m,d;
+
+ calc_box_center(ecenter,box,box_center);
+
+ /* The product of matrix shm with a coordinate gives the shift vector
+ which is required determine the periodic cell position */
+ shm01 = box[1][0]/box[1][1];
+ shm02 = (box[1][1]*box[2][0] - box[2][1]*box[1][0])/(box[1][1]*box[2][2]);
+ shm12 = box[2][1]/box[2][2];
+
+ clear_rvec(shift_center);
+ for(d=0; d<DIM; d++)
+ rvec_inc(shift_center,box[d]);
+ svmul(0.5,shift_center,shift_center);
+ rvec_sub(box_center,shift_center,shift_center);
+
+ shift_center[0] = shm01*shift_center[1] + shm02*shift_center[2];
+ shift_center[1] = shm12*shift_center[2];
+ shift_center[2] = 0;
+
+ for(i=0; (i<natoms); i++)
+ for(m=DIM-1; m>=0; m--) {
+ shift = shift_center[m];
+ if (m == 0) {
+ shift += shm01*x[i][1] + shm02*x[i][2];
+ } else if (m == 1) {
+ shift += shm12*x[i][2];
+ }
+ while (x[i][m]-shift < 0)
+ for(d=0; d<=m; d++)
+ x[i][d] += box[m][d];
+ while (x[i][m]-shift >= box[m][m])
+ for(d=0; d<=m; d++)
+ x[i][d] -= box[m][d];
+ }
+}
+
+const char *
+put_atoms_in_compact_unitcell(int ePBC,int ecenter,matrix box,
+ int natoms,rvec x[])
+ {
+ t_pbc pbc;
+ rvec box_center,dx;
+ int i;
+
+ set_pbc(&pbc,ePBC,box);
+ calc_box_center(ecenter,box,box_center);
+ for(i=0; i<natoms; i++) {
+ pbc_dx(&pbc,x[i],box_center,dx);
+ rvec_add(box_center,dx,x[i]);
+ }
+
+ return pbc.bLimitDistance ?
+ "WARNING: Could not put all atoms in the compact unitcell\n"
+ : NULL;
+ }
+
--- /dev/null
- if ((nodeid != hid[ai+1]) ||
- (nodeid != hid[ai+2]))
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <assert.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "sysstuff.h"
+#include "macros.h"
+#include "smalloc.h"
+#include "typedefs.h"
+#include "mshift.h"
+#include "invblock.h"
+#include "txtdump.h"
+#include <math.h>
+#include "gmx_fatal.h"
+#include "splitter.h"
+
+typedef struct {
+ int nr;
+ t_iatom *ia;
+} t_sf;
+
+static t_sf *init_sf(int nr)
+{
+ t_sf *sf;
+ int i;
+
+ snew(sf,nr);
+ for(i=0; (i<nr); i++) {
+ sf[i].nr=0;
+ sf[i].ia=NULL;
+ }
+
+ return sf;
+}
+
+static void done_sf(int nr,t_sf *sf)
+{
+ int i;
+
+ for(i=0; (i<nr); i++) {
+ sf[i].nr=0;
+ sfree(sf[i].ia);
+ sf[i].ia=NULL;
+ }
+ sfree(sf);
+}
+
+static void push_sf(t_sf *sf,int nr,t_iatom ia[])
+{
+ int i;
+
+ srenew(sf->ia,sf->nr+nr);
+ for(i=0; (i<nr); i++)
+ sf->ia[sf->nr+i]=ia[i];
+ sf->nr+=nr;
+}
+
+static int min_nodeid(int nr,atom_id list[],int hid[])
+{
+ int i,nodeid,minnodeid;
+
+ if (nr <= 0)
+ gmx_incons("Invalid node number");
+ minnodeid=hid[list[0]];
+ for (i=1; (i<nr); i++)
+ if ((nodeid=hid[list[i]]) < minnodeid)
+ minnodeid=nodeid;
+
+ return minnodeid;
+}
+
+
+static void split_force2(t_inputrec *ir, int nnodes,int hid[],int ftype,t_ilist *ilist,
+ int *multinr,
+ int *constr_min_nodeid,int * constr_max_nodeid,
+ int *left_range, int *right_range)
+{
+ int i,j,k,type,nodeid,nratoms,tnr;
+ int nvsite_constr;
+ t_iatom ai,aj;
+ int node_low_ai,node_low_aj,node_high_ai,node_high_aj;
+ int node_low,node_high;
+ int nodei,nodej;
+ t_sf *sf;
+ int nextra;
+
+ sf = init_sf(nnodes);
+
+ node_high = node_low = 0;
+ nextra = 0;
+
+ /* Walk along all the bonded forces, find the appropriate node
+ * to calc it on, and add it to that nodes list.
+ */
+ for (i=0; i<ilist->nr; i+=(1+nratoms))
+ {
+ type = ilist->iatoms[i];
+ nratoms = interaction_function[ftype].nratoms;
+
+ if (ftype == F_CONSTR)
+ {
+ ai = ilist->iatoms[i+1];
+ aj = ilist->iatoms[i+2];
+
+ nodei = hid[ai];
+ nodej = hid[aj];
+ nodeid = nodei;
+
+ if(ir->eConstrAlg == econtLINCS)
+ {
+ node_low_ai = constr_min_nodeid[ai];
+ node_low_aj = constr_min_nodeid[aj];
+ node_high_ai = constr_max_nodeid[ai];
+ node_high_aj = constr_max_nodeid[aj];
+
+ node_low = min(node_low_ai,node_low_aj);
+ node_high = max(node_high_ai,node_high_aj);
+
+ if (node_high-nodei > 1 || nodei-node_low > 1 ||
+ node_high-nodej > 1 || nodej-node_low > 1 )
+ {
+ gmx_fatal(FARGS,"Constraint dependencies further away than next-neighbor\n"
+ "in particle decomposition. Constraint between atoms %d--%d evaluated\n"
+ "on node %d and %d, but atom %d has connections within %d bonds (lincs_order)\n"
+ "of node %d, and atom %d has connections within %d bonds of node %d.\n"
+ "Reduce the # nodes, lincs_order, or\n"
+ "try domain decomposition.",ai,aj,nodei,nodej,ai,ir->nProjOrder,node_low,aj,ir->nProjOrder,node_high);
+ }
+
+ if (node_low < nodei || node_low < nodej)
+ {
+ right_range[node_low] = max(right_range[node_low],aj);
+ }
+ if (node_high > nodei || node_high > nodej)
+ {
+ left_range[node_high] = min(left_range[node_high],ai);
+ }
+ }
+ else
+ {
+ /* Shake */
+ if (hid[ilist->iatoms[i+2]] != nodei)
+ gmx_fatal(FARGS,"Shake block crossing node boundaries\n"
+ "constraint between atoms (%d,%d) (try LINCS instead!)",
+ ilist->iatoms[i+1]+1,ilist->iatoms[i+2]+1);
+ }
+ }
+ else if (ftype == F_SETTLE)
+ {
+ /* Only the first particle is stored for settles ... */
+ ai=ilist->iatoms[i+1];
+ nodeid=hid[ai];
- "constraint between atoms (%d-%d)",ai+1,ai+2+1);
++ if (nodeid != hid[ilist->iatoms[i+2]] ||
++ nodeid != hid[ilist->iatoms[i+3]])
+ gmx_fatal(FARGS,"Settle block crossing node boundaries\n"
- g0=g->start;
++ "constraint between atoms %d, %d, %d)",
++ ai,ilist->iatoms[i+2],ilist->iatoms[i+3]);
+ }
+ else if(interaction_function[ftype].flags & IF_VSITE)
+ {
+ /* Virtual sites are special, since we need to pre-communicate
+ * their coordinates to construct vsites before then main
+ * coordinate communication.
+ * Vsites can have constructing atoms both larger and smaller than themselves.
+ * To minimize communication and book-keeping, each vsite is constructed on
+ * the home node of the atomnr of the vsite.
+ * Since the vsite coordinates too have to be communicated to the next node,
+ * we need to
+ *
+ * 1. Pre-communicate coordinates of constructing atoms
+ * 2. Construct the vsite
+ * 3. Perform main coordinate communication
+ *
+ * Note that this has change from gromacs 4.0 and earlier, where the vsite
+ * was constructed on the home node of the lowest index of any of the constructing
+ * atoms and the vsite itself.
+ */
+
+ if (ftype==F_VSITE2)
+ nvsite_constr=2;
+ else if(ftype==F_VSITE4FD || ftype==F_VSITE4FDN)
+ nvsite_constr=4;
+ else
+ nvsite_constr=3;
+
+ /* Vsites are constructed on the home node of the actual site to save communication
+ * and simplify the book-keeping.
+ */
+ nodeid=hid[ilist->iatoms[i+1]];
+
+ for(k=2;k<nvsite_constr+2;k++)
+ {
+ if(hid[ilist->iatoms[i+k]]<(nodeid-1) ||
+ hid[ilist->iatoms[i+k]]>(nodeid+1))
+ gmx_fatal(FARGS,"Virtual site %d and its constructing"
+ " atoms are not on the same or adjacent\n"
+ " nodes. This is necessary to avoid a lot\n"
+ " of extra communication. The easiest way"
+ " to ensure this is to place virtual sites\n"
+ " close to the constructing atoms.\n"
+ " Sorry, but you will have to rework your topology!\n",
+ ilist->iatoms[i+1]);
+ }
+ }
+ else
+ {
+ nodeid=min_nodeid(nratoms,&ilist->iatoms[i+1],hid);
+ }
+
+ if (ftype == F_CONSTR && ir->eConstrAlg == econtLINCS)
+ {
+ push_sf(&(sf[nodeid]),nratoms+1,&(ilist->iatoms[i]));
+
+ if(node_low<nodeid)
+ {
+ push_sf(&(sf[node_low]),nratoms+1,&(ilist->iatoms[i]));
+ nextra+=nratoms+1;
+ }
+ if(node_high>nodeid)
+ {
+ push_sf(&(sf[node_high]),nratoms+1,&(ilist->iatoms[i]));
+ nextra+=nratoms+1;
+ }
+ }
+ else
+ {
+ push_sf(&(sf[nodeid]),nratoms+1,&(ilist->iatoms[i]));
+ }
+ }
+
+ if(nextra>0)
+ {
+ ilist->nr += nextra;
+ srenew(ilist->iatoms,ilist->nr);
+ }
+
+ tnr=0;
+ for(nodeid=0; (nodeid<nnodes); nodeid++) {
+ for (i=0; (i<sf[nodeid].nr); i++)
+ ilist->iatoms[tnr++]=sf[nodeid].ia[i];
+
+ multinr[nodeid]=(nodeid==0) ? 0 : multinr[nodeid-1];
+ multinr[nodeid]+=sf[nodeid].nr;
+ }
+
+ if (tnr != ilist->nr)
+ gmx_incons("Splitting forces over processors");
+
+ done_sf(nnodes,sf);
+}
+
+static int *home_index(int nnodes,t_block *cgs,int *multinr)
+{
+ /* This routine determines the node id for each particle */
+ int *hid;
+ int nodeid,j0,j1,j,k;
+
+ snew(hid,cgs->index[cgs->nr]);
+ /* Initiate to -1 to make it possible to check afterwards,
+ * all hid's should be set in the loop below
+ */
+ for(k=0; (k<cgs->index[cgs->nr]); k++)
+ hid[k]=-1;
+
+ /* loop over nodes */
+ for(nodeid=0; (nodeid<nnodes); nodeid++) {
+ j0 = (nodeid==0) ? 0 : multinr[nodeid-1];
+ j1 = multinr[nodeid];
+
+ /* j0 and j1 are the boundariesin the index array */
+ for(j=j0; (j<j1); j++) {
+ for(k=cgs->index[j]; (k<cgs->index[j+1]); k++) {
+ hid[k]=nodeid;
+ }
+ }
+ }
+ /* Now verify that all hid's are not -1 */
+ for(k=0; (k<cgs->index[cgs->nr]); k++)
+ if (hid[k] == -1)
+ gmx_fatal(FARGS,"hid[%d] = -1, cgs->nr = %d, natoms = %d",
+ k,cgs->nr,cgs->index[cgs->nr]);
+
+ return hid;
+}
+
+typedef struct {
+ int atom,ic,is;
+} t_border;
+
+void set_bor(t_border *b,int atom,int ic,int is)
+{
+ if (debug)
+ fprintf(debug,"border @ atom %5d [ ic = %5d, is = %5d ]\n",atom,ic,is);
+ b->atom = atom;
+ b->ic = ic;
+ b->is = is;
+}
+
+static gmx_bool is_bor(atom_id ai[],int i)
+{
+ return ((ai[i] != ai[i-1]) || ((ai[i] == NO_ATID) && (ai[i-1] == NO_ATID)));
+}
+
+static t_border *mk_border(FILE *fp,int natom,atom_id *invcgs,
+ atom_id *invshk,int *nb)
+{
+ t_border *bor;
+ atom_id *sbor,*cbor;
+ int i,j,is,ic,ns,nc,nbor;
+
+ if (debug) {
+ for(i=0; (i<natom); i++) {
+ fprintf(debug,"atom: %6d cgindex: %6d shkindex: %6d\n",
+ i, invcgs[i], invshk[i]);
+ }
+ }
+
+ snew(sbor,natom+1);
+ snew(cbor,natom+1);
+ ns = nc = 1;
+ for(i=1; (i<natom); i++) {
+ if (is_bor(invcgs,i))
+ cbor[nc++] = i;
+ if (is_bor(invshk,i))
+ sbor[ns++] = i;
+ }
+ sbor[ns] = 0;
+ cbor[nc] = 0;
+ if (fp)
+ {
+ fprintf(fp,"There are %d charge group borders",nc);
+ if(invshk!=NULL)
+ {
+ fprintf(fp," and %d shake borders",ns);
+ }
+ fprintf(fp,".\n");
+ }
+ snew(bor,max(nc,ns));
+ ic = is = nbor = 0;
+ while ((ic < nc) || (is < ns)) {
+ if (sbor[is] == cbor[ic]) {
+ set_bor(&(bor[nbor]),cbor[ic],ic,is);
+ nbor++;
+ if (ic < nc) ic++;
+ if (is < ns) is++;
+ }
+ else if (cbor[ic] > sbor[is]) {
+ if (is == ns) {
+ set_bor(&(bor[nbor]),cbor[ic],ic,is);
+ nbor++;
+ if (ic < nc) ic++;
+ }
+ else if (is < ns)
+ is++;
+ }
+ else if (ic < nc)
+ ic++;
+ else
+ is++;/*gmx_fatal(FARGS,"Can't happen is=%d, ic=%d (%s, %d)",
+ is,ic,__FILE__,__LINE__);*/
+ }
+ if (fp)
+ fprintf(fp,"There are %d total borders\n",nbor);
+
+ if (debug) {
+ fprintf(debug,"There are %d actual bor entries\n",nbor);
+ for(i=0; (i<nbor); i++)
+ fprintf(debug,"bor[%5d] = atom: %d ic: %d is: %d\n",i,
+ bor[i].atom,bor[i].ic,bor[i].is);
+ }
+
+ *nb = nbor;
+
+ return bor;
+}
+
+static void split_blocks(FILE *fp,t_inputrec *ir, int nnodes,
+ t_block *cgs,t_blocka *sblock,real capacity[],
+ int *multinr_cgs)
+{
+ int natoms,*maxatom;
+ int i,ii,ai,b0,b1;
+ int nodeid,last_shk,nbor;
+ t_border *border;
+ double tload,tcap;
+
+ gmx_bool bSHK;
+ atom_id *shknum,*cgsnum;
+
+ natoms = cgs->index[cgs->nr];
+
+ if (NULL != debug) {
+ pr_block(debug,0,"cgs",cgs,TRUE);
+ pr_blocka(debug,0,"sblock",sblock,TRUE);
+ fflush(debug);
+ }
+
+ cgsnum = make_invblock(cgs,natoms+1);
+ shknum = make_invblocka(sblock,natoms+1);
+ border = mk_border(fp,natoms,cgsnum,shknum,&nbor);
+
+ snew(maxatom,nnodes);
+ tload = capacity[0]*natoms;
+ tcap = 1.0;
+ nodeid = 0;
+ /* Start at bor is 1, to force the first block on the first processor */
+ for(i=0; (i<nbor) && (tload < natoms); i++) {
+ if(i<(nbor-1))
+ b1=border[i+1].atom;
+ else
+ b1=natoms;
+
+ b0 = border[i].atom;
+
+ if ((fabs(b0-tload)<fabs(b1-tload))) {
+ /* New nodeid time */
+ multinr_cgs[nodeid] = border[i].ic;
+ maxatom[nodeid] = b0;
+ tcap -= capacity[nodeid];
+ nodeid++;
+
+ /* Recompute target load */
+ tload = b0 + (natoms-b0)*capacity[nodeid]/tcap;
+
+ if(debug)
+ printf("tload: %g tcap: %g nodeid: %d\n",tload,tcap,nodeid);
+ }
+ }
+ /* Now the last one... */
+ while (nodeid < nnodes) {
+ multinr_cgs[nodeid] = cgs->nr;
+ /* Store atom number, see above */
+ maxatom[nodeid] = natoms;
+ nodeid++;
+ }
+ if (nodeid != nnodes) {
+ gmx_fatal(FARGS,"nodeid = %d, nnodes = %d, file %s, line %d",
+ nodeid,nnodes,__FILE__,__LINE__);
+ }
+
+ for(i=nnodes-1; (i>0); i--)
+ maxatom[i]-=maxatom[i-1];
+
+ if (fp) {
+ fprintf(fp,"Division over nodes in atoms:\n");
+ for(i=0; (i<nnodes); i++)
+ fprintf(fp," %7d",maxatom[i]);
+ fprintf(fp,"\n");
+ }
+
+ sfree(maxatom);
+ sfree(shknum);
+ sfree(cgsnum);
+ sfree(border);
+}
+
+typedef struct {
+ int atom,sid;
+} t_sid;
+
+static int sid_comp(const void *a,const void *b)
+{
+ t_sid *sa,*sb;
+ int dd;
+
+ sa=(t_sid *)a;
+ sb=(t_sid *)b;
+
+ dd=sa->sid-sb->sid;
+ if (dd == 0)
+ return (sa->atom-sb->atom);
+ else
+ return dd;
+}
+
+static int mk_grey(int nnodes,egCol egc[],t_graph *g,int *AtomI,
+ int maxsid,t_sid sid[])
+{
+ int j,ng,ai,aj,g0;
+
+ ng=0;
+ ai=*AtomI;
+
- g0=g->start;
++ g0=g->at_start;
+ /* Loop over all the bonds */
+ for(j=0; (j<g->nedge[ai]); j++) {
+ aj=g->edge[ai][j]-g0;
+ /* If there is a white one, make it gray and set pbc */
+ if (egc[aj] == egcolWhite) {
+ if (aj < *AtomI)
+ *AtomI = aj;
+ egc[aj] = egcolGrey;
+
+ /* Check whether this one has been set before... */
+ range_check(aj+g0,0,maxsid);
+ range_check(ai+g0,0,maxsid);
+ if (sid[aj+g0].sid != -1)
+ gmx_fatal(FARGS,"sid[%d]=%d, sid[%d]=%d, file %s, line %d",
+ ai,sid[ai+g0].sid,aj,sid[aj+g0].sid,__FILE__,__LINE__);
+ else {
+ sid[aj+g0].sid = sid[ai+g0].sid;
+ sid[aj+g0].atom = aj+g0;
+ }
+ ng++;
+ }
+ }
+ return ng;
+}
+
+static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[])
+/* Return the first node with colour Col starting at fC.
+ * return -1 if none found.
+ */
+{
+ int i;
+
+ for(i=fC; (i<g->nnodes); i++)
+ if ((g->nedge[i] > 0) && (egc[i]==Col))
+ return i;
+
+ return -1;
+}
+
+static int mk_sblocks(FILE *fp,t_graph *g,int maxsid,t_sid sid[])
+{
+ int ng,nnodes;
+ int nW,nG,nB; /* Number of Grey, Black, White */
+ int fW,fG; /* First of each category */
+ egCol *egc=NULL; /* The colour of each node */
+ int g0,nblock;
+
+ if (!g->nbound)
+ return 0;
+
+ nblock=0;
+
+ nnodes=g->nnodes;
+ snew(egc,nnodes);
+
++ g0=g->at_start;
+ nW=g->nbound;
+ nG=0;
+ nB=0;
+
+ fW=0;
+
+ /* We even have a loop invariant:
+ * nW+nG+nB == g->nbound
+ */
+
+ if (fp)
+ fprintf(fp,"Walking down the molecule graph to make constraint-blocks\n");
+
+ while (nW > 0) {
+ /* Find the first white, this will allways be a larger
+ * number than before, because no nodes are made white
+ * in the loop
+ */
+ if ((fW=first_colour(fW,egcolWhite,g,egc)) == -1)
+ gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW);
+
+ /* Make the first white node grey, and set the block number */
+ egc[fW] = egcolGrey;
+ range_check(fW+g0,0,maxsid);
+ sid[fW+g0].sid = nblock++;
+ nG++;
+ nW--;
+
+ /* Initial value for the first grey */
+ fG=fW;
+
+ if (debug)
+ fprintf(debug,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
+ nW,nG,nB,nW+nG+nB);
+
+ while (nG > 0) {
+ if ((fG=first_colour(fG,egcolGrey,g,egc)) == -1)
+ gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG);
+
+ /* Make the first grey node black */
+ egc[fG]=egcolBlack;
+ nB++;
+ nG--;
+
+ /* Make all the neighbours of this black node grey
+ * and set their block number
+ */
+ ng=mk_grey(nnodes,egc,g,&fG,maxsid,sid);
+ /* ng is the number of white nodes made grey */
+ nG+=ng;
+ nW-=ng;
+ }
+ }
+ sfree(egc);
+
+ if (debug)
+ fprintf(debug,"Found %d shake blocks\n",nblock);
+
+ return nblock;
+}
+
+
+typedef struct {
+ int first,last,sid;
+} t_merge_sid;
+
+static int ms_comp(const void *a, const void *b)
+{
+ t_merge_sid *ma = (t_merge_sid *)a;
+ t_merge_sid *mb = (t_merge_sid *)b;
+ int d;
+
+ d = ma->first-mb->first;
+ if (d == 0)
+ return ma->last-mb->last;
+ else
+ return d;
+}
+
+static int merge_sid(int i0,int at_start,int at_end,int nsid,t_sid sid[],
+ t_blocka *sblock)
+{
+ int i,j,k,n,isid,ndel;
+ t_merge_sid *ms;
+ int nChanged;
+
+ /* We try to remdy the following problem:
+ * Atom: 1 2 3 4 5 6 7 8 9 10
+ * Sid: 0 -1 1 0 -1 1 1 1 1 1
+ */
+
+ /* Determine first and last atom in each shake ID */
+ snew(ms,nsid);
+
+ for(k=0; (k<nsid); k++) {
+ ms[k].first = at_end+1;
+ ms[k].last = -1;
+ ms[k].sid = k;
+ }
+ for(i=at_start; (i<at_end); i++) {
+ isid = sid[i].sid;
+ range_check(isid,-1,nsid);
+ if (isid >= 0) {
+ ms[isid].first = min(ms[isid].first,sid[i].atom);
+ ms[isid].last = max(ms[isid].last,sid[i].atom);
+ }
+ }
+ qsort(ms,nsid,sizeof(ms[0]),ms_comp);
+
+ /* Now merge the overlapping ones */
+ ndel = 0;
+ for(k=0; (k<nsid); ) {
+ for(j=k+1; (j<nsid); ) {
+ if (ms[j].first <= ms[k].last) {
+ ms[k].last = max(ms[k].last,ms[j].last);
+ ms[k].first = min(ms[k].first,ms[j].first);
+ ms[j].sid = -1;
+ ndel++;
+ j++;
+ }
+ else {
+ k = j;
+ j = k+1;
+ }
+ }
+ if (j == nsid)
+ k++;
+ }
+ for(k=0; (k<nsid); k++) {
+ while ((k < nsid-1) && (ms[k].sid == -1)) {
+ for(j=k+1; (j<nsid); j++) {
+ memcpy(&(ms[j-1]),&(ms[j]),sizeof(ms[0]));
+ }
+ nsid--;
+ }
+ }
+
+ for(k=at_start; (k<at_end); k++) {
+ sid[k].atom = k;
+ sid[k].sid = -1;
+ }
+ sblock->nr = nsid;
+ sblock->nalloc_index = sblock->nr+1;
+ snew(sblock->index,sblock->nalloc_index);
+ sblock->nra = at_end - at_start;
+ sblock->nalloc_a = sblock->nra;
+ snew(sblock->a,sblock->nalloc_a);
+ sblock->index[0] = 0;
+ for(k=n=0; (k<nsid); k++) {
+ sblock->index[k+1] = sblock->index[k] + ms[k].last - ms[k].first+1;
+ for(j=ms[k].first; (j<=ms[k].last); j++) {
+ range_check(n,0,sblock->nra);
+ sblock->a[n++] = j;
+ range_check(j,0,at_end);
+ if (sid[j].sid == -1)
+ sid[j].sid = k;
+ else
+ fprintf(stderr,"Double sids (%d, %d) for atom %d\n",sid[j].sid,k,j);
+ }
+ }
+ assert(k == nsid);
+ /* Removed 2007-09-04
+ sblock->index[k+1] = natoms;
+ for(k=0; (k<natoms); k++)
+ if (sid[k].sid == -1)
+ sblock->a[n++] = k;
+ assert(n == natoms);
+ */
+ sblock->nra = n;
+ assert(sblock->index[k] == sblock->nra);
+ sfree(ms);
+
+ return nsid;
+}
+
+void gen_sblocks(FILE *fp,int at_start,int at_end,
+ t_idef *idef,t_blocka *sblock,
+ gmx_bool bSettle)
+{
+ t_graph *g;
+ int i,i0,j,k,istart,n;
+ t_sid *sid;
+ int isid,nsid;
+
+ g=mk_graph(NULL,idef,at_start,at_end,TRUE,bSettle);
+ if (debug)
+ p_graph(debug,"Graaf Dracula",g);
+ snew(sid,at_end);
+ for(i=at_start; (i<at_end); i++) {
+ sid[i].atom = i;
+ sid[i].sid = -1;
+ }
+ nsid=mk_sblocks(fp,g,at_end,sid);
+
+ if (!nsid)
+ return;
+
+ /* Now sort the shake blocks... */
+ qsort(sid+at_start,at_end-at_start,(size_t)sizeof(sid[0]),sid_comp);
+
+ if (debug) {
+ fprintf(debug,"Sorted shake block\n");
+ for(i=at_start; (i<at_end); i++)
+ fprintf(debug,"sid[%5d] = atom:%5d sid:%5d\n",i,sid[i].atom,sid[i].sid);
+ }
+ /* Now check how many are NOT -1, i.e. how many have to be shaken */
+ for(i0=at_start; (i0<at_end); i0++)
+ if (sid[i0].sid > -1)
+ break;
+
+ /* Now we have the sids that have to be shaken. We'll check the min and
+ * max atom numbers and this determines the shake block. DvdS 2007-07-19.
+ * For the purpose of making boundaries all atoms in between need to be
+ * part of the shake block too. There may be cases where blocks overlap
+ * and they will have to be merged.
+ */
+ nsid = merge_sid(i0,at_start,at_end,nsid,sid,sblock);
+ /* Now sort the shake blocks again... */
+ /*qsort(sid,natoms,(size_t)sizeof(sid[0]),sid_comp);*/
+
+ /* Fill the sblock struct */
+ /* sblock->nr = nsid;
+ sblock->nra = natoms;
+ srenew(sblock->a,sblock->nra);
+ srenew(sblock->index,sblock->nr+1);
+
+ i = i0;
+ isid = sid[i].sid;
+ n = k = 0;
+ sblock->index[n++]=k;
+ while (i < natoms) {
+ istart = sid[i].atom;
+ while ((i<natoms-1) && (sid[i+1].sid == isid))
+ i++;*/
+ /* After while: we found a new block, or are thru with the atoms */
+ /* for(j=istart; (j<=sid[i].atom); j++,k++)
+ sblock->a[k]=j;
+ sblock->index[n] = k;
+ if (i < natoms-1)
+ n++;
+ if (n > nsid)
+ gmx_fatal(FARGS,"Death Horror: nsid = %d, n= %d",nsid,n);
+ i++;
+ isid = sid[i].sid;
+ }
+ */
+ sfree(sid);
+ /* Due to unknown reason this free generates a problem sometimes */
+ done_graph(g);
+ sfree(g);
+ if (debug)
+ fprintf(debug,"Done gen_sblocks\n");
+}
+
+static t_blocka block2blocka(t_block *block)
+{
+ t_blocka blocka;
+ int i;
+
+ blocka.nr = block->nr;
+ blocka.nalloc_index = blocka.nr + 1;
+ snew(blocka.index,blocka.nalloc_index);
+ for(i=0; i<=block->nr; i++)
+ blocka.index[i] = block->index[i];
+ blocka.nra = block->index[block->nr];
+ blocka.nalloc_a = blocka.nra;
+ snew(blocka.a,blocka.nalloc_a);
+ for(i=0; i<blocka.nra; i++)
+ blocka.a[i] = i;
+
+ return blocka;
+}
+
+typedef struct
+{
+ int nconstr;
+ int index[10];
+} pd_constraintlist_t;
+
+
+static void
+find_constraint_range_recursive(pd_constraintlist_t * constraintlist,
+ int thisatom,
+ int depth,
+ int * min_atomid,
+ int * max_atomid)
+{
+ int i,j;
+ int nconstr;
+
+ for(i=0;i<constraintlist[thisatom].nconstr;i++)
+ {
+ j = constraintlist[thisatom].index[i];
+
+ *min_atomid = (j<*min_atomid) ? j : *min_atomid;
+ *max_atomid = (j>*max_atomid) ? j : *max_atomid;
+
+ if(depth>0)
+ {
+ find_constraint_range_recursive(constraintlist,j,depth-1,min_atomid,max_atomid);
+ }
+ }
+}
+
+static void
+pd_determine_constraints_range(t_inputrec * ir,
+ int natoms,
+ t_ilist * ilist,
+ int hid[],
+ int * min_nodeid,
+ int * max_nodeid)
+{
+ int i,j,k;
+ int nratoms;
+ int depth;
+ int ai,aj;
+ int min_atomid,max_atomid;
+ pd_constraintlist_t *constraintlist;
+
+ nratoms = interaction_function[F_CONSTR].nratoms;
+ depth = ir->nProjOrder;
+
+ snew(constraintlist,natoms);
+
+ /* Make a list of all the connections */
+ for(i=0;i<ilist->nr;i+=nratoms+1)
+ {
+ ai=ilist->iatoms[i+1];
+ aj=ilist->iatoms[i+2];
+ constraintlist[ai].index[constraintlist[ai].nconstr++]=aj;
+ constraintlist[aj].index[constraintlist[aj].nconstr++]=ai;
+ }
+
+ for(i=0;i<natoms;i++)
+ {
+ min_atomid = i;
+ max_atomid = i;
+
+ find_constraint_range_recursive(constraintlist,i,depth,&min_atomid,&max_atomid);
+
+ min_nodeid[i] = hid[min_atomid];
+ max_nodeid[i] = hid[max_atomid];
+ }
+ sfree(constraintlist);
+}
+
+
+void split_top(FILE *fp,int nnodes,gmx_localtop_t *top,t_inputrec *ir,t_block *mols,
+ real *capacity,int *multinr_cgs,int **multinr_nre, int *left_range, int * right_range)
+{
+ int natoms,i,j,k,mj,atom,maxatom,sstart,send,bstart,nodeid;
+ t_blocka sblock;
+ int *homeind;
+ int ftype,nvsite_constr,nra,nrd;
+ t_iatom *ia;
+ int minhome,ihome,minidx;
+ int *constr_min_nodeid;
+ int *constr_max_nodeid;
+
+ if (nnodes <= 1)
+ return;
+
+ natoms = mols->index[mols->nr];
+
+ if (fp)
+ fprintf(fp,"splitting topology...\n");
+
+/*#define MOL_BORDER*/
+/*Removed the above to allow splitting molecules with h-bond constraints
+ over processors. The results in DP are the same. */
+ init_blocka(&sblock);
+ if(ir->eConstrAlg != econtLINCS)
+ {
+#ifndef MOL_BORDER
+ /* Make a special shake block that includes settles */
+ gen_sblocks(fp,0,natoms,&top->idef,&sblock,TRUE);
+#else
+ sblock = block2blocka(mols);
+#endif
+ }
+
+ split_blocks(fp,ir,nnodes,&top->cgs,&sblock,capacity,multinr_cgs);
+
+ homeind=home_index(nnodes,&top->cgs,multinr_cgs);
+
+ snew(constr_min_nodeid,natoms);
+ snew(constr_max_nodeid,natoms);
+
+ if(top->idef.il[F_CONSTR].nr>0)
+ {
+ pd_determine_constraints_range(ir,natoms,&top->idef.il[F_CONSTR],homeind,constr_min_nodeid,constr_max_nodeid);
+ }
+ else
+ {
+ /* Not 100% necessary, but it is a bad habit to have uninitialized arrays around... */
+ for(i=0;i<natoms;i++)
+ {
+ constr_min_nodeid[i] = constr_max_nodeid[i] = homeind[i];
+ }
+ }
+
+ /* Default limits (no communication) for PD constraints */
+ left_range[0] = 0;
+ for(i=1;i<nnodes;i++)
+ {
+ left_range[i] = top->cgs.index[multinr_cgs[i-1]];
+ right_range[i-1] = left_range[i]-1;
+ }
+ right_range[nnodes-1] = top->cgs.index[multinr_cgs[nnodes-1]]-1;
+
+ for(j=0; (j<F_NRE); j++)
+ split_force2(ir, nnodes,homeind,j,&top->idef.il[j],multinr_nre[j],constr_min_nodeid,constr_max_nodeid,
+ left_range,right_range);
+
+ sfree(constr_min_nodeid);
+ sfree(constr_max_nodeid);
+
+ sfree(homeind);
+ done_blocka(&sblock);
+}
+
--- /dev/null
- nsettle = settle->nr/2;
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * GROwing Monsters And Cloning Shrimps
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "confio.h"
+#include "constr.h"
+#include "copyrite.h"
+#include "invblock.h"
+#include "main.h"
+#include "mdrun.h"
+#include "nrnb.h"
+#include "smalloc.h"
+#include "vec.h"
+#include "physics.h"
+#include "names.h"
+#include "txtdump.h"
+#include "domdec.h"
+#include "pdbio.h"
+#include "partdec.h"
+#include "splitter.h"
+#include "mtop_util.h"
+#include "gmxfio.h"
+#include "macros.h"
+
+typedef struct gmx_constr {
+ int ncon_tot; /* The total number of constraints */
+ int nflexcon; /* The number of flexible constraints */
+ int n_at2con_mt; /* The size of at2con = #moltypes */
+ t_blocka *at2con_mt; /* A list of atoms to constraints */
+ gmx_lincsdata_t lincsd; /* LINCS data */
+ gmx_shakedata_t shaked; /* SHAKE data */
+ gmx_settledata_t settled; /* SETTLE data */
+ int nblocks; /* The number of SHAKE blocks */
+ int *sblock; /* The SHAKE blocks */
+ int sblock_nalloc;/* The allocation size of sblock */
+ real *lagr; /* Lagrange multipliers for SHAKE */
+ int lagr_nalloc; /* The allocation size of lagr */
+ int maxwarn; /* The maximum number of warnings */
+ int warncount_lincs;
+ int warncount_settle;
+ gmx_edsam_t ed; /* The essential dynamics data */
+
+ gmx_mtop_t *warn_mtop; /* Only used for printing warnings */
+} t_gmx_constr;
+
+typedef struct {
+ atom_id iatom[3];
+ atom_id blocknr;
+} t_sortblock;
+
+static void *init_vetavars(t_vetavars *vars,
+ gmx_bool constr_deriv,
+ real veta,real vetanew, t_inputrec *ir, gmx_ekindata_t *ekind, gmx_bool bPscal)
+{
+ double g;
+ int i;
+
+ /* first, set the alpha integrator variable */
+ if ((ir->opts.nrdf[0] > 0) && bPscal)
+ {
+ vars->alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
+ } else {
+ vars->alpha = 1.0;
+ }
+ g = 0.5*veta*ir->delta_t;
+ vars->rscale = exp(g)*series_sinhx(g);
+ g = -0.25*vars->alpha*veta*ir->delta_t;
+ vars->vscale = exp(g)*series_sinhx(g);
+ vars->rvscale = vars->vscale*vars->rscale;
+ vars->veta = vetanew;
+
+ if (constr_deriv)
+ {
+ snew(vars->vscale_nhc,ir->opts.ngtc);
+ if ((ekind==NULL) || (!bPscal))
+ {
+ for (i=0;i<ir->opts.ngtc;i++)
+ {
+ vars->vscale_nhc[i] = 1;
+ }
+ }
+ else
+ {
+ for (i=0;i<ir->opts.ngtc;i++)
+ {
+ vars->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
+ }
+ }
+ }
+ else
+ {
+ vars->vscale_nhc = NULL;
+ }
+
+ return vars;
+}
+
+static void free_vetavars(t_vetavars *vars)
+{
+ if (vars->vscale_nhc != NULL)
+ {
+ sfree(vars->vscale_nhc);
+ }
+}
+
+static int pcomp(const void *p1, const void *p2)
+{
+ int db;
+ atom_id min1,min2,max1,max2;
+ t_sortblock *a1=(t_sortblock *)p1;
+ t_sortblock *a2=(t_sortblock *)p2;
+
+ db=a1->blocknr-a2->blocknr;
+
+ if (db != 0)
+ return db;
+
+ min1=min(a1->iatom[1],a1->iatom[2]);
+ max1=max(a1->iatom[1],a1->iatom[2]);
+ min2=min(a2->iatom[1],a2->iatom[2]);
+ max2=max(a2->iatom[1],a2->iatom[2]);
+
+ if (min1 == min2)
+ return max1-max2;
+ else
+ return min1-min2;
+}
+
+static int icomp(const void *p1, const void *p2)
+{
+ atom_id *a1=(atom_id *)p1;
+ atom_id *a2=(atom_id *)p2;
+
+ return (*a1)-(*a2);
+}
+
+int n_flexible_constraints(struct gmx_constr *constr)
+{
+ int nflexcon;
+
+ if (constr)
+ nflexcon = constr->nflexcon;
+ else
+ nflexcon = 0;
+
+ return nflexcon;
+}
+
+void too_many_constraint_warnings(int eConstrAlg,int warncount)
+{
+ const char *abort="- aborting to avoid logfile runaway.\n"
+ "This normally happens when your system is not sufficiently equilibrated,"
+ "or if you are changing lambda too fast in free energy simulations.\n";
+
+ gmx_fatal(FARGS,
+ "Too many %s warnings (%d)\n"
+ "If you know what you are doing you can %s"
+ "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
+ "but normally it is better to fix the problem",
+ (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE",warncount,
+ (eConstrAlg == econtLINCS) ?
+ "adjust the lincs warning threshold in your mdp file\nor " : "\n");
+}
+
+static void write_constr_pdb(const char *fn,const char *title,
+ gmx_mtop_t *mtop,
+ int start,int homenr,t_commrec *cr,
+ rvec x[],matrix box)
+{
+ char fname[STRLEN],format[STRLEN];
+ FILE *out;
+ int dd_ac0=0,dd_ac1=0,i,ii,resnr;
+ gmx_domdec_t *dd;
+ char *anm,*resnm;
+
+ dd = NULL;
+ if (PAR(cr))
+ {
+ sprintf(fname,"%s_n%d.pdb",fn,cr->sim_nodeid);
+ if (DOMAINDECOMP(cr))
+ {
+ dd = cr->dd;
+ dd_get_constraint_range(dd,&dd_ac0,&dd_ac1);
+ start = 0;
+ homenr = dd_ac1;
+ }
+ }
+ else
+ {
+ sprintf(fname,"%s.pdb",fn);
+ }
+ sprintf(format,"%s\n",get_pdbformat());
+
+ out = gmx_fio_fopen(fname,"w");
+
+ fprintf(out,"TITLE %s\n",title);
+ gmx_write_pdb_box(out,-1,box);
+ for(i=start; i<start+homenr; i++)
+ {
+ if (dd != NULL)
+ {
+ if (i >= dd->nat_home && i < dd_ac0)
+ {
+ continue;
+ }
+ ii = dd->gatindex[i];
+ }
+ else
+ {
+ ii = i;
+ }
+ gmx_mtop_atominfo_global(mtop,ii,&anm,&resnr,&resnm);
+ fprintf(out,format,"ATOM",(ii+1)%100000,
+ anm,resnm,' ',resnr%10000,' ',
+ 10*x[i][XX],10*x[i][YY],10*x[i][ZZ]);
+ }
+ fprintf(out,"TER\n");
+
+ gmx_fio_fclose(out);
+}
+
+static void dump_confs(FILE *fplog,gmx_large_int_t step,gmx_mtop_t *mtop,
+ int start,int homenr,t_commrec *cr,
+ rvec x[],rvec xprime[],matrix box)
+{
+ char buf[256],buf2[22];
+
+ char *env=getenv("GMX_SUPPRESS_DUMP");
+ if (env)
+ return;
+
+ sprintf(buf,"step%sb",gmx_step_str(step,buf2));
+ write_constr_pdb(buf,"initial coordinates",
+ mtop,start,homenr,cr,x,box);
+ sprintf(buf,"step%sc",gmx_step_str(step,buf2));
+ write_constr_pdb(buf,"coordinates after constraining",
+ mtop,start,homenr,cr,xprime,box);
+ if (fplog)
+ {
+ fprintf(fplog,"Wrote pdb files with previous and current coordinates\n");
+ }
+ fprintf(stderr,"Wrote pdb files with previous and current coordinates\n");
+}
+
+static void pr_sortblock(FILE *fp,const char *title,int nsb,t_sortblock sb[])
+{
+ int i;
+
+ fprintf(fp,"%s\n",title);
+ for(i=0; (i<nsb); i++)
+ fprintf(fp,"i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
+ i,sb[i].iatom[0],sb[i].iatom[1],sb[i].iatom[2],
+ sb[i].blocknr);
+}
+
+gmx_bool constrain(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
+ struct gmx_constr *constr,
+ t_idef *idef,t_inputrec *ir,gmx_ekindata_t *ekind,
+ t_commrec *cr,
+ gmx_large_int_t step,int delta_step,
+ t_mdatoms *md,
+ rvec *x,rvec *xprime,rvec *min_proj,matrix box,
+ real lambda,real *dvdlambda,
+ rvec *v,tensor *vir,
+ t_nrnb *nrnb,int econq,gmx_bool bPscal,real veta, real vetanew)
+{
+ gmx_bool bOK,bDump;
+ int start,homenr,nrend;
+ int i,j,d;
+ int ncons,error;
+ tensor rmdr;
+ rvec *vstor;
+ real invdt,vir_fac,t;
+ t_ilist *settle;
+ int nsettle;
+ t_pbc pbc;
+ char buf[22];
+ t_vetavars vetavar;
+
+ if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI))
+ {
+ gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted");
+ }
+
+ bOK = TRUE;
+ bDump = FALSE;
+
+ start = md->start;
+ homenr = md->homenr;
+ nrend = start+homenr;
+
+ /* set constants for pressure control integration */
+ init_vetavars(&vetavar,econq!=econqCoord,
+ veta,vetanew,ir,ekind,bPscal);
+
+ if (ir->delta_t == 0)
+ {
+ invdt = 0;
+ }
+ else
+ {
+ invdt = 1/ir->delta_t;
+ }
+
+ if (ir->efep != efepNO && EI_DYNAMICS(ir->eI))
+ {
+ /* Set the constraint lengths for the step at which this configuration
+ * is meant to be. The invmasses should not be changed.
+ */
+ lambda += delta_step*ir->delta_lambda;
+ }
+
+ if (vir != NULL)
+ {
+ clear_mat(rmdr);
+ }
+
+ where();
+ if (constr->lincsd)
+ {
+ bOK = constrain_lincs(fplog,bLog,bEner,ir,step,constr->lincsd,md,cr,
+ x,xprime,min_proj,box,lambda,dvdlambda,
+ invdt,v,vir!=NULL,rmdr,
+ econq,nrnb,
+ constr->maxwarn,&constr->warncount_lincs);
+ if (!bOK && constr->maxwarn >= 0)
+ {
+ if (fplog != NULL)
+ {
+ fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
+ econstr_names[econtLINCS],gmx_step_str(step,buf));
+ }
+ bDump = TRUE;
+ }
+ }
+
+ if (constr->nblocks > 0)
+ {
+ switch (econq) {
+ case (econqCoord):
+ bOK = bshakef(fplog,constr->shaked,
+ homenr,md->invmass,constr->nblocks,constr->sblock,
+ idef,ir,box,x,xprime,nrnb,
+ constr->lagr,lambda,dvdlambda,
+ invdt,v,vir!=NULL,rmdr,constr->maxwarn>=0,econq,
+ &vetavar);
+ break;
+ case (econqVeloc):
+ bOK = bshakef(fplog,constr->shaked,
+ homenr,md->invmass,constr->nblocks,constr->sblock,
+ idef,ir,box,x,min_proj,nrnb,
+ constr->lagr,lambda,dvdlambda,
+ invdt,NULL,vir!=NULL,rmdr,constr->maxwarn>=0,econq,
+ &vetavar);
+ break;
+ default:
+ gmx_fatal(FARGS,"Internal error, SHAKE called for constraining something else than coordinates");
+ break;
+ }
+
+ if (!bOK && constr->maxwarn >= 0)
+ {
+ if (fplog != NULL)
+ {
+ fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
+ econstr_names[econtSHAKE],gmx_step_str(step,buf));
+ }
+ bDump = TRUE;
+ }
+ }
+
+ settle = &idef->il[F_SETTLE];
+ if (settle->nr > 0)
+ {
- step,ddglatnr(cr->dd,settle->iatoms[error*2+1]));
++ nsettle = settle->nr/4;
+
+ switch (econq)
+ {
+ case econqCoord:
+ csettle(constr->settled,
+ nsettle,settle->iatoms,x[0],xprime[0],
+ invdt,v[0],vir!=NULL,rmdr,&error,&vetavar);
+ inc_nrnb(nrnb,eNR_SETTLE,nsettle);
+ if (v != NULL)
+ {
+ inc_nrnb(nrnb,eNR_CONSTR_V,nsettle*3);
+ }
+ if (vir != NULL)
+ {
+ inc_nrnb(nrnb,eNR_CONSTR_VIR,nsettle*3);
+ }
+
+ bOK = (error < 0);
+ if (!bOK && constr->maxwarn >= 0)
+ {
+ char buf[256];
+ sprintf(buf,
+ "\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
+ "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
- iH = settle->iatoms[1]+1;
++ step,ddglatnr(cr->dd,settle->iatoms[error*4+1]));
+ if (fplog)
+ {
+ fprintf(fplog,"%s",buf);
+ }
+ fprintf(stderr,"%s",buf);
+ constr->warncount_settle++;
+ if (constr->warncount_settle > constr->maxwarn)
+ {
+ too_many_constraint_warnings(-1,constr->warncount_settle);
+ }
+ bDump = TRUE;
+ break;
+ case econqVeloc:
+ case econqDeriv:
+ case econqForce:
+ case econqForceDispl:
+ settle_proj(fplog,constr->settled,econq,
+ nsettle,settle->iatoms,x,
+ xprime,min_proj,vir!=NULL,rmdr,&vetavar);
+ /* This is an overestimate */
+ inc_nrnb(nrnb,eNR_SETTLE,nsettle);
+ break;
+ case econqDeriv_FlexCon:
+ /* Nothing to do, since the are no flexible constraints in settles */
+ break;
+ default:
+ gmx_incons("Unknown constraint quantity for settle");
+ }
+ }
+ }
+
+ free_vetavars(&vetavar);
+
+ if (vir != NULL)
+ {
+ switch (econq)
+ {
+ case econqCoord:
+ vir_fac = 0.5/(ir->delta_t*ir->delta_t);
+ break;
+ case econqVeloc:
+ vir_fac = 0.5/ir->delta_t;
+ break;
+ case econqForce:
+ case econqForceDispl:
+ vir_fac = 0.5;
+ break;
+ default:
+ vir_fac = 0;
+ gmx_incons("Unsupported constraint quantity for virial");
+ }
+
+ if (EI_VV(ir->eI))
+ {
+ vir_fac *= 2; /* only constraining over half the distance here */
+ }
+ for(i=0; i<DIM; i++)
+ {
+ for(j=0; j<DIM; j++)
+ {
+ (*vir)[i][j] = vir_fac*rmdr[i][j];
+ }
+ }
+ }
+
+ if (bDump)
+ {
+ dump_confs(fplog,step,constr->warn_mtop,start,homenr,cr,x,xprime,box);
+ }
+
+ if (econq == econqCoord)
+ {
+ if (ir->ePull == epullCONSTRAINT)
+ {
+ if (EI_DYNAMICS(ir->eI))
+ {
+ t = ir->init_t + (step + delta_step)*ir->delta_t;
+ }
+ else
+ {
+ t = ir->init_t;
+ }
+ set_pbc(&pbc,ir->ePBC,box);
+ pull_constraint(ir->pull,md,&pbc,cr,ir->delta_t,t,x,xprime,v,*vir);
+ }
+ if (constr->ed && delta_step > 0)
+ {
+ /* apply the essential dynamcs constraints here */
+ do_edsam(ir,step,md,cr,xprime,v,box,constr->ed);
+ }
+ }
+
+ return bOK;
+}
+
+real *constr_rmsd_data(struct gmx_constr *constr)
+{
+ if (constr->lincsd)
+ return lincs_rmsd_data(constr->lincsd);
+ else
+ return NULL;
+}
+
+real constr_rmsd(struct gmx_constr *constr,gmx_bool bSD2)
+{
+ if (constr->lincsd)
+ return lincs_rmsd(constr->lincsd,bSD2);
+ else
+ return 0;
+}
+
+static void make_shake_sblock_pd(struct gmx_constr *constr,
+ t_idef *idef,t_mdatoms *md)
+{
+ int i,j,m,ncons;
+ int bstart,bnr;
+ t_blocka sblocks;
+ t_sortblock *sb;
+ t_iatom *iatom;
+ atom_id *inv_sblock;
+
+ /* Since we are processing the local topology,
+ * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
+ */
+ ncons = idef->il[F_CONSTR].nr/3;
+
+ init_blocka(&sblocks);
+ gen_sblocks(NULL,md->start,md->start+md->homenr,idef,&sblocks,FALSE);
+
+ /*
+ bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
+ nblocks=blocks->multinr[idef->nodeid] - bstart;
+ */
+ bstart = 0;
+ constr->nblocks = sblocks.nr;
+ if (debug)
+ fprintf(debug,"ncons: %d, bstart: %d, nblocks: %d\n",
+ ncons,bstart,constr->nblocks);
+
+ /* Calculate block number for each atom */
+ inv_sblock = make_invblocka(&sblocks,md->nr);
+
+ done_blocka(&sblocks);
+
+ /* Store the block number in temp array and
+ * sort the constraints in order of the sblock number
+ * and the atom numbers, really sorting a segment of the array!
+ */
+#ifdef DEBUGIDEF
+ pr_idef(fplog,0,"Before Sort",idef);
+#endif
+ iatom=idef->il[F_CONSTR].iatoms;
+ snew(sb,ncons);
+ for(i=0; (i<ncons); i++,iatom+=3) {
+ for(m=0; (m<3); m++)
+ sb[i].iatom[m] = iatom[m];
+ sb[i].blocknr = inv_sblock[iatom[1]];
+ }
+
+ /* Now sort the blocks */
+ if (debug) {
+ pr_sortblock(debug,"Before sorting",ncons,sb);
+ fprintf(debug,"Going to sort constraints\n");
+ }
+
+ qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
+
+ if (debug) {
+ pr_sortblock(debug,"After sorting",ncons,sb);
+ }
+
+ iatom=idef->il[F_CONSTR].iatoms;
+ for(i=0; (i<ncons); i++,iatom+=3)
+ for(m=0; (m<3); m++)
+ iatom[m]=sb[i].iatom[m];
+#ifdef DEBUGIDEF
+ pr_idef(fplog,0,"After Sort",idef);
+#endif
+
+ j=0;
+ snew(constr->sblock,constr->nblocks+1);
+ bnr=-2;
+ for(i=0; (i<ncons); i++) {
+ if (sb[i].blocknr != bnr) {
+ bnr=sb[i].blocknr;
+ constr->sblock[j++]=3*i;
+ }
+ }
+ /* Last block... */
+ constr->sblock[j++] = 3*ncons;
+
+ if (j != (constr->nblocks+1)) {
+ fprintf(stderr,"bstart: %d\n",bstart);
+ fprintf(stderr,"j: %d, nblocks: %d, ncons: %d\n",
+ j,constr->nblocks,ncons);
+ for(i=0; (i<ncons); i++)
+ fprintf(stderr,"i: %5d sb[i].blocknr: %5u\n",i,sb[i].blocknr);
+ for(j=0; (j<=constr->nblocks); j++)
+ fprintf(stderr,"sblock[%3d]=%5d\n",j,(int)constr->sblock[j]);
+ gmx_fatal(FARGS,"DEATH HORROR: "
+ "sblocks does not match idef->il[F_CONSTR]");
+ }
+ sfree(sb);
+ sfree(inv_sblock);
+}
+
+static void make_shake_sblock_dd(struct gmx_constr *constr,
+ t_ilist *ilcon,t_block *cgs,
+ gmx_domdec_t *dd)
+{
+ int ncons,c,cg;
+ t_iatom *iatom;
+
+ if (dd->ncg_home+1 > constr->sblock_nalloc) {
+ constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
+ srenew(constr->sblock,constr->sblock_nalloc);
+ }
+
+ ncons = ilcon->nr/3;
+ iatom = ilcon->iatoms;
+ constr->nblocks = 0;
+ cg = 0;
+ for(c=0; c<ncons; c++) {
+ if (c == 0 || iatom[1] >= cgs->index[cg+1]) {
+ constr->sblock[constr->nblocks++] = 3*c;
+ while (iatom[1] >= cgs->index[cg+1])
+ cg++;
+ }
+ iatom += 3;
+ }
+ constr->sblock[constr->nblocks] = 3*ncons;
+}
+
+t_blocka make_at2con(int start,int natoms,
+ t_ilist *ilist,t_iparams *iparams,
+ gmx_bool bDynamics,int *nflexiblecons)
+{
+ int *count,ncon,con,con_tot,nflexcon,ftype,i,a;
+ t_iatom *ia;
+ t_blocka at2con;
+ gmx_bool bFlexCon;
+
+ snew(count,natoms);
+ nflexcon = 0;
+ for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
+ ncon = ilist[ftype].nr/3;
+ ia = ilist[ftype].iatoms;
+ for(con=0; con<ncon; con++) {
+ bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
+ iparams[ia[0]].constr.dB == 0);
+ if (bFlexCon)
+ nflexcon++;
+ if (bDynamics || !bFlexCon) {
+ for(i=1; i<3; i++) {
+ a = ia[i] - start;
+ count[a]++;
+ }
+ }
+ ia += 3;
+ }
+ }
+ *nflexiblecons = nflexcon;
+
+ at2con.nr = natoms;
+ at2con.nalloc_index = at2con.nr+1;
+ snew(at2con.index,at2con.nalloc_index);
+ at2con.index[0] = 0;
+ for(a=0; a<natoms; a++) {
+ at2con.index[a+1] = at2con.index[a] + count[a];
+ count[a] = 0;
+ }
+ at2con.nra = at2con.index[natoms];
+ at2con.nalloc_a = at2con.nra;
+ snew(at2con.a,at2con.nalloc_a);
+
+ /* The F_CONSTRNC constraints have constraint numbers
+ * that continue after the last F_CONSTR constraint.
+ */
+ con_tot = 0;
+ for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
+ ncon = ilist[ftype].nr/3;
+ ia = ilist[ftype].iatoms;
+ for(con=0; con<ncon; con++) {
+ bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
+ iparams[ia[0]].constr.dB == 0);
+ if (bDynamics || !bFlexCon) {
+ for(i=1; i<3; i++) {
+ a = ia[i] - start;
+ at2con.a[at2con.index[a]+count[a]++] = con_tot;
+ }
+ }
+ con_tot++;
+ ia += 3;
+ }
+ }
+
+ sfree(count);
+
+ return at2con;
+}
+
+void set_constraints(struct gmx_constr *constr,
+ gmx_localtop_t *top,t_inputrec *ir,
+ t_mdatoms *md,t_commrec *cr)
+{
+ t_idef *idef;
+ int ncons;
+ t_ilist *settle;
+ int iO,iH;
+
+ idef = &top->idef;
+
+ if (constr->ncon_tot > 0)
+ {
+ /* We are using the local topology,
+ * so there are only F_CONSTR constraints.
+ */
+ ncons = idef->il[F_CONSTR].nr/3;
+
+ /* With DD we might also need to call LINCS with ncons=0 for
+ * communicating coordinates to other nodes that do have constraints.
+ */
+ if (ir->eConstrAlg == econtLINCS)
+ {
+ set_lincs(idef,md,EI_DYNAMICS(ir->eI),cr,constr->lincsd);
+ }
+ if (ir->eConstrAlg == econtSHAKE)
+ {
+ if (cr->dd)
+ {
+ make_shake_sblock_dd(constr,&idef->il[F_CONSTR],&top->cgs,cr->dd);
+ }
+ else
+ {
+ make_shake_sblock_pd(constr,idef,md);
+ }
+ if (ncons > constr->lagr_nalloc)
+ {
+ constr->lagr_nalloc = over_alloc_dd(ncons);
+ srenew(constr->lagr,constr->lagr_nalloc);
+ }
+ }
+ }
+
+ if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL)
+ {
+ settle = &idef->il[F_SETTLE];
+ iO = settle->iatoms[1];
- for (i=0; i<ilist[F_SETTLE].nr; i+=2)
++ iH = settle->iatoms[2];
+ constr->settled =
+ settle_init(md->massT[iO],md->massT[iH],
+ md->invmass[iO],md->invmass[iH],
+ idef->iparams[settle->iatoms[0]].settle.doh,
+ idef->iparams[settle->iatoms[0]].settle.dhh);
+ }
+
+ /* Make a selection of the local atoms for essential dynamics */
+ if (constr->ed && cr->dd)
+ {
+ dd_make_local_ed_indices(cr->dd,constr->ed);
+ }
+}
+
+static void constr_recur(t_blocka *at2con,
+ t_ilist *ilist,t_iparams *iparams,gmx_bool bTopB,
+ int at,int depth,int nc,int *path,
+ real r0,real r1,real *r2max,
+ int *count)
+{
+ int ncon1;
+ t_iatom *ia1,*ia2;
+ int c,con,a1;
+ gmx_bool bUse;
+ t_iatom *ia;
+ real len,rn0,rn1;
+
+ (*count)++;
+
+ ncon1 = ilist[F_CONSTR].nr/3;
+ ia1 = ilist[F_CONSTR].iatoms;
+ ia2 = ilist[F_CONSTRNC].iatoms;
+
+ /* Loop over all constraints connected to this atom */
+ for(c=at2con->index[at]; c<at2con->index[at+1]; c++) {
+ con = at2con->a[c];
+ /* Do not walk over already used constraints */
+ bUse = TRUE;
+ for(a1=0; a1<depth; a1++) {
+ if (con == path[a1])
+ bUse = FALSE;
+ }
+ if (bUse) {
+ ia = constr_iatomptr(ncon1,ia1,ia2,con);
+ /* Flexible constraints currently have length 0, which is incorrect */
+ if (!bTopB)
+ len = iparams[ia[0]].constr.dA;
+ else
+ len = iparams[ia[0]].constr.dB;
+ /* In the worst case the bond directions alternate */
+ if (nc % 2 == 0) {
+ rn0 = r0 + len;
+ rn1 = r1;
+ } else {
+ rn0 = r0;
+ rn1 = r1 + len;
+ }
+ /* Assume angles of 120 degrees between all bonds */
+ if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max) {
+ *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
+ if (debug) {
+ fprintf(debug,"Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0,rn1,sqrt(*r2max));
+ for(a1=0; a1<depth; a1++)
+ fprintf(debug," %d %5.3f",
+ path[a1],
+ iparams[constr_iatomptr(ncon1,ia1,ia2,con)[0]].constr.dA);
+ fprintf(debug," %d %5.3f\n",con,len);
+ }
+ }
+ /* Limit the number of recursions to 1000*nc,
+ * so a call does not take more than a second,
+ * even for highly connected systems.
+ */
+ if (depth + 1 < nc && *count < 1000*nc) {
+ if (ia[1] == at)
+ a1 = ia[2];
+ else
+ a1 = ia[1];
+ /* Recursion */
+ path[depth] = con;
+ constr_recur(at2con,ilist,iparams,
+ bTopB,a1,depth+1,nc,path,rn0,rn1,r2max,count);
+ path[depth] = -1;
+ }
+ }
+ }
+}
+
+static real constr_r_max_moltype(FILE *fplog,
+ gmx_moltype_t *molt,t_iparams *iparams,
+ t_inputrec *ir)
+{
+ int natoms,nflexcon,*path,at,count;
+
+ t_blocka at2con;
+ real r0,r1,r2maxA,r2maxB,rmax,lam0,lam1;
+
+ if (molt->ilist[F_CONSTR].nr == 0 &&
+ molt->ilist[F_CONSTRNC].nr == 0) {
+ return 0;
+ }
+
+ natoms = molt->atoms.nr;
+
+ at2con = make_at2con(0,natoms,molt->ilist,iparams,
+ EI_DYNAMICS(ir->eI),&nflexcon);
+ snew(path,1+ir->nProjOrder);
+ for(at=0; at<1+ir->nProjOrder; at++)
+ path[at] = -1;
+
+ r2maxA = 0;
+ for(at=0; at<natoms; at++) {
+ r0 = 0;
+ r1 = 0;
+
+ count = 0;
+ constr_recur(&at2con,molt->ilist,iparams,
+ FALSE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxA,&count);
+ }
+ if (ir->efep == efepNO) {
+ rmax = sqrt(r2maxA);
+ } else {
+ r2maxB = 0;
+ for(at=0; at<natoms; at++) {
+ r0 = 0;
+ r1 = 0;
+ count = 0;
+ constr_recur(&at2con,molt->ilist,iparams,
+ TRUE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxB,&count);
+ }
+ lam0 = ir->init_lambda;
+ if (EI_DYNAMICS(ir->eI))
+ lam0 += ir->init_step*ir->delta_lambda;
+ rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
+ if (EI_DYNAMICS(ir->eI)) {
+ lam1 = ir->init_lambda + (ir->init_step + ir->nsteps)*ir->delta_lambda;
+ rmax = max(rmax,(1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
+ }
+ }
+
+ done_blocka(&at2con);
+ sfree(path);
+
+ return rmax;
+}
+
+real constr_r_max(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir)
+{
+ int mt;
+ real rmax;
+
+ rmax = 0;
+ for(mt=0; mt<mtop->nmoltype; mt++) {
+ rmax = max(rmax,
+ constr_r_max_moltype(fplog,&mtop->moltype[mt],
+ mtop->ffparams.iparams,ir));
+ }
+
+ if (fplog)
+ fprintf(fplog,"Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n",1+ir->nProjOrder,rmax);
+
+ return rmax;
+}
+
+gmx_constr_t init_constraints(FILE *fplog,
+ gmx_mtop_t *mtop,t_inputrec *ir,
+ gmx_edsam_t ed,t_state *state,
+ t_commrec *cr)
+{
+ int ncon,nset,nmol,settle_type,i,natoms,mt,nflexcon;
+ struct gmx_constr *constr;
+ char *env;
+ t_ilist *ilist;
+ gmx_mtop_ilistloop_t iloop;
+
+ ncon =
+ gmx_mtop_ftype_count(mtop,F_CONSTR) +
+ gmx_mtop_ftype_count(mtop,F_CONSTRNC);
+ nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
+
+ if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
+ {
+ return NULL;
+ }
+
+ snew(constr,1);
+
+ constr->ncon_tot = ncon;
+ constr->nflexcon = 0;
+ if (ncon > 0)
+ {
+ constr->n_at2con_mt = mtop->nmoltype;
+ snew(constr->at2con_mt,constr->n_at2con_mt);
+ for(mt=0; mt<mtop->nmoltype; mt++)
+ {
+ constr->at2con_mt[mt] = make_at2con(0,mtop->moltype[mt].atoms.nr,
+ mtop->moltype[mt].ilist,
+ mtop->ffparams.iparams,
+ EI_DYNAMICS(ir->eI),&nflexcon);
+ for(i=0; i<mtop->nmolblock; i++)
+ {
+ if (mtop->molblock[i].type == mt)
+ {
+ constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
+ }
+ }
+ }
+
+ if (constr->nflexcon > 0)
+ {
+ if (fplog)
+ {
+ fprintf(fplog,"There are %d flexible constraints\n",
+ constr->nflexcon);
+ if (ir->fc_stepsize == 0)
+ {
+ fprintf(fplog,"\n"
+ "WARNING: step size for flexible constraining = 0\n"
+ " All flexible constraints will be rigid.\n"
+ " Will try to keep all flexible constraints at their original length,\n"
+ " but the lengths may exhibit some drift.\n\n");
+ constr->nflexcon = 0;
+ }
+ }
+ if (constr->nflexcon > 0)
+ {
+ please_cite(fplog,"Hess2002");
+ }
+ }
+
+ if (ir->eConstrAlg == econtLINCS)
+ {
+ constr->lincsd = init_lincs(fplog,mtop,
+ constr->nflexcon,constr->at2con_mt,
+ DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
+ ir->nLincsIter,ir->nProjOrder);
+ }
+
+ if (ir->eConstrAlg == econtSHAKE) {
+ if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
+ {
+ gmx_fatal(FARGS,"SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
+ }
+ if (constr->nflexcon)
+ {
+ gmx_fatal(FARGS,"For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
+ }
+ please_cite(fplog,"Ryckaert77a");
+ if (ir->bShakeSOR)
+ {
+ please_cite(fplog,"Barth95a");
+ }
+
+ constr->shaked = shake_init();
+ }
+ }
+
+ if (nset > 0) {
+ please_cite(fplog,"Miyamoto92a");
+
+ /* Check that we have only one settle type */
+ settle_type = -1;
+ iloop = gmx_mtop_ilistloop_init(mtop);
+ while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
+ {
++ for (i=0; i<ilist[F_SETTLE].nr; i+=4)
+ {
+ if (settle_type == -1)
+ {
+ settle_type = ilist[F_SETTLE].iatoms[i];
+ }
+ else if (ilist[F_SETTLE].iatoms[i] != settle_type)
+ {
+ gmx_fatal(FARGS,
+ "The [molecules] section of your topology specifies more than one block of\n"
+ "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
+ "are trying to partition your solvent into different *groups* (e.g. for\n"
+ "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
+ "files specify groups. Otherwise, you may wish to change the least-used\n"
+ "block of molecules with SETTLE constraints into 3 normal constraints.");
+ }
+ }
+ }
+ }
+
+ constr->maxwarn = 999;
+ env = getenv("GMX_MAXCONSTRWARN");
+ if (env)
+ {
+ constr->maxwarn = 0;
+ sscanf(env,"%d",&constr->maxwarn);
+ if (fplog)
+ {
+ fprintf(fplog,
+ "Setting the maximum number of constraint warnings to %d\n",
+ constr->maxwarn);
+ }
+ if (MASTER(cr))
+ {
+ fprintf(stderr,
+ "Setting the maximum number of constraint warnings to %d\n",
+ constr->maxwarn);
+ }
+ }
+ if (constr->maxwarn < 0 && fplog)
+ {
+ fprintf(fplog,"maxwarn < 0, will not stop on constraint errors\n");
+ }
+ constr->warncount_lincs = 0;
+ constr->warncount_settle = 0;
+
+ /* Initialize the essential dynamics sampling.
+ * Put the pointer to the ED struct in constr */
+ constr->ed = ed;
+ if (ed != NULL)
+ {
+ init_edsam(mtop,ir,cr,ed,state->x,state->box);
+ }
+
+ constr->warn_mtop = mtop;
+
+ return constr;
+}
+
+t_blocka *atom2constraints_moltype(gmx_constr_t constr)
+{
+ return constr->at2con_mt;
+}
+
+
+gmx_bool inter_charge_group_constraints(gmx_mtop_t *mtop)
+{
+ const gmx_moltype_t *molt;
+ const t_block *cgs;
+ const t_ilist *il;
+ int mb;
+ int nat,*at2cg,cg,a,ftype,i;
+ gmx_bool bInterCG;
+
+ bInterCG = FALSE;
+ for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++) {
+ molt = &mtop->moltype[mtop->molblock[mb].type];
+
+ if (molt->ilist[F_CONSTR].nr > 0 ||
+ molt->ilist[F_CONSTRNC].nr > 0) {
+ cgs = &molt->cgs;
+ snew(at2cg,molt->atoms.nr);
+ for(cg=0; cg<cgs->nr; cg++) {
+ for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
+ at2cg[a] = cg;
+ }
+
+ for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
+ il = &molt->ilist[ftype];
+ for(i=0; i<il->nr && !bInterCG; i+=3) {
+ if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
+ bInterCG = TRUE;
+ }
+ }
+ sfree(at2cg);
+ }
+ }
+
+ return bInterCG;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff