set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
-set(CMAKE_INSTALL_PREFIX "/usr/local/gromacs" CACHE STRING "" FORCE)
+set(CMAKE_INSTALL_PREFIX "/usr/local/gromacs" CACHE STRING "Installation prefix (installation will need write permissions here)" FORCE)
endif()
if(NOT CMAKE_BUILD_TYPE)
########################################################################
# User input options #
########################################################################
-option(GMX_DOUBLE "Use double precision. Much slower, use only if needed!)" OFF)
+option(GMX_DOUBLE "Use double precision (much slower, use only if you really need it)" OFF)
option(GMX_MPI "Build a parallel (message-passing) version of GROMACS" OFF)
-option(GMX_THREADS "Build a parallel (thread-based) version of GROMACS)" ON)
+option(GMX_THREADS "Build a parallel (thread-based) version of GROMACS (cannot be combined with MPI yet)" ON)
option(GMX_SOFTWARE_INVSQRT "Use GROMACS software 1/sqrt" ON)
mark_as_advanced(GMX_SOFTWARE_INVSQRT)
option(GMX_FAHCORE "Build a library with mdrun functionality" OFF)
mark_as_advanced(GMX_FAHCORE)
-option(GMX_OPENMM "Accelerated execution on GPUs through the OpenMM library" OFF)
+option(GMX_OPENMM "Accelerated execution on GPUs through the OpenMM library (rerun cmake after changing to see relevant options)" OFF)
set(GMX_ACCELERATION "auto"
CACHE STRING "Accelerated kernels. Pick one of: auto, none, SSE, BlueGene, Power6, ia64, altivec")
option(USE_VERSION_H "Generate development version string/information" ON)
mark_as_advanced(USE_VERSION_H)
-option(GMX_DEFAULT_SUFFIX "Use default GROMACS suffixes" ON)
-set(GMX_BINARY_SUFFIX "" CACHE STRING "Suffix for GROMACS binaries (default: _d for double, _mpi for mpi, _mpi_d for mpi and double).")
+option(GMX_DEFAULT_SUFFIX "Use default suffixes for GROMACS binaries and libs (_d for double, _mpi for MPI; rerun cmake after changing to see relevant options)" ON)
+set(GMX_BINARY_SUFFIX "" CACHE STRING "Suffix for GROMACS binaries (default: _d for double, _mpi for MPI, _mpi_d for MPI and double).")
set(GMX_LIBS_SUFFIX ""
- CACHE STRING "Suffix for GROMACS libs (default: _d for double, _mpi for mpi, _mpi_d for mpi and double).")
+ CACHE STRING "Suffix for GROMACS libs (default: _d for double, _mpi for MPI, _mpi_d for MPI and double).")
if (GMX_DEFAULT_SUFFIX)
set(GMX_BINARY_SUFFIX "")
set(GMX_LIBS_SUFFIX "")
# this should at some point become VERSION_LESS
if(NOT Git_FOUND OR Git_VERSION STRLESS "1.5.1")
message("No compatible git version found, won't be able to generate proper development version information.")
+ set(USE_VERSION_H OFF)
endif()
endif()
else()
CFLAGS="$CFLAGS $PTHREAD_CFLAGS -I$srcdir/include"
CXXFLAGS="$CXXFLAGS $PTHREAD_CFLAGS -I$srcdir/include"
CC="$PTHREAD_CC "
- AC_DEFINE(THREAD_PTHREADS,,[Use pthreads for thread_mpi multithreading])
+ AC_DEFINE(THREAD_PTHREADS,,[Use pthreads for multithreading])
# profiling
AC_ARG_ENABLE(tmpi-profiling,
}
+static int gmx_mm_check_and_reset_overflow(void)
+{
+ int MXCSR;
+ int sse_overflow;
+
+ MXCSR = _mm_getcsr();
+ /* The overflow flag is bit 3 in the register */
+ if (MXCSR & 0x0008)
+ {
+ sse_overflow = 1;
+ /* Set the overflow flag to zero */
+ MXCSR = MXCSR & 0xFFF7;
+ _mm_setcsr(MXCSR);
+ }
+ else
+ {
+ sse_overflow = 0;
+ }
+
+ return sse_overflow;
+}
+
static inline __m128d
gmx_mm_invsqrt_pd(__m128d x)
}
+static int gmx_mm_check_and_reset_overflow(void)
+{
+ int MXCSR;
+ int sse_overflow;
+
+ MXCSR = _mm_getcsr();
+ /* The overflow flag is bit 3 in the register */
+ if (MXCSR & 0x0008)
+ {
+ sse_overflow = 1;
+ /* Set the overflow flag to zero */
+ MXCSR = MXCSR & 0xFFF7;
+ _mm_setcsr(MXCSR);
+ }
+ else
+ {
+ sse_overflow = 0;
+ }
+
+ return sse_overflow;
+}
+
+
/************************
* *
* Simple math routines *
These functions set/read the energyhistory_t structure
that is written to checkpoints in checkpoint.c */
-/* initialize the energyhistory_t data structure */
-void init_energyhistory(energyhistory_t * enerhist);
-
/* Set the energyhistory_t data structure from a mdebin structure */
void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin);
void init_mtop(gmx_mtop_t *mtop);
void init_top (t_topology *top);
void init_inputrec(t_inputrec *ir);
+void init_energyhistory(energyhistory_t * enerhist);
+void done_energyhistory(energyhistory_t * enerhist);
void init_gtc_state(t_state *state,int ngtc, int nnhpres, int nhchainlength);
void init_state(t_state *state,int natoms,int ngtc, int nnhpres, int nhchainlength);
aminoacids.n.tdb dna.rtp forcefield.itp spc.itp \
aminoacids.r2b ffbonded.itp gb.itp tip3p.itp \
aminoacids.rtp ffnabonded.itp ions.itp tip4p.itp \
-spce.itp tips3p.itp watermodels.dat
+spce.itp tips3p.itp watermodels.dat tip5p.itp
EXTRA_DIST = ${topol_DATA}
OWT4 15.9994 ; TIP4P
HWT4 1.0080 ; TIP4P
MWT4 0.0000 ; TIP4P
+OWT5 15.9994 ; TIP5P
+HWT5 1.0080 ; TIP5P
+MWT5 0.0000 ; TIP5P
MNH3 0.0000 ; vsite (rigid tetrahedrical NH3 group)
MNH2 0.0000 ; vsite
MCH3 0.0000 ; vsite (rigid CH3 group connected to carbons)
OWT4 8 9.951400 0.0 A 3.15365e-01 6.48520e-01 ; TIP4p O
HWT4 1 4.032000 0.52 A 0.0 0.0 ; TIP4p H
MWT4 0 0.000000 -1.04 A 0.0 0.0 ; TIP 4p vsite
+OWT5 8 9.951400 0.0 A 3.12000e-01 6.69440e-01 ; TIP5p O
+HWT5 1 4.032000 0.241 A 0.0 0.0 ; TIP5p H
+MWT5 0 0.000000 -0.241 A 0.0 0.0 ; TIP5p vsite
OW 8 9.951400 -0.82 A 3.16557e-01 6.50194e-01 ; SPC 0
HW 1 4.032000 0.41 A 0.0 0.0; SPC H
#else
OWT4 8 15.999400 0.0 A 3.15365e-01 6.48520e-01 ; TIP4p O
HWT4 1 1.008000 0.52 A 0.0 0.0 ; TIP4p H
MWT4 0 0.000000 -1.04 A 0.0 0.0 ; TIP4p vsite
+OWT5 8 15.999400 0.0 A 3.12000e-01 6.69440e-01 ; TIP5p O
+HWT5 1 1.008000 0.241 A 0.0 0.0 ; TIP5p H
+MWT5 0 0.000000 -0.241 A 0.0 0.0 ; TIP5p vsite
OW 8 15.999400 -0.82 A 3.16557e-01 6.50194e-01 ; SPC O
HW 1 1.008000 0.41 A 0.0 0.0 ; SPC H
#endif
--- /dev/null
+[ moleculetype ]
+; molname nrexcl
+SOL 2
+
+[ atoms ]
+; id at type res nr res name at name cg nr charge mass
+ 1 OWT5 1 SOL OW 1 0 15.99940
+ 2 HWT5 1 SOL HW1 1 0.241 1.00800
+ 3 HWT5 1 SOL HW2 1 0.241 1.00800
+ 4 MWT5 1 SOL LP1 1 -0.241 0.00000
+ 5 MWT5 1 SOL LP2 1 -0.241 0.00000
+
+#ifndef FLEXIBLE
+
+[ settles ]
+; i funct doh dhh
+1 1 0.09572 0.15139
+
+#else
+
+[ bonds ]
+; i j funct length force.c.
+1 2 1 0.09572 502416.0 0.09572 502416.0
+1 3 1 0.09572 502416.0 0.09572 502416.0
+
+[ angles ]
+; i j k funct angle force.c.
+2 1 3 1 104.52 628.02 104.52 628.02
+
+#endif
+
+
+[ virtual_sites3 ]
+; Vsite from funct a b c
+4 1 2 3 4 -0.344908262 -0.34490826 -6.4437903493
+5 1 2 3 4 -0.344908262 -0.34490826 6.4437903493
+
+
+[ exclusions ]
+1 2 3 4 5
+2 1 3 4 5
+3 1 2 4 5
+4 1 2 3 5
+5 1 2 3 4
+
+
+; The positions of the vsites are computed as follows:
+;
+; LP1 LP2
+;
+; O
+;
+; H1 H2
+;
+; angle A (H1-O-H2) = 104.52
+; angle B (M1-O-M2) = 109.47
+; dist C (H-O) = 0.09572 nm
+; dist D (M-O) = 0.070 nm
+;
+;atom x y z
+;O 0.0 0.0 0.0
+;H1 0.585882276 0.756950327 0.0
+;H2 0.585882276 -0.756950327 0.0
+;M1 -0.404151276 0.0 0.571543301
+;M2 -0.404151276 0.0 -0.571543301
+; Dummy pos x4 = x1 + a4*(x2-x1) + b4*(x3-x1) + c4*((x2-x1) x (x3-x1))
+; Dummy pos x5 = x1 + a5*(x2-x1) + b5*(x3-x1) + c5*((x2-x1) x (x3-x1))
+; a4 = b4 = a5 = b5 = (D*cos(B/2)) / (2*C*cos(A/2)) = -0.34490826
+; c5 = -c4 = (D * sin(B/2))/ (C^2 * sin(A)) = 6.4437903
}
/* remove a t_fileio into the list. We assume the fio is locked, and we leave
- it locked. */
-static void gmx_fio_remove(t_fileio *fio, gmx_bool global_lock)
+ it locked.
+ NOTE: We also assume that the open_file_mutex has been locked */
+static void gmx_fio_remove(t_fileio *fio)
{
t_fileio *prev;
-#ifdef GMX_THREADS
- /* first lock the big open_files mutex. */
- /* We don't want two processes operating on this list at the same time */
- if (global_lock)
- tMPI_Thread_mutex_lock(&open_file_mutex);
-#endif
-
/* lock prev, because we're changing it */
gmx_fio_lock(fio->prev);
/* and make sure we point nowhere in particular */
fio->next=fio->prev=fio;
-
-#ifdef GMX_THREADS
- /* now unlock the big open_files mutex. */
- if (global_lock)
- tMPI_Thread_mutex_unlock(&open_file_mutex);
-#endif
-
}
{
int rc = 0;
+#ifdef GMX_THREADS
+ /* first lock the big open_files mutex. */
+ /* We don't want two processes operating on the list at the same time */
+ tMPI_Thread_mutex_lock(&open_file_mutex);
+#endif
+
gmx_fio_lock(fio);
/* first remove it from the list */
- gmx_fio_remove(fio, TRUE);
+ gmx_fio_remove(fio);
rc=gmx_fio_close_locked(fio);
gmx_fio_unlock(fio);
sfree(fio);
+#ifdef GMX_THREADS
+ tMPI_Thread_mutex_unlock(&open_file_mutex);
+#endif
+
return rc;
}
if (cur->fp == fp)
{
rc=gmx_fio_close_locked(cur);
- gmx_fio_remove(cur,FALSE);
+ gmx_fio_remove(cur);
gmx_fio_stop_getting_next(cur);
break;
}
grid_index(gmx_ana_nbsearch_t *d, const ivec cell)
{
return cell[XX] + cell[YY] * d->ncelldim[XX]
- + cell[ZZ] * d->ncelldim[YY] * d->ncelldim[ZZ];
+ + cell[ZZ] * d->ncelldim[XX] * d->ncelldim[YY];
}
/*! \brief
eks->mvcos = 0;
}
-static void init_energyhistory(energyhistory_t *enh)
+void init_energyhistory(energyhistory_t * enerhist)
{
- enh->ener_ave = NULL;
- enh->ener_sum = NULL;
- enh->ener_sum_sim = NULL;
- enh->nener = 0;
+ enerhist->nener = 0;
+
+ enerhist->ener_ave = NULL;
+ enerhist->ener_sum = NULL;
+ enerhist->ener_sum_sim = NULL;
+ enerhist->dht = NULL;
+
+ enerhist->nsteps = 0;
+ enerhist->nsum = 0;
+ enerhist->nsteps_sim = 0;
+ enerhist->nsum_sim = 0;
+
+ enerhist->dht = NULL;
+}
+
+static void done_delta_h_history(delta_h_history_t *dht)
+{
+ int i;
+
+ for(i=0; i<dht->nndh; i++)
+ {
+ sfree(dht->dh[i]);
+ }
+ sfree(dht->dh);
+ sfree(dht->ndh);
+}
+
+void done_energyhistory(energyhistory_t * enerhist)
+{
+ sfree(enerhist->ener_ave);
+ sfree(enerhist->ener_sum);
+ sfree(enerhist->ener_sum_sim);
+ sfree(enerhist->dht);
+
+ if (enerhist->dht != NULL)
+ {
+ done_delta_h_history(enerhist->dht);
+ sfree(enerhist->dht);
+ }
}
void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
}
}
-static char *fgets3(FILE *fp,char ptr[],int *len)
+/* reads a line into ptr, adjusting len and renewing ptr if neccesary */
+static char *fgets3(FILE *fp,char **ptr,int *len, int maxlen)
{
char *p;
- int slen;
+ int len_remaining=*len; /* remaining amount of allocated bytes in buf */
+ int curp=0; /* current position in buf to read into */
- if (fgets(ptr,*len-1,fp) == NULL)
+ do
{
- return NULL;
- }
- p = ptr;
- while ((strchr(ptr,'\n') == NULL) && (!feof(fp)))
- {
- /* This line is longer than len characters, let's increase len! */
- *len += STRLEN;
- p += STRLEN;
- srenew(ptr,*len);
- if (fgets(p-1,STRLEN,fp) == NULL)
+ if (len_remaining < 2)
{
- break;
+ if ( *len + STRLEN < maxlen)
+ {
+ /* This line is longer than len characters, let's increase len! */
+ *len += STRLEN;
+ len_remaining += STRLEN;
+ srenew(*ptr, *len);
+ }
+ else
+ {
+ /*something is wrong, we'll just keep reading and return NULL*/
+ len_remaining = STRLEN;
+ curp=0;
+ }
}
+ if (fgets(*ptr + curp, len_remaining, fp)==NULL)
+ {
+ /* if last line, skip */
+ return NULL;
+ }
+ curp += len_remaining-1; /* overwrite the nul char in next iteration */
+ len_remaining = 1;
+ }
+ while ((strchr(*ptr,'\n') == NULL) && (!feof(fp)));
+
+ if ( *len + STRLEN >= maxlen )
+ {
+ return NULL; /* this line was too long */
}
+
if (feof(fp))
{
/* We reached EOF before '\n', skip this last line. */
return NULL;
}
- slen = strlen(ptr);
- if (ptr[slen-1] == '\n')
{
- ptr[slen-1] = '\0';
+ /* now remove newline */
+ int slen = strlen(*ptr);
+ if ((*ptr)[slen-1] == '\n')
+ {
+ (*ptr)[slen-1] = '\0';
+ }
}
- return ptr;
+ return *ptr;
}
static int wordcount(char *ptr)
*legend = NULL;
}
- while ((ptr = fgets3(fp,tmpbuf,&len)) != NULL && ptr[0]!='&') {
+ while ((ptr = fgets3(fp,&tmpbuf,&len,10*STRLEN)) != NULL && ptr[0]!='&') {
line++;
trim(ptr);
if (ptr[0] == '@') {
if (MASTER(cr))
{
- /* Update mdebin with energy history if appending to output files */
- if ( Flags & MD_APPENDFILES )
+ if (opt2bSet("-cpi",nfile,fnm))
{
- restore_energyhistory_from_state(mdebin,&state_global->enerhist);
+ /* Update mdebin with energy history if appending to output files */
+ if ( Flags & MD_APPENDFILES )
+ {
+ restore_energyhistory_from_state(mdebin,&state_global->enerhist);
+ }
+ else
+ {
+ /* We might have read an energy history from checkpoint,
+ * free the allocated memory and reset the counts.
+ */
+ done_energyhistory(&state_global->enerhist);
+ init_energyhistory(&state_global->enerhist);
+ }
}
- /* Set the initial energy history in state to zero by updating once */
+ /* Set the initial energy history in state by updating once */
update_energyhistory(&state_global->enerhist,mdebin);
}
static int factorize(int n,int **fac,int **mfac)
{
- int d,ndiv;
-
- /* Decompose n in factors */
- snew(*fac,n/2);
- snew(*mfac,n/2);
- d = 2;
- ndiv = 0;
- while (n > 1)
- {
- while (n % d == 0)
- {
- if (ndiv == 0 || (*fac)[ndiv-1] != d)
- {
- ndiv++;
- (*fac)[ndiv-1] = d;
- }
- (*mfac)[ndiv-1]++;
- n /= d;
- }
- d++;
- }
-
- return ndiv;
-}
+ int d,ndiv;
-static gmx_bool is_prime(int n)
-{
- int i;
-
- i = 2;
- while (i*i <= n)
+ /* Decompose n in factors */
+ snew(*fac,n/2);
+ snew(*mfac,n/2);
+ d = 2;
+ ndiv = 0;
+ while (n > 1)
{
- if (n % i == 0)
+ while (n % d == 0)
{
- return FALSE;
+ if (ndiv == 0 || (*fac)[ndiv-1] != d)
+ {
+ ndiv++;
+ (*fac)[ndiv-1] = d;
+ }
+ (*mfac)[ndiv-1]++;
+ n /= d;
}
- i++;
+ d++;
}
+
+ return ndiv;
+}
- return TRUE;
+static gmx_bool largest_divisor(int n)
+{
+ int ndiv,*div,*mdiv,ldiv;
+
+ ndiv = factorize(n,&div,&mdiv);
+ ldiv = div[ndiv-1];
+ sfree(div);
+ sfree(mdiv);
+
+ return ldiv;
}
static int lcd(int n1,int n2)
gmx_bool bInterCGBondeds,gmx_bool bInterCGMultiBody)
{
int npme,nkx,nky;
+ int ldiv;
real limit;
if (MASTER(cr))
{
- if (cr->nnodes > 12 && is_prime(cr->nnodes))
+ if (cr->nnodes > 12)
{
- gmx_fatal(FARGS,"The number of nodes you selected (%d) is a large prime. In most cases this will lead to bad performance. Choose a non-prime number, or set the decomposition (option -dd) manually.",cr->nnodes);
+ ldiv = largest_divisor(cr->nnodes);
+ /* Check if the largest divisor is more than nnodes^2/3 */
+ if (ldiv*ldiv*ldiv > cr->nnodes*cr->nnodes)
+ {
+ gmx_fatal(FARGS,"The number of nodes you selected (%d) contains a large prime factor %d. In most cases this will lead to bad performance. Choose a number with smaller prime factors or set the decomposition (option -dd) manually.",
+ cr->nnodes,ldiv);
+ }
}
if (EEL_PME(ir->coulombtype))
}
-void init_energyhistory(energyhistory_t * enerhist)
-{
- enerhist->nener = 0;
-
- enerhist->ener_ave = NULL;
- enerhist->ener_sum = NULL;
- enerhist->ener_sum_sim = NULL;
- enerhist->dht = NULL;
-
- enerhist->nsteps = 0;
- enerhist->nsum = 0;
- enerhist->nsteps_sim = 0;
- enerhist->nsum_sim = 0;
-}
-
void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
{
int i;
pme->nky <= pme->pme_order*(pme->nnodes_minor > 1 ? 2 : 1) ||
pme->nkz <= pme->pme_order)
{
- gmx_fatal(FARGS,"The pme grid dimensions need to be larger than pme_order (%d) and in parallel larger than 2*pme_ordern for x and/or y",pme->pme_order);
+ gmx_fatal(FARGS,"The pme grid dimensions need to be larger than pme_order (%d) and in parallel larger than 2*pme_order for x and/or y",pme->pme_order);
}
if (pme->nnodes > 1) {
snew(td->v,td->nx);
snew(td->f,td->nx);
}
- for(i=td->nx0; (i<td->nx); i++)
+ for(i=0; (i<td->nx); i++)
td->x[i] = i/tabscale;
}
td->f[i] = Ftab;
}
+ /* Continue the table linearly from nx0 to 0.
+ * These values are only required for energy minimization with overlap or TPI.
+ */
+ for(i=td->nx0-1; i>=0; i--) {
+ td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
+ td->f[i] = td->f[i+1];
+ }
+
#ifdef DEBUG_SWITCH
gmx_fio_fclose(fp);
#endif
#include "pme.h"
#include "gbutil.h"
+#if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_X86_64_SSE2) )
+#if defined(GMX_DOUBLE)
+#include "gmx_sse2_double.h"
+#else
+#include "gmx_sse2_single.h"
+#endif
+#endif
+
static void global_max(t_commrec *cr,int *n)
{
int nbin;
double invbinw,*bin,refvolshift,logV,bUlogV;
real dvdl,prescorr,enercorr,dvdlcorr;
+ gmx_bool bEnergyOutOfBounds;
const char *tpid_leg[2]={"direct","reweighted"};
/* Since numerical problems can lead to extreme negative energies
enerd->term[F_DISPCORR] = enercorr;
enerd->term[F_EPOT] += enercorr;
enerd->term[F_PRES] += prescorr;
- enerd->term[F_DVDL] += dvdlcorr;
-
+ enerd->term[F_DVDL] += dvdlcorr;
+
+ epot = enerd->term[F_EPOT];
+ bEnergyOutOfBounds = FALSE;
+#if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_X86_64_SSE2) )
+ /* With SSE the energy can overflow, check for this */
+ if (gmx_mm_check_and_reset_overflow())
+ {
+ if (debug)
+ {
+ fprintf(debug,"Found an SSE overflow, assuming the energy is out of bounds\n");
+ }
+ bEnergyOutOfBounds = TRUE;
+ }
+#endif
/* If the compiler doesn't optimize this check away
* we catch the NAN energies. With tables extreme negative
* energies might occur close to r=0.
*/
- epot = enerd->term[F_EPOT];
if (epot != epot || epot*beta < bU_neg_limit)
+ {
+ bEnergyOutOfBounds = TRUE;
+ }
+ if (bEnergyOutOfBounds)
{
if (debug)
{
if (MASTER(cr))
{
- /* Update mdebin with energy history if appending to output files */
- if ( Flags & MD_APPENDFILES )
+ if (opt2bSet("-cpi",nfile,fnm))
{
- restore_energyhistory_from_state(mdebin,&state_global->enerhist);
+ /* Update mdebin with energy history if appending to output files */
+ if ( Flags & MD_APPENDFILES )
+ {
+ restore_energyhistory_from_state(mdebin,&state_global->enerhist);
+ }
+ else
+ {
+ /* We might have read an energy history from checkpoint,
+ * free the allocated memory and reset the counts.
+ */
+ done_energyhistory(&state_global->enerhist);
+ init_energyhistory(&state_global->enerhist);
+ }
}
- /* Set the initial energy history in state to zero by updating once */
+ /* Set the initial energy history in state by updating once */
update_energyhistory(&state_global->enerhist,mdebin);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff