# Binary and library suffix options
include(gmxManageSuffixes)
-##################################################################
+################################################################
# Shared library settings
-##################################################################
-if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin")
+################################################################
+if((NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin") OR ((CMAKE_SYSTEM_VERSION VERSION_GREATER 8.0) AND (CMAKE_VERSION VERSION_GREATER 2.8.11)))
if(GMX_LIB_INSTALL_DIR STREQUAL "lib")
set(CMAKE_BUILD_WITH_INSTALL_RPATH TRUE)
endif()
- set(CMAKE_INSTALL_RPATH "\$ORIGIN/../${GMX_LIB_INSTALL_DIR}")
+ if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin")
+ set(CMAKE_INSTALL_RPATH "\$ORIGIN/../${GMX_LIB_INSTALL_DIR}")
+ else()
+ set(CMAKE_INSTALL_RPATH "@executable_path/../${GMX_LIB_INSTALL_DIR}")
+ endif()
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
+ set(CMAKE_MACOSX_RPATH 1)
else()
+ # We are on Darwin/OSX, and cmake cannot handle proper RPATHs
if(CMAKE_SYSTEM_VERSION VERSION_GREATER 8.0) #rpath supported for >10.4
set(CMAKE_INSTALL_NAME_DIR "@rpath")
set(GMX_EXE_LINKER_FLAGS ${GMX_EXE_LINKER_FLAGS} "-Wl,-rpath,@executable_path/../${GMX_LIB_INSTALL_DIR}")
else()
find_package(CUDA 4.0 ${FIND_CUDA_QUIETLY})
endif()
+ # Cmake 2.8.12 (and CMake 3.0) introduced a new bug where the cuda
+ # library dir is added twice as an rpath on APPLE, which in turn causes
+ # the install_name_tool to wreck the binaries when it tries to remove this
+ # path. Since this is set inside the cuda module, we remove the extra rpath
+ # added in the library string - an rpath is not a library anyway, and at
+ # least for Gromacs this works on all CMake versions. This should be
+ # reasonably future-proof, since newer versions of CMake appear to handle
+ # the rpath automatically based on the provided library path, meaning
+ # the explicit rpath specification is no longer needed.
+ if(APPLE AND (CMAKE_VERSION VERSION_GREATER 2.8.11))
+ foreach(elem ${CUDA_LIBRARIES})
+ if(elem MATCHES "-Wl,.*")
+ list(REMOVE_ITEM CUDA_LIBRARIES ${elem})
+ endif()
+ endforeach(elem)
+ endif()
endif()
# Depending on the current vale of GMX_GPU and GMX_GPU_AUTO:
mark_as_advanced(CUDA_HOST_COMPILER CUDA_HOST_COMPILER_OPTIONS)
endif()
+ if(APPLE AND CMAKE_C_COMPILER_ID MATCHES "GNU")
+ # Some versions of gcc-4.8 and gcc-4.9 produce errors (in particular on OS X)
+ # if we do not use -D__STRICT_ANSI__. It is harmless, so we might as well add it for all versions.
+ set(CUDA_HOST_COMPILER_OPTIONS "${CUDA_HOST_COMPILER_OPTIONS}-D__STRICT_ANSI__;")
+ endif()
+
# on Linux we need to add -fPIC when building shared gmx libs
# Note: will add -fPIC for any compiler that supports it as it shouldn't hurt
if(BUILD_SHARED_LIBS)
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald, recommended
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tip4pew TIP4P-Ew TIP 4-point optimized with Ewald
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
; rtp residue to rtp building block table
;GMX Force-field
HISD HSD
+HIS1 HSD
HISE HSE
HISH HSP
LYSN LSN
MNH2 0.0000 ; vsite
MCH3 0.0000 ; vsite (rigid CH3 group connected to carbons)
MCH3S 0.0000 ; vsite (rigid CH3 group connected to S)
+MW 0.0 ; Dummy mass in rigid tyrosine rings
; DNA and RNA section
P 30.974 ; Phosphorus ### For DNA
P2 30.974000 ; pyrophosphate phosphorus (see toppar_all27_na_nad_ppi.str) ### For DNA
MNH2 0 0.000000 0.00 A 0.0 0.0
MCH3 0 0.000000 0.00 A 0.0 0.0
MCH3S 0 0.000000 0.00 A 0.0 0.0
+MW 0 0.000000 0.00 A 0.0 0.0
; Ions and noble gases (useful for tutorials)
Cu2+ 29 63.54600 2.00 A 2.08470e-01 4.76976e+00
Ar 18 39.94800 0.00 A 3.41000e-01 2.74580e-02
MNH2 0 0 0 0 0 ; dummy mass
MCH3 0 0 0 0 0 ; dummy mass
MCH3S 0 0 0 0 0 ; dummy mass
+ MW 0 0 0 0 0
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
tip3p TIP3P TIP 3-point, recommended
tip4p TIP4P TIP 4-point
tips3p TIPS3P CHARMM TIP 3-point with LJ on H's (note: twice as slow in GROMACS)
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
spc SPC simple point charge
spce SPC/E extended simple point charge
HB1 opls_140 0.060 2
HB2 opls_140 0.060 2
CG opls_267 0.520 3
- OD1 opls_269 -0.530 3
- OD2 opls_268 -0.440 4
+ OD1 opls_269 -0.440 3
+ OD2 opls_268 -0.530 4
HD2 opls_270 0.450 4
C opls_235 0.500 5
O opls_236 -0.500 5
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: http://redmine.gromacs.org/issues/1348
+
[ moleculetype ]
; molname nrexcl
SOL 2
4 opls_120 1 SOL LP1 1 -0.241
5 opls_120 1 SOL LP2 1 -0.241
+#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 ]
; The position of the virtual site is computed as follows:
;
tip4p TIP4P TIP 4-point, recommended
tip3p TIP3P TIP 3-point
-tip5p TIP5P TIP 5-point
+tip5p TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
tip5pe TIP5P TIP 5-point improved for Ewald sums
spc SPC simple point charge
spce SPC/E extended simple point charge
if (bSeekForwardOnly)
{
- low = gmx_ftell(fp);
+ low = gmx_ftell(fp)-header_size;
}
if (gmx_fseek(fp, 0, SEEK_END))
{
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * Copyright (c) 2013, by the GROMACS development team, led by
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
int natoms; /* number of atoms (atoms, x, v, f) */
real t0; /* time of the first frame, needed *
* for skipping frames with -dt */
+ real tf; /* internal frame time - DO NOT CHANGE */
real tpf; /* time of the previous frame, not */
/* the read, but real file frames */
real tppf; /* time of two frames ago */
fr->bDouble = FALSE;
fr->natoms = -1;
fr->t0 = 0;
+ fr->tf = 0;
fr->tpf = 0;
fr->tppf = 0;
fr->title = NULL;
int ftp;
bRet = FALSE;
- pt = fr->time;
+ pt = fr->tf;
do
{
clear_trxframe(fr, FALSE);
fr->tppf = fr->tpf;
- fr->tpf = fr->time;
+ fr->tpf = fr->tf;
if (status->tng)
{
/* DvdS 2005-05-31: this has been fixed along with the increased
* accuracy of the control over -b and -e options.
*/
- if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN)))
+ if (bTimeSet(TBEGIN) && (fr->tf < rTimeValue(TBEGIN)))
{
if (xtc_seek_time(status->fio, rTimeValue(TBEGIN), fr->natoms, TRUE))
{
gmx_fio_getname(status->fio));
#endif
}
+ fr->tf = fr->time;
if (bRet)
{
#endif
break;
}
+ fr->tf = fr->time;
/* Return FALSE if we read a frame that's past the set ending time. */
if (!bFirst && (!(fr->flags & TRX_DONT_SKIP) && check_times(fr->time) > 0))
int i, j, Dih, Chi;
j = 0;
- for (Dih = 0; (Dih < NONCHI+maxchi); Dih++)
+ /* NONCHI points to chi1, therefore we have to start counting there. */
+ for (Dih = NONCHI; (Dih < NONCHI+maxchi); Dih++)
{
for (i = 0; (i < nlist); i++)
{
sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
fprintf(stderr, " and %s ", hisfile);
fp = xvgropen(hisfile, histitle, "number", "", oenv);
- fprintf(fp, "@ xaxis tick on\n");
- fprintf(fp, "@ xaxis tick major 1\n");
- fprintf(fp, "@ type xy\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ xaxis tick on\n");
+ fprintf(fp, "@ xaxis tick major 1\n");
+ fprintf(fp, "@ type xy\n");
+ }
for (k = 0; (k < nbin); k++)
{
if (bNormalize)
fprintf(fp, "%5d %10d\n", k, chi_prhist[k]);
}
}
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
gmx_ffclose(fp);
}
{
if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
{
- if (output_env_get_print_xvgr_codes(oenv))
- {
- fprintf(out, "&\n");
- }
+ fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
fprintf(out, "%10.4f %10.5f\n",
x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
}
- if (output_env_get_print_xvgr_codes(oenv))
- {
- fprintf(out, "&\n");
- }
+ fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
}
gmx_ffclose(out);
out = xvgropen(outfile, "Subspace overlap",
"Eigenvectors of trajectory 2", "Overlap", oenv);
- fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
- noutvec);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
+ }
overlap = 0;
for (x = 0; x < nvec2; x++)
{
{
if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
{
- fprintf(xvgrout, "&\n");
+ fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
}
}
if (s < nset-1)
{
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
}
gmx_ffclose(fp);
}
fprintf(stdout, "Set %3d: err.est. %g a %g tau1 %g tau2 %g\n",
s+1, ee, a, tau1, tau2);
- fprintf(fp, "@ legend string %d \"av %f\"\n", 2*s, av[s]);
- fprintf(fp, "@ legend string %d \"ee %6g\"\n",
- 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt));
+ if (output_env_get_xvg_format(oenv) == exvgXMGR)
+ {
+ fprintf(fp, "@ legend string %d \"av %f\"\n", 2*s, av[s]);
+ fprintf(fp, "@ legend string %d \"ee %6g\"\n", 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt));
+ }
+ else if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
+ {
+ fprintf(fp, "@ s%d legend \"av %f\"\n", 2*s, av[s]);
+ fprintf(fp, "@ s%d legend \"ee %6g\"\n", 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt));
+ }
for (i = 0; i < nbs; i++)
{
fprintf(fp, "%g %g %g\n", tbs[i], sig[s]*sqrt(ybs[i]/(n*dt)),
s+1, sig[s]*anal_ee_inf(ac_fit, n*dt),
ac_fit[1], ac_fit[0], ac_fit[2]);
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
for (i = 0; i < nbs; i++)
{
fprintf(fp, "%g %g\n", tbs[i],
}
if (s < nset-1)
{
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
}
sfree(fitsig);
}
if (s < nset-1)
{
- fprintf(out, "&\n");
+ fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
}
gmx_ffclose(out);
{
maxstat = max(maxstat, angstat[i]*norm_fac);
}
- fprintf(out, "@with g0\n");
- fprintf(out, "@ world xmin -180\n");
- fprintf(out, "@ world xmax 180\n");
- fprintf(out, "@ world ymin 0\n");
- fprintf(out, "@ world ymax %g\n", maxstat*1.05);
- fprintf(out, "@ xaxis tick major 60\n");
- fprintf(out, "@ xaxis tick minor 30\n");
- fprintf(out, "@ yaxis tick major 0.005\n");
- fprintf(out, "@ yaxis tick minor 0.0025\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "@with g0\n");
+ fprintf(out, "@ world xmin -180\n");
+ fprintf(out, "@ world xmax 180\n");
+ fprintf(out, "@ world ymin 0\n");
+ fprintf(out, "@ world ymax %g\n", maxstat*1.05);
+ fprintf(out, "@ xaxis tick major 60\n");
+ fprintf(out, "@ xaxis tick minor 30\n");
+ fprintf(out, "@ yaxis tick major 0.005\n");
+ fprintf(out, "@ yaxis tick minor 0.0025\n");
+ }
}
for (i = first; (i <= last); i++)
{
strcpy(hhisfile, hisfile);
strcat(hhisfile, ".xvg");
fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
- fprintf(fp, "@ with g0\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ with g0\n");
+ }
xvgr_world(fp, -180, 0, 180, 0.1, oenv);
- fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
- fprintf(fp, "@ xaxis tick on\n");
- fprintf(fp, "@ xaxis tick major 90\n");
- fprintf(fp, "@ xaxis tick minor 30\n");
- fprintf(fp, "@ xaxis ticklabel prec 0\n");
- fprintf(fp, "@ yaxis tick off\n");
- fprintf(fp, "@ yaxis ticklabel off\n");
- fprintf(fp, "@ type xy\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
+ fprintf(fp, "@ xaxis tick on\n");
+ fprintf(fp, "@ xaxis tick major 90\n");
+ fprintf(fp, "@ xaxis tick minor 30\n");
+ fprintf(fp, "@ xaxis ticklabel prec 0\n");
+ fprintf(fp, "@ yaxis tick off\n");
+ fprintf(fp, "@ yaxis ticklabel off\n");
+ fprintf(fp, "@ type xy\n");
+ }
if (bSSHisto)
{
for (k = 0; (k < 3); k++)
}
}
}
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
gmx_ffclose(fp);
if (bSSHisto)
{
for (k = 0; (k < 3); k++)
{
- fprintf(ssfp[k], "&\n");
+ fprintf(ssfp[k], "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
gmx_ffclose(ssfp[k]);
}
}
FILE *fp;
fp = xvgropen(fn, title, xaxis, yaxis, oenv);
- fprintf(fp, "@ with g0\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ with g0\n");
+ }
xvgr_world(fp, -180, -180, 180, 180, oenv);
- fprintf(fp, "@ xaxis tick on\n");
- fprintf(fp, "@ xaxis tick major 90\n");
- fprintf(fp, "@ xaxis tick minor 30\n");
- fprintf(fp, "@ xaxis ticklabel prec 0\n");
- fprintf(fp, "@ yaxis tick on\n");
- fprintf(fp, "@ yaxis tick major 90\n");
- fprintf(fp, "@ yaxis tick minor 30\n");
- fprintf(fp, "@ yaxis ticklabel prec 0\n");
- fprintf(fp, "@ s0 type xy\n");
- fprintf(fp, "@ s0 symbol 2\n");
- fprintf(fp, "@ s0 symbol size 0.410000\n");
- fprintf(fp, "@ s0 symbol fill 1\n");
- fprintf(fp, "@ s0 symbol color 1\n");
- fprintf(fp, "@ s0 symbol linewidth 1\n");
- fprintf(fp, "@ s0 symbol linestyle 1\n");
- fprintf(fp, "@ s0 symbol center false\n");
- fprintf(fp, "@ s0 symbol char 0\n");
- fprintf(fp, "@ s0 skip 0\n");
- fprintf(fp, "@ s0 linestyle 0\n");
- fprintf(fp, "@ s0 linewidth 1\n");
- fprintf(fp, "@ type xy\n");
-
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ xaxis tick on\n");
+ fprintf(fp, "@ xaxis tick major 90\n");
+ fprintf(fp, "@ xaxis tick minor 30\n");
+ fprintf(fp, "@ xaxis ticklabel prec 0\n");
+ fprintf(fp, "@ yaxis tick on\n");
+ fprintf(fp, "@ yaxis tick major 90\n");
+ fprintf(fp, "@ yaxis tick minor 30\n");
+ fprintf(fp, "@ yaxis ticklabel prec 0\n");
+ fprintf(fp, "@ s0 type xy\n");
+ fprintf(fp, "@ s0 symbol 2\n");
+ fprintf(fp, "@ s0 symbol size 0.410000\n");
+ fprintf(fp, "@ s0 symbol fill 1\n");
+ fprintf(fp, "@ s0 symbol color 1\n");
+ fprintf(fp, "@ s0 symbol linewidth 1\n");
+ fprintf(fp, "@ s0 symbol linestyle 1\n");
+ fprintf(fp, "@ s0 symbol center false\n");
+ fprintf(fp, "@ s0 symbol char 0\n");
+ fprintf(fp, "@ s0 skip 0\n");
+ fprintf(fp, "@ s0 linestyle 0\n");
+ fprintf(fp, "@ s0 linewidth 1\n");
+ fprintf(fp, "@ type xy\n");
+ }
return fp;
}
if (clustidfn)
{
fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
- fprintf(fp, "@ s0 symbol 2\n");
- fprintf(fp, "@ s0 symbol size 0.2\n");
- fprintf(fp, "@ s0 linestyle 0\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ s0 symbol 2\n");
+ fprintf(fp, "@ s0 symbol size 0.2\n");
+ fprintf(fp, "@ s0 linestyle 0\n");
+ }
for (i = 0; i < nf; i++)
{
fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
if (sizefn)
{
fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
- fprintf(fp, "@g%d type %s\n", 0, "bar");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@g%d type %s\n", 0, "bar");
+ }
}
snew(structure, nf);
fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
DD = pow(10.0, 0.125);
fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);
xvgr_legend(fp, asize(leg), leg, oenv);
- fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
- fprintf(fp, "@ xaxes scale Logarithmic\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
+ fprintf(fp, "@ xaxes scale Logarithmic\n");
+ }
for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
{
f = calc_fluidicity(Delta, toler);
oenv);
xvgr_legend(out, 0, NULL, oenv);
j = 0;
- for (m = 0; (m < egNR+egSP); m++)
+ if (output_env_get_print_xvgr_codes(oenv))
{
- if (egrp_use[m])
+ char str1[STRLEN], str2[STRLEN];
+ if (output_env_get_xvg_format(oenv) == exvgXMGR)
{
- fprintf(out, "@ legend string %d \"%s\"\n", j++, egrp_nm[m]);
+ sprintf(str1, "@ legend string ");
+ sprintf(str2, " ");
}
- }
- if (bFree)
- {
- fprintf(out, "@ legend string %d \"%s\"\n", j++, "Free");
- }
- if (bFree)
- {
- fprintf(out, "@ legend string %d \"%s\"\n", j++, "Diff");
- }
- fprintf(out, "@TYPE xy\n");
- fprintf(out, "#%3s", "grp");
- for (m = 0; (m < egNR+egSP); m++)
- {
- if (egrp_use[m])
+ else
{
- fprintf(out, " %9s", egrp_nm[m]);
+ sprintf(str1, "@ s");
+ sprintf(str2, " legend ");
}
+
+ for (m = 0; (m < egNR+egSP); m++)
+ {
+ if (egrp_use[m])
+ {
+ fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
+ }
+ }
+ if (bFree)
+ {
+ fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
+ }
+ if (bFree)
+ {
+ fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
+ }
+ fprintf(out, "@TYPE xy\n");
+ fprintf(out, "#%3s", "grp");
+
+ for (m = 0; (m < egNR+egSP); m++)
+ {
+ if (egrp_use[m])
+ {
+ fprintf(out, " %9s", egrp_nm[m]);
+ }
+ }
+ if (bFree)
+ {
+ fprintf(out, " %9s", "Free");
+ }
+ if (bFree)
+ {
+ fprintf(out, " %9s", "Diff");
+ }
+ fprintf(out, "\n");
}
- if (bFree)
- {
- fprintf(out, " %9s", "Free");
- }
- if (bFree)
- {
- fprintf(out, " %9s", "Diff");
- }
- fprintf(out, "\n");
for (i = 0; (i < ngroups); i++)
{
fprintf(out, "%3.0f", groupnr[i]);
{
fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
"Time (ps)", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(fort, "%s", orinst_sub);
}
fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
"Orientation restraint deviation",
"Time (ps)", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(fodt, "%s", orinst_sub);
}
out = xvgropen(opt2fn("-ora", NFILE, fnm),
"Average calculated orientations",
"Restraint label", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "%s", orinst_sub);
}
out = xvgropen(opt2fn("-oda", NFILE, fnm),
"Average restraint deviation",
"Restraint label", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "%s", orinst_sub);
}
out = xvgropen(opt2fn("-odr", NFILE, fnm),
"RMS orientation restraint deviations",
"Restraint label", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "%s", orinst_sub);
}
hi = d+dd;
fprintf(out, "%5d %5d %1d %5d %10d %10g %10g %10g %10g\n",
ind_grp[i]+1, ind_grp[j]+1, 1, k, 1,
- lo, hi, hi+1, 1.0);
+ lo, hi, hi+disre_up2, 1.0);
}
}
}
if (bCalcN)
{
+ char **legend;
+
+ snew(legend, 5);
+ for (i = 0; i < 5; i++)
+ {
+ snew(legend[i], STRLEN);
+ }
tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);
fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
"Increase in number of contacts", "Residue", "Ratio", oenv);
- fprintf(fp, "@ legend on\n");
- fprintf(fp, "@ legend box on\n");
- fprintf(fp, "@ legend loctype view\n");
- fprintf(fp, "@ legend 0.75, 0.8\n");
- fprintf(fp, "@ legend string 0 \"Total/mean\"\n");
- fprintf(fp, "@ legend string 1 \"Total\"\n");
- fprintf(fp, "@ legend string 2 \"Mean\"\n");
- fprintf(fp, "@ legend string 3 \"# atoms\"\n");
- fprintf(fp, "@ legend string 4 \"Mean/# atoms\"\n");
- fprintf(fp, "#%3s %8s %3s %8s %3s %8s\n",
- "res", "ratio", "tot", "mean", "natm", "mean/atm");
+ sprintf(legend[0], "Total/mean");
+ sprintf(legend[1], "Total");
+ sprintf(legend[2], "Mean");
+ sprintf(legend[3], "# atoms");
+ sprintf(legend[4], "Mean/# atoms");
+ xvgr_legend(fp, 5, (const char**)legend, oenv);
for (i = 0; (i < nres); i++)
{
if (mean_n[i] == 0)
#include "gromacs/fileio/trxio.h"
#include "rmpbc.h"
#include "gmx_ana.h"
+#include "names.h"
-static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
+static void periodic_dist(int ePBC,
+ matrix box, rvec x[], int n, atom_id index[],
real *rmin, real *rmax, int *min_ind)
{
-#define NSHIFT 26
- int sx, sy, sz, i, j, s;
+#define NSHIFT_MAX 26
+ int nsz, nshift, sx, sy, sz, i, j, s;
real sqr_box, r2min, r2max, r2;
- rvec shift[NSHIFT], d0, d;
+ rvec shift[NSHIFT_MAX], d0, d;
- sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));
+ sqr_box = min(norm2(box[XX]), norm2(box[YY]));
+ if (ePBC == epbcXYZ)
+ {
+ sqr_box = min(sqr_box, norm2(box[ZZ]));
+ nsz = 1;
+ }
+ else if (ePBC == epbcXY)
+ {
+ nsz = 0;
+ }
+ else
+ {
+ gmx_fatal(FARGS, "pbc = %s is not supported by g_mindist",
+ epbc_names[ePBC]);
+ nsz = 0; /* Keep compilers quiet */
+ }
- s = 0;
- for (sz = -1; sz <= 1; sz++)
+ nshift = 0;
+ for (sz = -nsz; sz <= nsz; sz++)
{
for (sy = -1; sy <= 1; sy++)
{
{
for (i = 0; i < DIM; i++)
{
- shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
+ shift[nshift][i] =
+ sx*box[XX][i] + sy*box[YY][i] + sz*box[ZZ][i];
}
- s++;
+ nshift++;
}
}
}
{
r2max = r2;
}
- for (s = 0; s < NSHIFT; s++)
+ for (s = 0; s < nshift; s++)
{
rvec_add(d0, shift[s], d);
r2 = norm2(d);
gmx_rmpbc(gpbc, natoms, box, x);
}
- periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
+ periodic_dist(ePBC, box, x, n, index, &rmin, &rmax, ind_min);
if (rmin < rmint)
{
rmint = rmin;
}
if (bSplit && !bFirst && fabs(t/output_env_get_time_factor(oenv)) < 1e-5)
{
- fprintf(out, "&\n");
+ fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2]));
{
if (bSplit && !bFirst && fabs(t/output_env_get_time_factor(oenv)) < 1e-5)
{
- fprintf(dist, "&\n");
+ fprintf(dist, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
if (num)
{
- fprintf(num, "&\n");
+ fprintf(num, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
if (atm)
{
- fprintf(atm, "&\n");
+ fprintf(atm, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
}
fprintf(dist, "%12e", output_env_conv_time(oenv, t));
/* Handle printing of internal distances. */
if (outi)
{
- fprintf(outi, "@ xaxes scale Logarithmic\n");
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(outi, "@ xaxes scale Logarithmic\n");
+ }
ymax = -1;
ymin = 1e300;
j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
xvgr_line_props(out, 0, elNone, ecFrank, oenv);
xvgr_view(out, 0.2, 0.2, 0.8, 0.8, oenv);
xvgr_world(out, -180, -180, 180, 180, oenv);
- fprintf(out, "@ xaxis tick on\n@ xaxis tick major 60\n@ xaxis tick minor 30\n");
- fprintf(out, "@ yaxis tick on\n@ yaxis tick major 60\n@ yaxis tick minor 30\n");
- fprintf(out, "@ s0 symbol 2\n@ s0 symbol size 0.4\n@ s0 symbol fill 1\n");
-
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "@ xaxis tick on\n@ xaxis tick major 60\n@ xaxis tick minor 30\n");
+ fprintf(out, "@ yaxis tick on\n@ yaxis tick major 60\n@ yaxis tick minor 30\n");
+ fprintf(out, "@ s0 symbol 2\n@ s0 symbol size 0.4\n@ s0 symbol fill 1\n");
+ }
j = 0;
do
{
bPrev = (prev > 0);
if (bPrev)
{
+ fprintf(stderr, "WARNING: using option -prev with large trajectories will\n"
+ " require a lot of memory and could lead to crashes\n");
prev = abs(prev);
if (freq != 1)
{
if (bSplit && i > 0 &&
fabs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
{
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]);
for (j = 0; (j < nrms); j++)
{
if (bSplit && i > 0 && fabs(time[i]) < 1e-5)
{
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
fprintf(fp, "%12.7f", time[i]);
for (j = 0; (j < nrms); j++)
tcafc[kc][i] /= tcafc[kc][0];
fprintf(fp_cub, "%g %g\n", i*dt, tcafc[kc][i]);
}
- fprintf(fp_cub, "&\n");
+ fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
tcafc[kc][0] = 1.0;
}
}
fp_vk = xvgropen(fn_vk, "Fits", "k (nm\\S-1\\N)",
"\\8h\\4 (10\\S-3\\N kg m\\S-1\\N s\\S-1\\N)", oenv);
- fprintf(fp_vk, "@ s0 symbol 2\n");
- fprintf(fp_vk, "@ s0 symbol color 1\n");
- fprintf(fp_vk, "@ s0 linestyle 0\n");
- if (fn_cub)
+ if (output_env_get_print_xvgr_codes(oenv))
{
- fprintf(fp_vk, "@ s1 symbol 3\n");
- fprintf(fp_vk, "@ s1 symbol color 2\n");
+ fprintf(fp_vk, "@ s0 symbol 2\n");
+ fprintf(fp_vk, "@ s0 symbol color 1\n");
+ fprintf(fp_vk, "@ s0 linestyle 0\n");
+ if (fn_cub)
+ {
+ fprintf(fp_vk, "@ s1 symbol 3\n");
+ fprintf(fp_vk, "@ s1 symbol color 2\n");
+ }
}
fp = xvgropen(fn_tcf, "TCAF Fits", "Time (ps)", "", oenv);
for (k = 0; k < nk; k++)
{
fprintf(fp, "%g %g\n", i*dt, fit_function(effnVAC, fitparms, i*dt));
}
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
gmx_ffclose(fp);
do_view(oenv, fn_tcf, "-nxy");
if (fn_cub)
{
fprintf(stdout, "Averaged over k-vectors:\n");
- fprintf(fp_vk, "&\n");
+ fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
for (k = 0; k < nkc; k++)
{
tcafc[k][0] = 1.0;
{
fprintf(fp_cub, "%g %g\n", i*dt, fit_function(effnVAC, fitparms, i*dt));
}
- fprintf(fp_cub, "&\n");
+ fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
- fprintf(fp_vk, "&\n");
+ fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
gmx_ffclose(fp_cub);
do_view(oenv, fn_cub, "-nxy");
}
if (orfile)
{
fp = xvgropen(orfile, "Van Hove function", "r (nm)", "G (nm\\S-1\\N)", oenv);
- fprintf(fp, "@ subtitle \"for particles in group %s\"\n", grpname);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ subtitle \"for particles in group %s\"\n", grpname);
+ }
snew(legend, nr);
for (fbin = 0; fbin < nr; fbin++)
{
{
sprintf(buf, "Probability of moving less than %g nm", rint);
fp = xvgropen(otfile, buf, "t (ps)", "", oenv);
- fprintf(fp, "@ subtitle \"for particles in group %s\"\n", grpname);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ subtitle \"for particles in group %s\"\n", grpname);
+ }
for (f = 0; f <= ftmax; f++)
{
fprintf(fp, "%g %g\n", f*dt, (real)pt[f]/(tcount[f]*isize));
bsProfiles_av2[i] += tmp*tmp;
fprintf(fp, "%e\t%e\n", (i+0.5)*opt->dz+opt->min, tmp);
}
- fprintf(fp, "&\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
}
gmx_ffclose(fp);
/* write average and stddev */
fp = xvgropen(fnres, "Average and stddev from bootstrapping", "z", ylabel, opt->oenv);
- fprintf(fp, "@TYPE xydy\n");
+ if (output_env_get_print_xvgr_codes(opt->oenv))
+ {
+ fprintf(fp, "@TYPE xydy\n");
+ }
for (i = 0; i < opt->bins; i++)
{
bsProfiles_av [i] /= opt->nBootStrap;
{
fprintf(fpcorr, "%g %g\n", k*dt, corr[k]);
}
- fprintf(fpcorr, "&\n");
+ fprintf(fpcorr, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
}
/* esimate integrated correlation time, fitting is too unstable */
/* plot IACT along reaction coordinate */
fp = xvgropen(fn, "Integrated autocorrelation times", "z", "IACT [ps]", opt->oenv);
- fprintf(fp, "@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
- fprintf(fp, "# WIN tau(gr1) tau(gr2) ...\n");
- for (i = 0; i < nwins; i++)
+ if (output_env_get_print_xvgr_codes(opt->oenv))
{
- fprintf(fp, "# %3d ", i);
- for (ig = 0; ig < window[i].nPull; ig++)
+ fprintf(fp, "@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
+ fprintf(fp, "# WIN tau(gr1) tau(gr2) ...\n");
+ for (i = 0; i < nwins; i++)
{
- fprintf(fp, " %11g", window[i].tau[ig]);
+ fprintf(fp, "# %3d ", i);
+ for (ig = 0; ig < window[i].nPull; ig++)
+ {
+ fprintf(fp, " %11g", window[i].tau[ig]);
+ }
+ fprintf(fp, "\n");
}
- fprintf(fp, "\n");
}
for (i = 0; i < nwins; i++)
{
opt->sigSmoothIact);
/* smooth IACT along reaction coordinate and overwrite g=1+2tau */
smoothIact(window, nwins, opt);
- fprintf(fp, "&\n@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
- fprintf(fp, "@ s1 symbol color 2\n");
+ fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
+ if (output_env_get_print_xvgr_codes(opt->oenv))
+ {
+ fprintf(fp, "@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
+ fprintf(fp, "@ s1 symbol color 2\n");
+ }
for (i = 0; i < nwins; i++)
{
for (ig = 0; ig < window[i].nPull; ig++)
* a shell connected to a dummy with spring constant that differ in the
* three spatial dimensions in the molecular frame.
*/
- int i, m, aO, aH1, aH2, aD, aS, type, type0;
+ int i, m, aO, aH1, aH2, aD, aS, type, type0, ki;
+ ivec dt;
rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj;
#ifdef DEBUG
rvec df;
aS = forceatoms[i+5];
/* Compute vectors describing the water frame */
- rvec_sub(x[aH1], x[aO], dOH1);
- rvec_sub(x[aH2], x[aO], dOH2);
- rvec_sub(x[aH2], x[aH1], dHH);
- rvec_sub(x[aD], x[aO], dOD);
- rvec_sub(x[aS], x[aD], dDS);
+ pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1);
+ pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2);
+ pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH);
+ pbc_rvec_sub(pbc, x[aD], x[aO], dOD);
+ ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS);
cprod(dOH1, dOH2, nW);
/* Compute inverse length of normal vector
kdx[YY] = kk[YY]*dx[YY];
kdx[ZZ] = kk[ZZ]*dx[ZZ];
vtot += iprod(dx, kdx);
+
+ if (g)
+ {
+ ivec_sub(SHIFT_IVEC(g, aS), SHIFT_IVEC(g, aD), dt);
+ ki = IVEC2IS(dt);
+ }
+
for (m = 0; (m < DIM); m++)
{
/* This is a tensor operation but written out for speed */
#ifdef DEBUG
df[m] = fij;
#endif
- f[aS][m] += fij;
- f[aD][m] -= fij;
+ f[aS][m] += fij;
+ f[aD][m] -= fij;
+ fshift[ki][m] += fij;
+ fshift[CENTRAL][m] -= fij;
}
#ifdef DEBUG
if (debug)
vtot += 0.5*kk*dx[m]*dx[m];
*dvdlambda +=
0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
- -fm*dpdl[m];
+ + fm*dpdl[m];
/* Here we correct for the pbc_dx which included rdist */
if (bForceValid)
ivec dd_nc, ivec dd_nc_f)
{
int npp;
- gmx_bool mm;
-
- mm = FALSE;
+ gmx_bool mm = FALSE;
+ gmx_bool patchlevel_differs = FALSE;
+ gmx_bool version_differs = FALSE;
check_string(fplog, "Version", gmx_version(), version, &mm);
+ patchlevel_differs = mm;
+
+ if (patchlevel_differs)
+ {
+ /* Gromacs should be able to continue from checkpoints between
+ * different patch level versions, but we do not guarantee
+ * compatibility between different major/minor versions - check this.
+ */
+ int gmx_major, gmx_minor;
+ int cpt_major, cpt_minor;
+ sscanf(gmx_version(), "VERSION %d.%d", &gmx_major, &gmx_minor);
+ sscanf(version, "VERSION %d.%d", &cpt_major, &cpt_minor);
+ version_differs = (gmx_major != cpt_major || gmx_minor != cpt_minor);
+ }
+
check_string(fplog, "Build time", BUILD_TIME, btime, &mm);
check_string(fplog, "Build user", BUILD_USER, buser, &mm);
check_string(fplog, "Build host", BUILD_HOST, bhost, &mm);
if (mm)
{
- fprintf(stderr,
- "Gromacs binary or parallel settings not identical to previous run.\n"
- "Continuation is exact, but is not guaranteed to be binary identical%s.\n\n",
- fplog ? ",\n see the log file for details" : "");
+ const char msg_version_difference[] =
+ "The current Gromacs major & minor version are not identical to those that\n"
+ "generated the checkpoint file. In principle Gromacs does not support\n"
+ "continuation from checkpoints between different versions, so we advise\n"
+ "against this. If you still want to try your luck we recommend that you use\n"
+ "the -noappend flag to keep your output files from the two versions separate.\n"
+ "This might also work around errors where the output fields in the energy\n"
+ "file have changed between the different major & minor versions.\n";
- if (fplog)
+ const char msg_mismatch_notice[] =
+ "Gromacs patchlevel, binary or parallel settings differ from previous run.\n"
+ "Continuation is exact, but not guaranteed to be binary identical.\n";
+
+ const char msg_logdetails[] =
+ "See the log file for details.\n";
+
+ if (version_differs)
{
- fprintf(fplog,
- "Gromacs binary or parallel settings not identical to previous run.\n"
- "Continuation is exact, but is not guaranteed to be binary identical.\n\n");
+ fprintf(stderr, "%s%s\n", msg_version_difference, fplog ? msg_logdetails : "");
+
+ if (fplog)
+ {
+ fprintf(fplog, "%s\n", msg_version_difference);
+ }
+ }
+ else
+ {
+ /* Major & minor versions match at least, but something is different. */
+ fprintf(stderr, "%s%s\n", msg_mismatch_notice, fplog ? msg_logdetails : "");
+ if (fplog)
+ {
+ fprintf(fplog, "%s\n", msg_mismatch_notice);
+ }
}
}
}
}
fprintf(stderr, "\n");
- gmx_fatal(FARGS, "File appending requested, but only %d of the %d output files are present", nexist, nfiles);
+ gmx_fatal(FARGS, "File appending requested, but %d of the %d output files are not present or are named differently", nfiles-nexist, nfiles);
}
}
#include "calc_verletbuf.h"
#include "tomorse.h"
#include "gromacs/imd/imd.h"
+#include "gromacs/utility/cstringutil.h"
static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[])
}
}
+static void check_shells_inputrec(gmx_mtop_t *mtop,
+ t_inputrec *ir,
+ warninp_t wi)
+{
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ int a, nshells = 0;
+ char warn_buf[STRLEN];
+
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
+ {
+ if (atom->ptype == eptShell ||
+ atom->ptype == eptBond)
+ {
+ nshells++;
+ }
+ }
+ if (IR_TWINRANGE(*ir) && (nshells > 0))
+ {
+ snprintf(warn_buf, STRLEN,
+ "The combination of using shells and a twin-range cut-off is not supported");
+ warning_error(wi, warn_buf);
+ }
+ if ((nshells > 0) && (ir->nstcalcenergy != 1))
+ {
+ set_warning_line(wi, "unknown", -1);
+ snprintf(warn_buf, STRLEN,
+ "There are %d shells, changing nstcalcenergy from %d to 1",
+ nshells, ir->nstcalcenergy);
+ ir->nstcalcenergy = 1;
+ warning(wi, warn_buf);
+ }
+}
+
static gmx_bool nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype)
{
int nint, mb;
check_vel(sys, state.v);
}
+ /* check for shells and inpurecs */
+ check_shells_inputrec(sys, ir, wi);
+
/* check masses */
check_mol(sys, wi);
{
pgrp = &pull->group[g];
+ if (strcmp(pgnames[g], "") == 0)
+ {
+ gmx_fatal(FARGS, "Group pull_group%d required by grompp was undefined.", g);
+ }
+
ig = search_string(pgnames[g], grps->nr, gnames);
pgrp->nat = grps->index[ig+1] - grps->index[ig];
break;
case eCOMB_ARITHMETIC:
- /* c0 and c1 are epsilon and sigma */
+ /* c0 and c1 are sigma and epsilon */
for (i = k = 0; (i < nr); i++)
{
for (j = 0; (j < nr); j++, k++)
cj0 = get_atomtype_nbparam(j, 0, atype);
ci1 = get_atomtype_nbparam(i, 1, atype);
cj1 = get_atomtype_nbparam(j, 1, atype);
- plist->param[k].c[0] = (ci0+cj0)*0.5;
+ plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5;
+ /* Negative sigma signals that c6 should be set to zero later,
+ * so we need to propagate that through the combination rules.
+ */
+ if (ci0 < 0 || cj0 < 0)
+ {
+ plist->param[k].c[0] *= -1;
+ }
plist->param[k].c[1] = sqrt(ci1*cj1);
}
}
break;
case eCOMB_GEOM_SIG_EPS:
- /* c0 and c1 are epsilon and sigma */
+ /* c0 and c1 are sigma and epsilon */
for (i = k = 0; (i < nr); i++)
{
for (j = 0; (j < nr); j++, k++)
cj0 = get_atomtype_nbparam(j, 0, atype);
ci1 = get_atomtype_nbparam(i, 1, atype);
cj1 = get_atomtype_nbparam(j, 1, atype);
- plist->param[k].c[0] = sqrt(ci0*cj0);
+ plist->param[k].c[0] = sqrt(fabs(ci0*cj0));
+ /* Negative sigma signals that c6 should be set to zero later,
+ * so we need to propagate that through the combination rules.
+ */
+ if (ci0 < 0 || cj0 < 0)
+ {
+ plist->param[k].c[0] *= -1;
+ }
plist->param[k].c[1] = sqrt(ci1*cj1);
}
}
/* return time conversion factor from ps (i.e. 1e-3 for ps->ns) */
real output_env_get_time_invfactor(const output_env_t oenv);
-/* return inverse time conversion factor from ps (i.e. 1e3 for ps->ns) */
+/* return inverse time conversion factor to ps (i.e. 1e3 for ns->ps) */
real output_env_conv_time(const output_env_t oenv, real time);
/* return converted time */
* If x!=NULL, the shells are predict for the global coordinates x.
*/
gmx_shellfc_t init_shell_flexcon(FILE *fplog,
- gmx_bool bCutoffSchemeIsVerlet,
gmx_mtop_t *mtop, int nflexcon,
rvec *x);
real *table_V,
real *table_FDV0,
int ntab,
- real dx,
+ double dx,
real beta,
real_space_grid_contribution_computer v_lr);
/* Fill tables of ntab points with spacing dr with the ewald long-range
/* Now loop over all the thread data blocks that contribute
* to the grid region we (our thread) are operating on.
*/
- /* Note that ffy_nx/y is equal to the number of grid points
+ /* Note that fft_nx/y is equal to the number of grid points
* between the first point of our node grid and the one of the next node.
*/
for (sx = 0; sx >= -pmegrids->nthread_comm[XX]; sx--)
}
pmegrid_g = &pmegrids->grid_th[fx*pmegrids->nc[YY]*pmegrids->nc[ZZ]];
ox += pmegrid_g->offset[XX];
- if (!bCommX)
- {
- tx1 = min(ox + pmegrid_g->n[XX], ne[XX]);
- }
- else
- {
- tx1 = min(ox + pmegrid_g->n[XX], pme->pme_order);
- }
+ /* Determine the end of our part of the source grid */
+ tx1 = min(ox + pmegrid_g->n[XX], ne[XX]);
for (sy = 0; sy >= -pmegrids->nthread_comm[YY]; sy--)
{
}
pmegrid_g = &pmegrids->grid_th[fy*pmegrids->nc[ZZ]];
oy += pmegrid_g->offset[YY];
- if (!bCommY)
- {
- ty1 = min(oy + pmegrid_g->n[YY], ne[YY]);
- }
- else
- {
- ty1 = min(oy + pmegrid_g->n[YY], pme->pme_order);
- }
+ /* Determine the end of our part of the source grid */
+ ty1 = min(oy + pmegrid_g->n[YY], ne[YY]);
for (sz = 0; sz >= -pmegrids->nthread_comm[ZZ]; sz--)
{
}
gmx_shellfc_t init_shell_flexcon(FILE *fplog,
- gmx_bool bCutoffSchemeIsVerlet,
gmx_mtop_t *mtop, int nflexcon,
rvec *x)
{
return NULL;
}
- if (bCutoffSchemeIsVerlet)
- {
- gmx_fatal(FARGS, "The shell code does not work with the Verlet cut-off scheme.\n");
- }
-
snew(shfc, 1);
shfc->nflexcon = nflexcon;
{
fprintf(fplog, "\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n");
}
+ /* Prediction improves performance, so we should implement either:
+ * 1. communication for the atoms needed for prediction
+ * 2. prediction using the velocities of shells; currently the
+ * shell velocities are zeroed, it's a bit tricky to keep
+ * track of the shell displacements and thus the velocity.
+ */
shfc->bPredict = FALSE;
}
}
force[i] = shfc->f[i];
}
- /* When we had particle decomposition, this code only worked with
- * PD when all particles involved with each shell were in the same
- * charge group. Not sure if this is still relevant. */
if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr))
{
/* This is the only time where the coordinates are used
* before do_force is called, which normally puts all
* charge groups in the box.
*/
- cg0 = 0;
- cg1 = top->cgs.nr;
- put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box,
- &(top->cgs), state->x, fr->cg_cm);
+ if (inputrec->cutoff_scheme == ecutsVERLET)
+ {
+ put_atoms_in_box_omp(fr->ePBC, state->box, md->homenr, state->x);
+ }
+ else
+ {
+ cg0 = 0;
+ cg1 = top->cgs.nr;
+ put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box,
+ &(top->cgs), state->x, fr->cg_cm);
+ }
+
if (graph)
{
mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
}
}
- /* After this all coordinate arrays will contain whole molecules */
+ /* After this all coordinate arrays will contain whole charge groups */
if (graph)
{
shift_self(graph, state->box, state->x);
real *table_v,
real *table_fdv0,
int ntab,
- real dx,
+ double dx,
real beta,
real_space_grid_contribution_computer v_lr)
{
}
if (md->nPerturbed && md->bPerturbed[n])
{
- dekindl += 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
+ /* The minus sign here might be confusing.
+ * The kinetic contribution from dH/dl doesn't come from
+ * d m(l)/2 v^2 / dl, but rather from d p^2/2m(l) / dl,
+ * where p are the momenta. The difference is only a minus sign.
+ */
+ dekindl -= 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
}
}
ekind->dekindl = dekindl;
}
else
{
- lamfac = 0;
+ lamfac = lambda;
type = md->typeB;
}
}
for (c = 0; c < pull->ncoord; c++)
{
pcrd = &pull->coord[c];
+ ref = pcrd->init + pcrd->rate*t;
low_get_pull_coord_dr(pull, pcrd, pbc, t,
rnew[pcrd->group[1]],
}
if (cmp_bool(fp, "bF", -1, fr1->bF, fr2->bF))
{
- cmp_rvecs_rmstol(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, ftol, abstol);
+ if (bRMSD)
+ {
+ cmp_rvecs(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, bRMSD, ftol, abstol);
+ }
+ else
+ {
+ cmp_rvecs_rmstol(fp, "f", min(fr1->natoms, fr2->natoms), fr1->f, fr2->f, ftol, abstol);
+ }
}
if (cmp_bool(fp, "bBox", -1, fr1->bBox, fr2->bBox))
{
debug_gmx();
/* Check for polarizable models and flexible constraints */
- shellfc = init_shell_flexcon(fplog, fr->cutoff_scheme == ecutsVERLET,
+ shellfc = init_shell_flexcon(fplog,
top_global, n_flexible_constraints(constr),
(ir->bContinuation ||
(DOMAINDECOMP(cr) && !MASTER(cr))) ?
NULL : state_global->x);
-
+ if (shellfc && ir->nstcalcenergy != 1)
+ {
+ gmx_fatal(FARGS, "You have nstcalcenergy set to a value (%d) that is different from 1.\nThis is not supported in combinations with shell particles.\nPlease make a new tpr file.", ir->nstcalcenergy);
+ }
+ if (shellfc && DOMAINDECOMP(cr))
+ {
+ gmx_fatal(FARGS, "Shell particles are not implemented with domain decomposition, use a single rank");
+ }
if (shellfc && ir->eI == eiNM)
{
/* Currently shells don't work with Normal Modes */
"[PAR]",
"When PME is used with domain decomposition, separate nodes can",
"be assigned to do only the PME mesh calculation;",
- "this is computationally more efficient starting at about 12 nodes.",
+ "this is computationally more efficient starting at about 12 nodes",
+ "or even fewer when OpenMP parallelization is used.",
"The number of PME nodes is set with option [TT]-npme[tt],",
"this can not be more than half of the nodes.",
"By default [TT]mdrun[tt] makes a guess for the number of PME",
- "nodes when the number of nodes is larger than 11 or performance wise",
- "not compatible with the PME grid x dimension.",
- "But the user should optimize npme. Performance statistics on this issue",
+ "nodes when the number of nodes is larger than 16. With GPUs,",
+ "PME nodes are not selected automatically, since the optimal setup",
+ "depends very much on the details of the hardware.",
+ "In all cases you might gain performance by optimizing [TT]-npme[tt].",
+ "Performance statistics on this issue",
"are written at the end of the log file.",
"For good load balancing at high parallelization, the PME grid x and y",
"dimensions should be divisible by the number of PME nodes",
{
if (re->q[re->type][re->ind[j]] < re->q[re->type][re->ind[i]])
{
+ /* Unordered replicas are supposed to work, but there
+ * is still an issues somewhere.
+ * Note that at this point still re->ind[i]=i.
+ */
+ gmx_fatal(FARGS, "Replicas with indices %d < %d have %ss %g > %g, please order your replicas on increasing %s",
+ i, j,
+ erename[re->type],
+ re->q[re->type][i], re->q[re->type][j],
+ erename[re->type]);
+
k = re->ind[i];
re->ind[i] = re->ind[j];
re->ind[j] = k;
dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
if (bPrint)
{
- fprintf(fplog, " dpV = %10.3e d = %10.3e\nb", dpV, delta + dpV);
+ fprintf(fplog, " dpV = %10.3e d = %10.3e\n", dpV, delta + dpV);
}
delta += dpV;
}
gmx_rng_t rng;
bMultiEx = (re->nex > 1); /* multiple exchanges at each state */
- fprintf(fplog, "Replica exchange at step " "%"GMX_PRId64 " time %g\n", step, time);
+ fprintf(fplog, "Replica exchange at step " "%"GMX_PRId64 " time %.5f\n", step, time);
if (re->bNPT)
{
}
static void
-prepare_to_do_exchange(FILE *fplog,
- const int *destinations,
- const int replica_id,
- const int nrepl,
- int *maxswap,
- int **order,
- int **cyclic,
- int *incycle,
- gmx_bool *bThisReplicaExchanged)
+prepare_to_do_exchange(FILE *fplog,
+ struct gmx_repl_ex *re,
+ const int replica_id,
+ int *maxswap,
+ gmx_bool *bThisReplicaExchanged)
{
int i, j;
/* Hold the cyclic decomposition of the (multiple) replica
gmx_bool bAnyReplicaExchanged = FALSE;
*bThisReplicaExchanged = FALSE;
- for (i = 0; i < nrepl; i++)
+ for (i = 0; i < re->nrepl; i++)
{
- if (destinations[i] != i)
+ if (re->destinations[i] != re->ind[i])
{
/* only mark as exchanged if the index has been shuffled */
bAnyReplicaExchanged = TRUE;
if (bAnyReplicaExchanged)
{
/* reinitialize the placeholder arrays */
- for (i = 0; i < nrepl; i++)
+ for (i = 0; i < re->nrepl; i++)
{
- for (j = 0; j < nrepl; j++)
+ for (j = 0; j < re->nrepl; j++)
{
- cyclic[i][j] = -1;
- order[i][j] = -1;
+ re->cyclic[i][j] = -1;
+ re->order[i][j] = -1;
}
}
/* Identify the cyclic decomposition of the permutation (very
* fast if neighbor replica exchange). */
- cyclic_decomposition(destinations, cyclic, incycle, nrepl, maxswap);
+ cyclic_decomposition(re->destinations, re->cyclic, re->incycle, re->nrepl, maxswap);
/* Now translate the decomposition into a replica exchange
* order at each step. */
- compute_exchange_order(fplog, cyclic, order, nrepl, *maxswap);
+ compute_exchange_order(fplog, re->cyclic, re->order, re->nrepl, *maxswap);
/* Did this replica do any exchange at any point? */
for (j = 0; j < *maxswap; j++)
{
- if (replica_id != order[replica_id][j])
+ if (replica_id != re->order[replica_id][j])
{
*bThisReplicaExchanged = TRUE;
break;
{
replica_id = re->repl;
test_for_replica_exchange(fplog, cr->ms, re, enerd, det(state_local->box), step, time);
- prepare_to_do_exchange(fplog, re->destinations, replica_id, re->nrepl, &maxswap,
- re->order, re->cyclic, re->incycle, &bThisReplicaExchanged);
+ prepare_to_do_exchange(fplog, re, replica_id, &maxswap, &bThisReplicaExchanged);
}
/* Do intra-simulation broadcast so all processors belonging to
* each simulation know whether they need to participate in
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff