# NOTE: when releasing the "-dev" suffix needs to be stripped off!
# REGRESSIONTEST_VERSION and REGRESSIONTEST_BRANCH should always be
# defined.
-set(PROJECT_VERSION "5.0-rc1-dev")
+set(PROJECT_VERSION "5.1-dev")
# If this is a released tarball, "-dev" will not be present in
# PROJECT_VERSION, and REGRESSIONTEST_VERSION specifies the version
# number of the regressiontest tarball against which the code tarball
# PROJECT_VERSION, and REGRESSIONTEST_BRANCH specifies the name of the
# gerrit.gromacs.org branch whose HEAD can test this code, *if* this
# code contains all recent fixes from the corresponding code branch.
-set(REGRESSIONTEST_BRANCH "refs/heads/release-5-0")
+set(REGRESSIONTEST_BRANCH "refs/heads/master")
set(CUSTOM_VERSION_STRING ""
CACHE STRING "Custom version string (if empty, use hard-coded default)")
set(LIBRARY_VERSION ${LIBRARY_SOVERSION}.0.0)
# It is a bit irritating, but this has to be set separately for now!
SET(CPACK_PACKAGE_VERSION_MAJOR "5")
-SET(CPACK_PACKAGE_VERSION_MINOR "0")
+SET(CPACK_PACKAGE_VERSION_MINOR "1")
#SET(CPACK_PACKAGE_VERSION_PATCH "0")
# The numerical gromacs version. It is 40600 for 4.6.0.
mark_as_advanced(GMX_INSTALL_PREFIX)
include(gmxBuildTypeReference)
+include(gmxBuildTypeProfile)
include(gmxBuildTypeTSAN)
include(gmxBuildTypeASAN)
include(gmxBuildTypeReleaseWithAssert)
if(NOT CMAKE_BUILD_TYPE)
- set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel Reference RelWithAssert." FORCE)
+ set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel Reference RelWithAssert Profile." FORCE)
# There's no need to offer a user the choice of ThreadSanitizer
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release"
- "MinSizeRel" "RelWithDebInfo" "Reference" "RelWithAssert")
+ "MinSizeRel" "RelWithDebInfo" "Reference" "RelWithAssert" "Profile")
endif()
if(CMAKE_CONFIGURATION_TYPES)
# Add appropriate GROMACS-specific build types for the Visual
"List of configuration types"
FORCE)
endif()
-set(build_types_with_explicit_flags RELEASE DEBUG RELWITHDEBUGINFO RELWITHASSERT MINSIZEREL)
+set(build_types_with_explicit_flags RELEASE DEBUG RELWITHDEBUGINFO RELWITHASSERT MINSIZEREL PROFILE)
enable_language(C)
enable_language(CXX)
# 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}")
#######################
## uninstall target
#######################
- CONFIGURE_FILE(
- "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_uninstall.cmake.in"
- "${CMAKE_CURRENT_BINARY_DIR}/cmake/cmake_uninstall.cmake"
- IMMEDIATE @ONLY)
+ CONFIGURE_FILE( "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_uninstall.cmake.in"
+ "${CMAKE_CURRENT_BINARY_DIR}/cmake/cmake_uninstall.cmake"
+ IMMEDIATE @ONLY)
###########################
ADD_CUSTOM_TARGET(uninstall
"${CMAKE_COMMAND}" -P
"${CMAKE_CURRENT_BINARY_DIR}/cmake/cmake_uninstall.cmake")
###########################
+ set_directory_properties(PROPERTIES
+ ADDITIONAL_MAKE_CLEAN_FILES "install_manifest.txt")
########################################################################
# Manual #
# be set up elsewhere and passed to this function, but it is
# inconvenient in CMake to pass more than one list, and such a
# list is only used here.
- foreach(build_type RELWITHDEBUGINFO RELWITHASSERT MINSIZEREL)
+ foreach(build_type RELWITHDEBUGINFO RELWITHASSERT MINSIZEREL PROFILE)
set(GMXC_${language}FLAGS_${build_type} "${GMXC_${language}FLAGS_RELEASE}")
endforeach()
# Copy the flags that are only used by the real Release build
# Since 4.8 on by default. For previous version disabling is a no-op. Only disabling for Release because with assert
# the warnings are OK.
GMX_TEST_CFLAG(CFLAGS_WARN_REL "-Wno-array-bounds" GMXC_CFLAGS_RELEASE_ONLY)
- # Since gcc 4.8 strict - false postives with old gmx_fatal. TODO: Remove in master
- GMX_TEST_CFLAG(CFLAGS_WARN_UNINIT "-Wno-maybe-uninitialized" GMXC_CFLAGS)
# new in gcc 4.5
GMX_TEST_CFLAG(CFLAGS_EXCESS_PREC "-fexcess-precision=fast" GMXC_CFLAGS_RELEASE)
- GMX_TEST_CFLAG(CFLAGS_COPT "-fomit-frame-pointer -funroll-all-loops"
+ GMX_TEST_CFLAG(CFLAGS_COPT "-funroll-all-loops"
GMXC_CFLAGS_RELEASE)
GMX_TEST_CFLAG(CFLAGS_NOINLINE "-fno-inline" GMXC_CFLAGS_DEBUG)
endif()
GMX_TEST_CFLAG(CXXFLAGS_WARN_REL "-Wno-array-bounds" GMXC_CXXFLAGS_RELEASE_ONLY)
# new in gcc 4.5
GMX_TEST_CXXFLAG(CXXFLAGS_EXCESS_PREC "-fexcess-precision=fast" GMXC_CXXFLAGS_RELEASE)
- GMX_TEST_CXXFLAG(CXXFLAGS_COPT "-fomit-frame-pointer -funroll-all-loops"
+ GMX_TEST_CXXFLAG(CXXFLAGS_COPT "-funroll-all-loops"
GMXC_CXXFLAGS_RELEASE)
GMX_TEST_CXXFLAG(CXXFLAGS_NOINLINE "-fno-inline" GMXC_CXXFLAGS_DEBUG)
endif()
# unreferenced local variable (only C)
# conversion from 'size_t' to 'int', possible loss of data
# conversion from 'const char*' to 'void*', different 'const' qualifiers (only C)
- GMX_TEST_CFLAG(CFLAGS_WARN "/wd4800 /wd4355 /wd4996 /wd4305 /wd4244 /wd4101 /wd4267 /wd4090" GMXC_CFLAGS)
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/wd4800 /wd4355 /wd4996 /wd4305 /wd4244 /wd4267" GMXC_CXXFLAGS)
+ if(NOT CMAKE_CONFIGURATION_TYPES)
+ GMX_TEST_CFLAG(CFLAGS_WARN "/wd4800 /wd4355 /wd4996 /wd4305 /wd4244 /wd4101 /wd4267 /wd4090" GMXC_CFLAGS)
+ GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/wd4800 /wd4355 /wd4996 /wd4305 /wd4244 /wd4267" GMXC_CXXFLAGS)
+ else() #Projects only use the C++ flags
+ GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/wd4800 /wd4355 /wd4996 /wd4305 /wd4244 /wd4101 /wd4267 /wd4090" GMXC_CXXFLAGS)
+ endif()
endif()
if (CMAKE_C_COMPILER_ID MATCHES "Clang")
@Article{PSmith93c,
author = {P. E. Smith and W. F. van Gunsteren},
title = {{The Viscosity of SPC and SPC/E Water}},
- journal = BTcpc,
+ journal = BTcpl,
year = 1993,
volume = 215,
pages = {315--318}
pages = {2887--2899},
year = 2013
}
+
+ @article {Mac2000,
+ author = {MacKerell, Alexander D. and Banavali, Nilesh K.},
+ title = {All-atom empirical force field for nucleic acids: {II.} Application to molecular dynamics simulations of {DNA} and {RNA} in solution},
+ journal = BTjcc,
+ volume = {21},
+ number = {2},
+ publisher = {John Wiley & Sons, Inc.},
+ issn = {1096-987X},
+ url = {http://dx.doi.org/10.1002/(SICI)1096-987X(20000130)21:2<105::AID-JCC3>3.0.CO;2-P},
+ doi = {10.1002/(SICI)1096-987X(20000130)21:2<105::AID-JCC3>3.0.CO;2-P},
+ pages = {105--120},
+ keywords = {CHARMM, force field, molecular dynamics, parametrization, DNA, RNA},
+ year = {2000},
+ }
add_subdirectory(random)
add_subdirectory(onlinehelp)
add_subdirectory(options)
+add_subdirectory(pbcutil)
add_subdirectory(timing)
+add_subdirectory(topology)
add_subdirectory(utility)
add_subdirectory(fileio)
add_subdirectory(swap)
endif()
set_source_files_properties(selection/scanner.cpp PROPERTIES COMPILE_FLAGS "${_scanner_cpp_compiler_flags}")
- target_link_libraries(libgromacs ${GMX_GPU_LIBRARIES}
+ target_link_libraries(libgromacs
+ ${EXTRAE_LIBRARIES}
+ ${GMX_GPU_LIBRARIES}
${GMX_EXTRA_LIBRARIES}
${GMX_TNG_LIBRARIES}
- ${EXTRAE_LIBRARIES}
${FFT_LIBRARIES} ${LINEAR_ALGEBRA_LIBRARIES}
${XML_LIBRARIES}
${THREAD_LIB} ${GMX_SHARED_LINKER_FLAGS})
#include "xdrf.h"
#include "xdr_datatype.h"
-#include "futil.h"
+#include "gromacs/utility/futil.h"
/* This is just for clarity - it can never be anything but 4! */
#define XDR_INT_SIZE 4
if (bSeekForwardOnly)
{
- low = gmx_ftell(fp);
+ low = gmx_ftell(fp)-header_size;
}
if (gmx_fseek(fp, 0, SEEK_END))
{
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
+#include "pdbio.h"
+
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <ctype.h>
+#include <stdlib.h>
#include <string.h>
-#include "sysstuff.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/types/ifunc.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/utility/futil.h"
+
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/topology/residuetypes.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/topology/topology.h"
#include "gromacs/utility/cstringutil.h"
-#include "vec.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
-#include "typedefs.h"
-#include "symtab.h"
-#include "pdbio.h"
-#include "vec.h"
-#include "copyrite.h"
-#include "futil.h"
-#include "atomprop.h"
-#include "physics.h"
-#include "pbc.h"
-#include "gmxfio.h"
-
-#include "gromacs/legacyheaders/gmx_fatal.h"
typedef struct {
int ai, aj;
};
- /* this is not very good,
- but these are only used in gmx_trjconv and gmx_editconv */
- static gmx_bool bWideFormat = FALSE;
#define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
- void set_pdb_wide_format(gmx_bool bSet)
- {
- bWideFormat = bSet;
- }
-
static void xlate_atomname_pdb2gmx(char *name)
{
int i, length;
gmx_conect conect, gmx_bool bTerSepChains)
{
gmx_conect_t *gc = (gmx_conect_t *)conect;
- char resnm[6], nm[6], pdbform[128], pukestring[100];
+ char resnm[6], nm[6], pukestring[100];
atom_id i, ii;
- int resind, resnr, type;
+ int resind, resnr;
+ enum PDB_record type;
unsigned char resic, ch;
+ char altloc;
real occup, bfac;
gmx_bool bOccup;
int nlongname = 0;
int chainnum, lastchainnum;
int lastresind, lastchainresind;
- gmx_residuetype_t rt;
+ gmx_residuetype_t*rt;
const char *p_restype;
const char *p_lastrestype;
bromacs(pukestring, 99);
fprintf(out, "TITLE %s\n", (title && title[0]) ? title : pukestring);
- if (bWideFormat)
- {
- fprintf(out, "REMARK This file does not adhere to the PDB standard\n");
- fprintf(out, "REMARK As a result of, some programs may not like it\n");
- }
if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) )
{
gmx_write_pdb_box(out, ePBC, box);
}
strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
+ resnm[sizeof(resnm)-1] = 0;
strncpy(nm, *atoms->atomname[i], sizeof(nm)-1);
+ nm[sizeof(nm)-1] = 0;
+
/* rename HG12 to 2HG1, etc. */
xlate_atomname_gmx2pdb(nm);
resnr = atoms->resinfo[resind].nr;
}
if (atoms->pdbinfo)
{
- type = atoms->pdbinfo[i].type;
+ type = (enum PDB_record)(atoms->pdbinfo[i].type);
+ altloc = atoms->pdbinfo[i].altloc;
+ if (!isalnum(altloc))
+ {
+ altloc = ' ';
+ }
occup = bOccup ? 1.0 : atoms->pdbinfo[i].occup;
bfac = atoms->pdbinfo[i].bfac;
}
else
{
- type = 0;
- occup = 1.0;
- bfac = 0.0;
- }
- if (bWideFormat)
- {
- strcpy(pdbform,
- "%-6s%5u %-4.4s %3.3s %c%4d%c %10.5f%10.5f%10.5f%8.4f%8.4f %2s\n");
- }
- else
- {
- /* Check whether atomname is an element name */
- if ((strlen(nm) < 4) && (gmx_strcasecmp(nm, atoms->atom[i].elem) != 0))
- {
- strcpy(pdbform, get_pdbformat());
- }
- else
- {
- strcpy(pdbform, get_pdbformat4());
- if (strlen(nm) > 4)
- {
- int maxwln = 20;
- if (nlongname < maxwln)
- {
- fprintf(stderr, "WARNING: Writing out atom name (%s) longer than 4 characters to .pdb file\n", nm);
- }
- else if (nlongname == maxwln)
- {
- fprintf(stderr, "WARNING: More than %d long atom names, will not write more warnings\n", maxwln);
- }
- nlongname++;
- }
- }
- strcat(pdbform, "%6.2f%6.2f %2s\n");
+ type = epdbATOM;
+ occup = 1.0;
+ bfac = 0.0;
+ altloc = ' ';
}
- fprintf(out, pdbform, pdbtp[type], (i+1)%100000, nm, resnm, ch, resnr,
- (resic == '\0') ? ' ' : resic,
- 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], occup, bfac, atoms->atom[i].elem);
+
+ gmx_fprintf_pdb_atomline(out,
+ type,
+ i+1,
+ nm,
+ altloc,
+ resnm,
+ ch,
+ resnr,
+ resic,
+ 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
+ occup,
+ bfac,
+ atoms->atom[i].elem);
+
if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
{
- fprintf(out, "ANISOU%5u %-4.4s%3.3s %c%4d%c %7d%7d%7d%7d%7d%7d\n",
+ fprintf(out, "ANISOU%5u %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n",
(i+1)%100000, nm, resnm, ch, resnr,
(resic == '\0') ? ' ' : resic,
atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
t_atom *atomn;
int j, k;
char nc = '\0';
- char anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12];
+ char anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12], elem[3];
char xc[12], yc[12], zc[12], occup[12], bfac[12];
unsigned char resic;
char chainid;
}
bfac[k] = nc;
+ /* 10 blanks */
+ j += 10;
+
+ /* Element name */
+ for (k = 0; (k < 2); k++, j++)
+ {
+ elem[k] = line[j];
+ }
+ elem[k] = nc;
+ trim(elem);
+
if (atoms->atom)
{
atomn = &(atoms->atom[natom]);
atomn->m = 0.0;
atomn->q = 0.0;
atomn->atomnumber = atomnumber;
- atomn->elem[0] = '\0';
+ strncpy(atomn->elem, elem, 4);
}
x[natom][XX] = strtod(xc, NULL)*0.1;
x[natom][YY] = strtod(yc, NULL)*0.1;
return gc;
}
- const char* get_pdbformat()
+ int
+ gmx_fprintf_pdb_atomline(FILE * fp,
+ enum PDB_record record,
+ int atom_seq_number,
+ const char * atom_name,
+ char alternate_location,
+ const char * res_name,
+ char chain_id,
+ int res_seq_number,
+ char res_insertion_code,
+ real x,
+ real y,
+ real z,
+ real occupancy,
+ real b_factor,
+ const char * element)
{
- static const char *pdbformat = "%-6s%5u %-4.4s%3.3s %c%4d%c %8.3f%8.3f%8.3f";
- return pdbformat;
- }
+ char tmp_atomname[6], tmp_resname[6];
+ gmx_bool start_name_in_col13;
+ int n;
- const char* get_pdbformat4()
- {
- static const char *pdbformat4 = "%-6s%5u %-4.4s %3.3s %c%4d%c %8.3f%8.3f%8.3f";
- return pdbformat4;
+ if (record != epdbATOM && record != epdbHETATM)
+ {
+ gmx_fatal(FARGS, "Can only print PDB atom lines as ATOM or HETATM records");
+ }
+
+ /* Format atom name */
+ if (atom_name != NULL)
+ {
+ /* If the atom name is an element name with two chars, it should start already in column 13.
+ * Otherwise it should start in column 14, unless the name length is 4 chars.
+ */
+ if ( (element != NULL) && (strlen(element) >= 2) && (gmx_strncasecmp(atom_name, element, 2) == 0) )
+ {
+ start_name_in_col13 = TRUE;
+ }
+ else
+ {
+ start_name_in_col13 = (strlen(atom_name) >= 4);
+ }
+ sprintf(tmp_atomname, start_name_in_col13 ? "" : " ");
+ strncat(tmp_atomname, atom_name, 4);
+ tmp_atomname[5] = '\0';
+ }
+ else
+ {
+ tmp_atomname[0] = '\0';
+ }
+
+ /* Format residue name */
+ strncpy(tmp_resname, (res_name != NULL) ? res_name : "", 4);
+ /* Make sure the string is terminated if strlen was > 4 */
+ tmp_resname[4] = '\0';
+ /* String is properly terminated, so now we can use strcat. By adding a
+ * space we can write it right-justified, and if the original name was
+ * three characters or less there will be a space added on the right side.
+ */
+ strcat(tmp_resname, " ");
+
+ /* Truncate integers so they fit */
+ atom_seq_number = atom_seq_number % 100000;
+ res_seq_number = res_seq_number % 10000;
+
+ n = fprintf(fp,
+ "%-6s%5d %-4.4s%c%4.4s%c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f %2s\n",
+ pdbtp[record],
+ atom_seq_number,
+ tmp_atomname,
+ alternate_location,
+ tmp_resname,
+ chain_id,
+ res_seq_number,
+ res_insertion_code,
+ x, y, z,
+ occupancy,
+ b_factor,
+ (element != NULL) ? element : "");
+
+ return n;
}
#ifndef GMX_FILEIO_PDBIO_H
#define GMX_FILEIO_PDBIO_H
-#include "../legacyheaders/sysstuff.h"
-#include "../legacyheaders/typedefs.h"
-#include "../legacyheaders/symtab.h"
-#include "../legacyheaders/atomprop.h"
+#include <stdio.h>
+
+#include "../legacyheaders/types/simple.h"
#ifdef __cplusplus
extern "C" {
#endif
+struct gmx_atomprop;
+struct t_atoms;
+struct t_topology;
+
typedef struct gmx_conect_t *gmx_conect;
- /* THE pdb format (for ATOM/HETATOM lines) */
- const char* get_pdbformat(void);
- const char* get_pdbformat4(void);
-
/* Enumerated type for pdb records. The other entries are ignored
* when reading a pdb file
*/
- enum {
+ enum PDB_record {
epdbATOM, epdbHETATM, epdbANISOU, epdbCRYST1, epdbCOMPND,
epdbMODEL, epdbENDMDL, epdbTER, epdbHEADER, epdbTITLE, epdbREMARK,
epdbCONECT, epdbNR
};
+ /* Write a PDB line with an ATOM or HETATM record directly to a file pointer.
+ *
+ * Returns the number of characters printed.
+ */
+ int
+ gmx_fprintf_pdb_atomline(FILE * fp,
+ enum PDB_record record,
+ int atom_seq_number,
+ const char * atom_name,
+ char alternate_location,
+ const char * res_name,
+ char chain_id,
+ int res_seq_number,
+ char res_insertion_code,
+ real x,
+ real y,
+ real z,
+ real occupancy,
+ real b_factor,
+ const char * element);
+
/* Enumerated value for indexing an uij entry (anisotropic temperature factors) */
enum {
U11, U22, U33, U12, U13, U23
};
- void set_pdb_wide_format(gmx_bool bSet);
- /* If bSet, use wider format for occupancy and bfactor */
-
void pdb_use_ter(gmx_bool bSet);
/* set read_pdbatoms to read upto 'TER' or 'ENDMDL' (default, bSet=FALSE).
This function is fundamentally broken as far as thread-safety is concerned.*/
* This function is fundamentally broken as far as thread-safety is concerned.
*/
-void write_pdbfile_indexed(FILE *out, const char *title, t_atoms *atoms,
+void write_pdbfile_indexed(FILE *out, const char *title, struct t_atoms *atoms,
rvec x[], int ePBC, matrix box, char chain,
int model_nr, atom_id nindex, const atom_id index[],
gmx_conect conect, gmx_bool bTerSepChains);
/* REALLY low level */
-void write_pdbfile(FILE *out, const char *title, t_atoms *atoms,
+void write_pdbfile(FILE *out, const char *title, struct t_atoms *atoms,
rvec x[], int ePBC, matrix box, char chain,
int model_nr, gmx_conect conect, gmx_bool bTerSepChains);
/* Low level pdb file writing routine.
* which may be useful for visualization purposes.
*/
-void get_pdb_atomnumber(t_atoms *atoms, gmx_atomprop_t aps);
+void get_pdb_atomnumber(struct t_atoms *atoms, struct gmx_atomprop *aps);
/* Routine to extract atomic numbers from the atom names */
int read_pdbfile(FILE *in, char *title, int *model_nr,
- t_atoms *atoms, rvec x[], int *ePBC, matrix box,
+ struct t_atoms *atoms, rvec x[], int *ePBC, matrix box,
gmx_bool bChange, gmx_conect conect);
/* Function returns number of atoms found.
* ePBC and gmx_conect structure may be NULL.
*/
void read_pdb_conf(const char *infile, char *title,
- t_atoms *atoms, rvec x[], int *ePBC, matrix box,
+ struct t_atoms *atoms, rvec x[], int *ePBC, matrix box,
gmx_bool bChange, gmx_conect conect);
/* Read a pdb file and extract ATOM and HETATM fields.
* Read a box from the CRYST1 line, return 0 box when no CRYST1 is found.
void gmx_conect_add(gmx_conect conect, int ai, int aj);
/* Add a connection between ai and aj (numbered from 0 to natom-1) */
-gmx_conect gmx_conect_generate(t_topology *top);
+gmx_conect gmx_conect_generate(struct t_topology *top);
/* Generate a conect structure from a topology */
gmx_conect gmx_conect_init(void);
#endif
#include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
-#include "gromacs/legacyheaders/main.h"
-#include "gromacs/legacyheaders/physics.h"
+#include "gromacs/legacyheaders/types/ifunc.h"
+
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/basenetwork.h"
#include "gromacs/utility/common.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/programcontext.h"
-#include "gmxfio.h"
static const char *modeToVerb(char mode)
{
gmx_int64_t numMolecules,
tng_molecule_t *tngMol)
{
+ tng_chain_t tngChain = NULL;
+ tng_residue_t tngRes = NULL;
+
if (tng_molecule_add(tng, moleculeName, tngMol) != TNG_SUCCESS)
{
gmx_file("Cannot add molecule to TNG molecular system.");
{
const t_resinfo *resInfo = &atoms->resinfo[at->resind];
char chainName[2] = {resInfo->chainid, 0};
- tng_chain_t tngChain = NULL;
- tng_residue_t tngRes = NULL;
tng_atom_t tngAtom = NULL;
+ t_atom *prevAtom;
- if (tng_molecule_chain_find (tng, *tngMol, chainName,
- (gmx_int64_t)-1, &tngChain) !=
- TNG_SUCCESS)
+ if (atomIndex > 0)
+ {
+ prevAtom = &atoms->atom[atomIndex - 1];
+ }
+ else
{
- tng_molecule_chain_add (tng, *tngMol, chainName,
- &tngChain);
+ prevAtom = 0;
}
- /* FIXME: When TNG supports both residue index and residue
- * number the latter should be used. Wait for TNG 2.0*/
- if (tng_chain_residue_find(tng, tngChain, *resInfo->name,
- at->resind + 1, &tngRes)
- != TNG_SUCCESS)
+ /* If this is the first atom or if the residue changed add the
+ * residue to the TNG molecular system. */
+ if (!prevAtom || resInfo != &atoms->resinfo[prevAtom->resind])
{
+ /* If this is the first atom or if the chain changed add
+ * the chain to the TNG molecular system. */
+ if (!prevAtom || resInfo->chainid !=
+ atoms->resinfo[prevAtom->resind].chainid)
+ {
+ tng_molecule_chain_add(tng, *tngMol, chainName,
+ &tngChain);
+ }
+ /* FIXME: When TNG supports both residue index and residue
+ * number the latter should be used. Wait for TNG 2.0*/
tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
}
tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
/* This file is completely threadsafe - keep it that way! */
+#include <stdlib.h>
#include <string.h>
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gmx_fatal.h"
#include "macros.h"
#include "names.h"
-#include "symtab.h"
-#include "futil.h"
+#include "gromacs/utility/futil.h"
#include "filenm.h"
#include "gmxfio.h"
#include "tpxio.h"
#include "txtdump.h"
#include "confio.h"
-#include "atomprop.h"
#include "copyrite.h"
-#include "vec.h"
-#include "mtop_util.h"
+
+#include "gromacs/math/vec.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
#define TPX_TAG_RELEASE "release"
tpxv_ComputationalElectrophysiology = 96, /**< support for ion/water position swaps (computational electrophysiology) */
tpxv_Use64BitRandomSeed, /**< change ld_seed from int to gmx_int64_t */
tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */
- tpxv_InteractiveMolecularDynamics /**< interactive molecular dynamics (IMD) */
+ tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */
+ tpxv_RemoveObsoleteParameters1 /**< remove optimize_fft, dihre_fc, nstcheckpoint */
};
/*! \brief Version number of the file format written to run input
*
* When developing a feature branch that needs to change the run input
* file format, change tpx_tag instead. */
- static const int tpx_version = tpxv_InteractiveMolecularDynamics;
+ static const int tpx_version = tpxv_RemoveObsoleteParameters1;
/* This number should only be increased when you edit the TOPOLOGY section
int i, j, k, *tmp, idum = 0;
real rdum, bd_temp;
rvec vdum;
- gmx_bool bSimAnn;
+ gmx_bool bSimAnn, bdum = 0;
real zerotemptime, finish_t, init_temp, finish_temp;
if (file_version != tpx_version)
gmx_fio_do_int(fio, idum);
ir->nsteps = idum;
}
+
if (file_version > 25)
{
if (file_version >= 62)
}
gmx_fio_do_int(fio, ir->ns_type);
gmx_fio_do_int(fio, ir->nstlist);
- gmx_fio_do_int(fio, ir->ndelta);
+ gmx_fio_do_int(fio, idum); /* used to be ndelta; not used anymore */
if (file_version < 41)
{
gmx_fio_do_int(fio, idum);
}
ir->nstcomm = abs(ir->nstcomm);
- if (file_version > 25)
- {
- gmx_fio_do_int(fio, ir->nstcheckpoint);
- }
- else
+ /* ignore nstcheckpoint */
+ if (file_version > 25 && file_version < tpxv_RemoveObsoleteParameters1)
{
- ir->nstcheckpoint = 0;
+ gmx_fio_do_int(fio, idum);
}
gmx_fio_do_int(fio, ir->nstcgsteep);
gmx_fio_do_real(fio, ir->epsilon_surface);
}
- gmx_fio_do_gmx_bool(fio, ir->bOptFFT);
+ /* ignore bOptFFT */
+ if (file_version < tpxv_RemoveObsoleteParameters1)
+ {
+ gmx_fio_do_gmx_bool(fio, bdum);
+ }
if (file_version >= 93)
{
ir->orires_tau = 0;
ir->nstorireout = 0;
}
+
+ /* ignore dihre_fc */
if (file_version >= 26 && file_version < 79)
{
- gmx_fio_do_real(fio, ir->dihre_fc);
+ gmx_fio_do_real(fio, rdum);
if (file_version < 56)
{
gmx_fio_do_real(fio, rdum);
gmx_fio_do_int(fio, idum);
}
}
- else
- {
- ir->dihre_fc = 0;
- }
gmx_fio_do_real(fio, ir->em_stepsize);
gmx_fio_do_real(fio, ir->em_tol);
gmx_fio_ndo_int(fio, block->a, block->nra);
}
+ /* This is a primitive routine to make it possible to translate atomic numbers
+ * to element names when reading TPR files, without making the Gromacs library
+ * directory a dependency on mdrun (which is the case if we need elements.dat).
+ */
+ static const char *
+ atomicnumber_to_element(int atomicnumber)
+ {
+ const char * p;
+
+ /* This does not have to be complete, so we only include elements likely
+ * to occur in PDB files.
+ */
+ switch (atomicnumber)
+ {
+ case 1: p = "H"; break;
+ case 5: p = "B"; break;
+ case 6: p = "C"; break;
+ case 7: p = "N"; break;
+ case 8: p = "O"; break;
+ case 9: p = "F"; break;
+ case 11: p = "Na"; break;
+ case 12: p = "Mg"; break;
+ case 15: p = "P"; break;
+ case 16: p = "S"; break;
+ case 17: p = "Cl"; break;
+ case 18: p = "Ar"; break;
+ case 19: p = "K"; break;
+ case 20: p = "Ca"; break;
+ case 25: p = "Mn"; break;
+ case 26: p = "Fe"; break;
+ case 28: p = "Ni"; break;
+ case 29: p = "Cu"; break;
+ case 30: p = "Zn"; break;
+ case 35: p = "Br"; break;
+ case 47: p = "Ag"; break;
+ default: p = ""; break;
+ }
+ return p;
+ }
+
+
static void do_atom(t_fileio *fio, t_atom *atom, int ngrp, gmx_bool bRead,
int file_version, gmx_groups_t *groups, int atnr)
{
- int i, myngrp;
+ int i, myngrp;
+ char * p_elem;
gmx_fio_do_real(fio, atom->m);
gmx_fio_do_real(fio, atom->q);
if (file_version >= 52)
{
gmx_fio_do_int(fio, atom->atomnumber);
+ if (bRead)
+ {
+ /* Set element string from atomic number if present.
+ * This routine returns an empty string if the name is not found.
+ */
+ strncpy(atom->elem, atomicnumber_to_element(atom->atomnumber), 4);
+ /* avoid warnings about potentially unterminated string */
+ atom->elem[3] = '\0';
+ }
}
else if (bRead)
{
static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead,
int file_version)
{
- int mt, mb, i;
- t_blocka dumb;
+ int mt, mb, i;
+ t_blocka dumb;
if (bRead)
{
#ifndef GMX_FILEIO_TRX_H
#define GMX_FILEIO_TRX_H
-#include "../legacyheaders/types/atoms.h"
+#include "../math/vectypes.h"
+#include "../utility/basedefinitions.h"
+#include "../utility/real.h"
#ifdef __cplusplus
extern "C" {
#endif
+struct t_atoms;
+
typedef struct gmxvmdplugin t_gmxvmdplugin;
-typedef struct trxframe
+typedef struct t_trxframe
{
int flags; /* flags for read_first/next_frame */
int not_ok; /* integrity flags */
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 */
real lambda; /* free energy perturbation lambda */
int fep_state; /* which fep state are we in? */
gmx_bool bAtoms;
- t_atoms *atoms; /* atoms struct (natoms) */
+ struct t_atoms *atoms; /* atoms struct (natoms) */
gmx_bool bPrec;
real prec; /* precision of x, fraction of 1 nm */
gmx_bool bX;
#include <assert.h>
#include <math.h>
-#include "sysstuff.h"
-#include "typedefs.h"
#ifdef GMX_USE_PLUGINS
#include "vmdio.h"
#endif
-#include "gromacs/utility/smalloc.h"
-#include "pbc.h"
#include "gmxfio.h"
#include "trxio.h"
#include "tpxio.h"
#include "tngio.h"
#include "tngio_for_tools.h"
#include "names.h"
-#include "vec.h"
-#include "futil.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
#include "xtcio.h"
#include "pdbio.h"
#include "confio.h"
#include "xdrf.h"
#include "gromacs/fileio/timecontrol.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/fileio/trx.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
/* defines for frame counter output */
#define SKIP1 10
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))
#include <math.h>
#include <stdio.h>
#include <string.h>
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
#include "txtdump.h"
#include "bondf.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
#include "typedefs.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "gstat.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/trxio.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/fatalerror.h"
void print_one(const output_env_t oenv, const char *base, const char *name,
const char *title, const char *ylabel, int nf, real time[],
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);
}
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+
#include <math.h>
+#include <stdlib.h>
#include <string.h>
+
#include "gromacs/commandline/pargs.h"
-#include "sysstuff.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "gmx_fatal.h"
-#include "vec.h"
-#include "pbc.h"
-#include "gromacs/fileio/futil.h"
-#include "index.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/topology/index.h"
#include "gromacs/fileio/pdbio.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/matio.h"
-#include "mshift.h"
-#include "xvgr.h"
-#include "rmpbc.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
+#include "gromacs/pbcutil/rmpbc.h"
#include "txtdump.h"
#include "eigio.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "gmx_ana.h"
#include "gromacs/math/do_fit.h"
{
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]);
}
rvec *x;
real *b = NULL;
matrix box;
- char *resnm, *atnm, pdbform[STRLEN];
+ char *resnm, *atnm;
gmx_bool bPDB, b4D;
FILE *out;
}
if ( ( b4D || bSplit ) && bPDB)
{
- strcpy(pdbform, get_pdbformat());
- strcat(pdbform, "%8.4f%8.4f\n");
-
out = gmx_ffopen(threedplotfile, "w");
fprintf(out, "HEADER %s\n", str);
if (b4D)
fprintf(out, "TER\n");
j = 0;
}
- fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1,
- 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i]);
+ gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
+ 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
if (j > 0)
{
fprintf(out, "CONECT%5d%5d\n", i, i+1);
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+
#include <math.h>
+#include <stdlib.h>
#include <string.h>
+
#include "gromacs/commandline/pargs.h"
-#include "sysstuff.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "gmx_fatal.h"
-#include "vec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/math/vec.h"
#include "copyrite.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "readinp.h"
#include "txtdump.h"
#include "gstat.h"
#include "gromacs/statistics/statistics.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
#include "gmx_ana.h"
#include "geminate.h"
}
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);
#include <math.h>
#include <string.h>
-#include "sysstuff.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/commandline/pargs.h"
#include "copyrite.h"
-#include "vec.h"
-#include "index.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/topology/index.h"
#include "macros.h"
-#include "gmx_fatal.h"
-#include "xvgr.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
#include "gstat.h"
#include "gromacs/fileio/trnio.h"
#include "gmx_ana.h"
{
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++)
{
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include <stdio.h>
+
#include <math.h>
+#include <stdio.h>
+#include <string.h>
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/pdbio.h"
-#include "gmx_fatal.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
#include "gstat.h"
#include "macros.h"
#include "gromacs/math/utilities.h"
-#include "physics.h"
-#include "index.h"
+#include "gromacs/math/units.h"
+#include "gromacs/topology/residuetypes.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/tpxio.h"
-#include <string.h>
-#include "sysstuff.h"
#include "txtdump.h"
#include "typedefs.h"
-#include "vec.h"
-#include "xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
#include "gromacs/fileio/matio.h"
#include "gmx_ana.h"
return j;
}
-static void histogramming(FILE *log, int nbin, gmx_residuetype_t rt,
+static void histogramming(FILE *log, int nbin, gmx_residuetype_t *rt,
int nf, int maxchi, real **dih,
int nlist, t_dlist dlist[],
atom_id index[],
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;
}
{
FILE *fp;
int nh[edMax];
- char buf[STRLEN];
int i, Dih, Xi;
real S2Max, S2Min;
x0 *= 10.0; /* nm -> angstrom */
y0 *= 10.0; /* nm -> angstrom */
z0 *= 10.0; /* nm -> angstrom */
- sprintf(buf, "%s%%6.f%%6.2f\n", get_pdbformat());
for (i = 0; (i < 10); i++)
{
- fprintf(fp, buf, "ATOM ", atoms->nr+1+i, "CA", "LEG", ' ',
- atoms->nres+1, ' ', x0, y0, z0+(1.2*i), 0.0, -0.1*i);
+ gmx_fprintf_pdb_atomline(fp, epdbATOM, atoms->nr+1+i, "CA", ' ', "LEG", ' ', atoms->nres+1, ' ',
+ x0, y0, z0+(1.2*i), 0.0, -0.1*i, "");
}
gmx_ffclose(fp);
}
gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc;
real dt = 0, traj_t_ns;
output_env_t oenv;
- gmx_residuetype_t rt;
+ gmx_residuetype_t *rt;
atom_id isize, *index;
int ndih, nactdih, nf;
#endif
#include <math.h>
+#include <stdlib.h>
#include <string.h>
#include "macros.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/utility/cstringutil.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "macros.h"
-#include "index.h"
+#include "gromacs/topology/index.h"
#include "gromacs/random/random.h"
-#include "pbc.h"
-#include "rmpbc.h"
-#include "xvgr.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/fileio/matio.h"
#include "cmat.h"
#include "gromacs/fileio/trnio.h"
#include "gromacs/linearalgebra/eigensolver.h"
#include "gromacs/math/do_fit.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
/* print to two file pointers at once (i.e. stderr and log) */
static gmx_inline
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",
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include <math.h>
-#include <ctype.h>
-#include "sysstuff.h"
+#include <ctype.h>
+#include <math.h>
+#include <stdlib.h>
#include <string.h>
+
#include "gromacs/utility/cstringutil.h"
#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
#include "macros.h"
#include "gstat.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "pbc.h"
-#include "gromacs/fileio/futil.h"
-#include "gromacs/commandline/pargs.h"
-#include "index.h"
+#include "viewit.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/index.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "gmx_ana.h"
#include "macros.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
typedef struct {
char *atomname;
return nr;
}
- void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis)
+ void center_coords(t_atoms *atoms, atom_id *index_center, int ncenter,
+ matrix box, rvec x0[])
{
- int i, m;
+ int i, k, m;
real tmass, mm;
rvec com, shift, box_center;
tmass = 0;
clear_rvec(com);
- for (i = 0; (i < atoms->nr); i++)
+ for (k = 0; (k < ncenter); k++)
{
+ i = index_center[k];
+ if (i >= atoms->nr)
+ {
+ gmx_fatal(FARGS, "Index %d refers to atom %d, which is larger than natoms (%d).",
+ k+1, i+1, atoms->nr);
+ }
mm = atoms->atom[i].m;
tmass += mm;
for (m = 0; (m < DIM); m++)
}
calc_box_center(ecenterDEF, box, box_center);
rvec_sub(box_center, com, shift);
- shift[axis] -= box_center[axis];
+ /* Important - while the center was calculated based on a group, we should move all atoms */
for (i = 0; (i < atoms->nr); i++)
{
rvec_dec(x0[i], shift);
int ePBC,
int axis, int nr_grps, real *slWidth,
t_electron eltab[], int nr, gmx_bool bCenter,
- const output_env_t oenv)
+ atom_id *index_center, int ncenter,
+ gmx_bool bRelative, const output_env_t oenv)
{
rvec *x0; /* coordinates without pbc */
matrix box; /* box (3x3) */
slice; /* current slice */
t_electron *found; /* found by bsearch */
t_electron sought; /* thingie thought by bsearch */
+ real boxSz, aveBox;
gmx_rmpbc_t gpbc = NULL;
real t,
gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
}
+ aveBox = 0;
+
if (!*nslices)
{
*nslices = (int)(box[axis][axis] * 10); /* default value */
+ fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
}
- fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
snew(*slDensity, nr_grps);
for (i = 0; i < nr_grps; i++)
{
gmx_rmpbc(gpbc, natoms, box, x0);
+ /* Translate atoms so the com of the center-group is in the
+ * box geometrical center.
+ */
if (bCenter)
{
- center_coords(&top->atoms, box, x0, axis);
+ center_coords(&top->atoms, index_center, ncenter, box, x0);
}
- *slWidth = box[axis][axis]/(*nslices);
invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
+ if (bRelative)
+ {
+ *slWidth = 1.0/(*nslices);
+ boxSz = 1.0;
+ }
+ else
+ {
+ *slWidth = box[axis][axis]/(*nslices);
+ boxSz = box[axis][axis];
+ }
+
+ aveBox += box[axis][axis];
+
for (n = 0; n < nr_grps; n++)
{
for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
z -= box[axis][axis];
}
+ if (bRelative)
+ {
+ z = z/box[axis][axis];
+ }
+
/* determine which slice atom is in */
- slice = (z / (*slWidth));
+ if (bCenter)
+ {
+ slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
+ }
+ else
+ {
+ slice = (z / (*slWidth));
+ }
sought.nr_el = 0;
sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
nr_frames);
+ if (bRelative)
+ {
+ aveBox /= nr_frames;
+ *slWidth = aveBox/(*nslices);
+ }
+
for (n = 0; n < nr_grps; n++)
{
for (i = 0; i < *nslices; i++)
void calc_density(const char *fn, atom_id **index, int gnx[],
double ***slDensity, int *nslices, t_topology *top, int ePBC,
int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
- const output_env_t oenv)
+ atom_id *index_center, int ncenter,
+ gmx_bool bRelative, const output_env_t oenv)
{
- rvec *x0; /* coordinates without pbc */
- matrix box; /* box (3x3) */
+ rvec *x0; /* coordinates without pbc */
+ matrix box; /* box (3x3) */
double invvol;
- int natoms; /* nr. atoms in trj */
+ int natoms; /* nr. atoms in trj */
t_trxstatus *status;
- int **slCount, /* nr. of atoms in one slice for a group */
- i, j, n, /* loop indices */
- teller = 0,
+ int **slCount, /* nr. of atoms in one slice for a group */
+ i, j, n, /* loop indices */
ax1 = 0, ax2 = 0,
nr_frames = 0, /* number of frames */
slice; /* current slice */
real t,
z;
+ real boxSz, aveBox;
char *buf; /* for tmp. keeping atomname */
gmx_rmpbc_t gpbc = NULL;
gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
}
+ aveBox = 0;
+
if (!*nslices)
{
*nslices = (int)(box[axis][axis] * 10); /* default value */
{
gmx_rmpbc(gpbc, natoms, box, x0);
+ /* Translate atoms so the com of the center-group is in the
+ * box geometrical center.
+ */
if (bCenter)
{
- center_coords(&top->atoms, box, x0, axis);
+ center_coords(&top->atoms, index_center, ncenter, box, x0);
}
- *slWidth = box[axis][axis]/(*nslices);
invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
- teller++;
+
+ if (bRelative)
+ {
+ *slWidth = 1.0/(*nslices);
+ boxSz = 1.0;
+ }
+ else
+ {
+ *slWidth = box[axis][axis]/(*nslices);
+ boxSz = box[axis][axis];
+ }
+
+ aveBox += box[axis][axis];
for (n = 0; n < nr_grps; n++)
{
z -= box[axis][axis];
}
+ if (bRelative)
+ {
+ z = z/box[axis][axis];
+ }
+
/* determine which slice atom is in */
- slice = (int)(z / (*slWidth));
+ if (bCenter)
+ {
+ slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
+ }
+ else
+ {
+ slice = floor(z / (*slWidth));
+ }
+
+ /* Slice should already be 0<=slice<nslices, but we just make
+ * sure we are not hit by IEEE rounding errors since we do
+ * math operations after applying PBC above.
+ */
+ if (slice < 0)
+ {
+ slice += *nslices;
+ }
+ else if (slice >= *nslices)
+ {
+ slice -= *nslices;
+ }
+
(*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
}
}
fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
nr_frames);
+ if (bRelative)
+ {
+ aveBox /= nr_frames;
+ *slWidth = aveBox/(*nslices);
+ }
+
for (n = 0; n < nr_grps; n++)
{
for (i = 0; i < *nslices; i++)
void plot_density(double *slDensity[], const char *afile, int nslices,
int nr_grps, char *grpname[], real slWidth,
const char **dens_opt,
- gmx_bool bSymmetrize, const output_env_t oenv)
+ gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
+ const output_env_t oenv)
{
FILE *den;
+ const char *title = NULL;
+ const char *xlabel = NULL;
const char *ylabel = NULL;
int slice, n;
real ddd;
+ real axispos;
+
+ title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
+
+ if (bCenter)
+ {
+ xlabel = bRelative ?
+ "Average relative position from center (nm)" :
+ "Relative position from center (nm)";
+ }
+ else
+ {
+ xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
+ }
switch (dens_opt[0][0])
{
case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
}
- den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
+ den = xvgropen(afile,
+ title, xlabel, ylabel, oenv);
xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
for (slice = 0; (slice < nslices); slice++)
{
- fprintf(den, "%12g ", slice * slWidth);
+ if (bCenter)
+ {
+ axispos = (slice - nslices/2.0 + 0.5) * slWidth;
+ }
+ else
+ {
+ axispos = (slice + 0.5) * slWidth;
+ }
+ fprintf(den, "%12g ", axispos);
for (n = 0; (n < nr_grps); n++)
{
if (bSymmetrize)
"[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
"For the total density of NPT simulations, use [gmx-energy] instead.",
"[PAR]",
+
+ "Option [TT]-center[tt] performs the histogram binning relative to the center",
+ "of an arbitrary group, in absolute box coordinates. If you are calculating",
+ "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
+ "bZ/2 if you center based on the entire system.",
+ "Note that this behaviour has changed in Gromacs 5.0; earlier versions",
+ "merely performed a static binning in (0,bZ) and shifted the output. Now",
+ "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
+
+ "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
+ "automatically turn on [TT]-center[tt] too.",
+
+ "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
+ "box coordinates, and scales the final output with the average box dimension",
+ "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
+
"Densities are in kg/m^3, and number densities or electron densities can also be",
"calculated. For electron densities, a file describing the number of",
"electrons for each type of atom should be provided using [TT]-ei[tt].",
"The first line contains the number of lines to read from the file.",
"There should be one line for each unique atom name in your system.",
"The number of electrons for each atom is modified by its atomic",
- "partial charge."
+ "partial charge.[PAR]",
+
+ "IMPORTANT CONSIDERATIONS FOR BILAYERS[PAR]",
+ "One of the most common usage scenarios is to calculate the density of various",
+ "groups across a lipid bilayer, typically with the z axis being the normal",
+ "direction. For short simulations, small systems, and fixed box sizes this",
+ "will work fine, but for the more general case lipid bilayers can be complicated.",
+ "The first problem that while both proteins and lipids have low volume",
+ "compressibility, lipids have quite high area compressiblity. This means the",
+ "shape of the box (thickness and area/lipid) will fluctuate substantially even",
+ "for a fully relaxed system. Since Gromacs places the box between the origin",
+ "and positive coordinates, this in turn means that a bilayer centered in the",
+ "box will move a bit up/down due to these fluctuations, and smear out your",
+ "profile. The easiest way to fix this (if you want pressure coupling) is",
+ "to use the [TT]-center[tt] option that calculates the density profile with",
+ "respect to the center of the box. Note that you can still center on the",
+ "bilayer part even if you have a complex non-symmetric system with a bilayer",
+ "and, say, membrane proteins - then our output will simply have more values",
+ "on one side of the (center) origin reference.[PAR]",
+
+ "Even the centered calculation will lead to some smearing out the output",
+ "profiles, as lipids themselves are compressed and expanded. In most cases",
+ "you probably want this (since it corresponds to macroscopic experiments),",
+ "but if you want to look at molecular details you can use the [TT]-relative[tt]",
+ "option to attempt to remove even more of the effects of volume fluctuations.[PAR]",
+
+ "Finally, large bilayers that are not subject to a surface tension will exhibit",
+ "undulatory fluctuations, where there are 'waves' forming in the system.",
+ "This is a fundamental property of the biological system, and if you are",
+ "comparing against experiments you likely want to include the undulation",
+ "smearing effect.",
+ "",
};
output_env_t oenv;
static int ngrps = 1; /* nr. of groups */
static gmx_bool bSymmetrize = FALSE;
static gmx_bool bCenter = FALSE;
+ static gmx_bool bRelative = FALSE;
+
t_pargs pa[] = {
{ "-d", FALSE, etSTR, {&axtitle},
"Take the normal on the membrane in direction X, Y or Z." },
"Density"},
{ "-ng", FALSE, etINT, {&ngrps},
"Number of groups of which to compute densities." },
- { "-symm", FALSE, etBOOL, {&bSymmetrize},
+ { "-center", FALSE, etBOOL, {&bCenter},
+ "Perform the binning relative to the center of the (changing) box. Useful for bilayers." },
+ { "-symm", FALSE, etBOOL, {&bSymmetrize},
"Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
- { "-center", FALSE, etBOOL, {&bCenter},
- "Shift the center of mass along the axis to zero. This means if your axis is Z and your box is bX, bY, bZ, the center of mass will be at bX/2, bY/2, 0."}
+ { "-relative", FALSE, etBOOL, {&bRelative},
+ "Use relative coordinates for changing boxes and scale output by average dimensions." }
};
const char *bugs[] = {
"When calculating electron densities, atomnames are used instead of types. This is bad.",
};
- double **density; /* density per slice */
- real slWidth; /* width of one slice */
- char **grpname; /* groupnames */
- int nr_electrons; /* nr. electrons */
- int *ngx; /* sizes of groups */
- t_electron *el_tab; /* tabel with nr. of electrons*/
- t_topology *top; /* topology */
+ double **density; /* density per slice */
+ real slWidth; /* width of one slice */
+ char *grpname_center; /* centering group name */
+ char **grpname; /* groupnames */
+ int nr_electrons; /* nr. electrons */
+ int ncenter; /* size of centering group */
+ int *ngx; /* sizes of groups */
+ t_electron *el_tab; /* tabel with nr. of electrons*/
+ t_topology *top; /* topology */
int ePBC;
- atom_id **index; /* indices for all groups */
+ atom_id *index_center; /* index for centering group */
+ atom_id **index; /* indices for all groups */
int i;
t_filenm fnm[] = { /* files for g_density */
snew(index, ngrps);
snew(ngx, ngrps);
+ if (bCenter)
+ {
+ fprintf(stderr,
+ "\nNote: that the center of mass is calculated inside the box without applying\n"
+ "any special periodicity. If necessary, it is your responsibility to first use\n"
+ "trjconv to make sure atoms in this group are placed in the right periodicity.\n\n"
+ "Select the group to center density profiles around:\n");
+ get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ncenter,
+ &index_center, &grpname_center);
+ }
+ else
+ {
+ ncenter = 0;
+ }
+
+ fprintf(stderr, "\nSelect %d group%s to calculate density for:\n", ngrps, (ngrps > 1) ? "s" : "");
get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
if (dens_opt[0][0] == 'e')
calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
&nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
- nr_electrons, bCenter, oenv);
+ nr_electrons, bCenter, index_center, ncenter,
+ bRelative, oenv);
}
else
{
calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
- ePBC, axis, ngrps, &slWidth, bCenter, oenv);
+ ePBC, axis, ngrps, &slWidth, bCenter, index_center, ncenter,
+ bRelative, oenv);
}
plot_density(density, opt2fn("-o", NFILE, fnm),
nslices, ngrps, grpname, slWidth, dens_opt,
- bSymmetrize, oenv);
+ bCenter, bRelative, bSymmetrize, oenv);
do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */
return 0;
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include <stdio.h>
+
#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
#include "gromacs/fileio/confio.h"
#include "copyrite.h"
-#include "gmx_fatal.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
#include "gstat.h"
#include "macros.h"
#include "gromacs/math/utilities.h"
-#include "physics.h"
-#include "index.h"
+#include "gromacs/math/units.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/commandline/pargs.h"
-#include <string.h>
-#include "sysstuff.h"
#include "txtdump.h"
#include "typedefs.h"
-#include "vec.h"
-#include "xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
#include "correl.h"
#include "gmx_ana.h"
#include "gromacs/fft/fft.h"
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);
#include "gromacs/fileio/pdbio.h"
#include "gromacs/fileio/confio.h"
-#include "symtab.h"
-#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/strdb.h"
-#include "index.h"
-#include "vec.h"
+#include "gromacs/topology/index.h"
#include "typedefs.h"
#include "gromacs/gmxlib/conformation-utilities.h"
-#include "physics.h"
-#include "atomprop.h"
+#include "gromacs/math/units.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "pbc.h"
#include "princ.h"
#include "txtdump.h"
#include "viewit.h"
-#include "rmpbc.h"
#include "gmx_ana.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
typedef struct
{
for (i = 0; i < nat; i++)
{
- fprintf(out, get_pdbformat(), "ATOM", a0 + i, "C", "BOX", 'K' + i
- / NCUCVERT, r0 + i, 10 * vert[i][XX], 10 * vert[i][YY], 10
- * vert[i][ZZ]);
- fprintf(out, "\n");
+ gmx_fprintf_pdb_atomline(out, epdbATOM, a0 + i, "C", ' ', "BOX", 'K' + i / NCUCVERT, r0 + i, ' ',
+ 10*vert[i][XX], 10*vert[i][YY], 10*vert[i][ZZ], 1.0, 0.0, "");
}
edge = compact_unitcell_edges();
{
for (x = 0; x <= 1; x++)
{
- fprintf(out, get_pdbformat(), "ATOM", a0 + i, "C", "BOX", 'K' + i
- / 8, r0 + i, x * 10 * box[XX][XX],
- y * 10 * box[YY][YY], z * 10 * box[ZZ][ZZ]);
- fprintf(out, "\n");
+ gmx_fprintf_pdb_atomline(out, epdbATOM, a0 + i, "C", ' ', "BOX", 'K' + i/8, r0+i, ' ',
+ x * 10 * box[XX][XX], y * 10 * box[YY][YY], z * 10 * box[ZZ][ZZ], 1.0, 0.0, "");
i++;
}
}
out = gmx_ffopen(outfile, "w");
if (bMead)
{
- set_pdb_wide_format(TRUE);
fprintf(out, "REMARK "
"The B-factors in this file hold atomic radii\n");
fprintf(out, "REMARK "
#include "gromacs/utility/cstringutil.h"
#include "typedefs.h"
-#include "gmx_fatal.h"
-#include "vec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/math/vec.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/commandline/pargs.h"
#include "names.h"
#include "macros.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
#include "gstat.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "gromacs/fileio/matio.h"
#include "gromacs/fileio/strdb.h"
#include "gmx_ana.h"
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]);
#endif
#include <math.h>
+#include <stdlib.h>
#include <string.h>
#include "typedefs.h"
-#include "gmx_fatal.h"
-#include "vec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/math/vec.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/fileio/enxio.h"
#include "names.h"
#include "copyrite.h"
#include "macros.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
#include "gstat.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "viewit.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "gmx_ana.h"
#include "mdebin.h"
{
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);
}
#endif
#include <math.h>
-#include "sysstuff.h"
-#include "gromacs/commandline/pargs.h"
#include <string.h>
+
+#include "gromacs/commandline/pargs.h"
#include "gromacs/utility/smalloc.h"
#include "typedefs.h"
#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "macros.h"
-#include "vec.h"
-#include "index.h"
-#include "gmx_fatal.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
#include "gmx_ana.h"
int gmx_genpr(int argc, char *argv[])
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);
}
}
}
#include <string.h>
#include "macros.h"
-#include "vec.h"
-#include "sysstuff.h"
+#include "gromacs/math/vec.h"
#include "typedefs.h"
#include "gromacs/fileio/filenm.h"
#include "gromacs/commandline/pargs.h"
-#include "gromacs/fileio/futil.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/fileio/matio.h"
-#include "xvgr.h"
-#include "index.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/topology/index.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "rmpbc.h"
-#include "pbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gmx_ana.h"
int **nmat, **totnmat;
real *mean_n;
int *tot_n;
- matrix box;
+ matrix box = {{0}};
output_env_t oenv;
gmx_rmpbc_t gpbc = NULL;
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 <math.h>
#include <stdlib.h>
-
-#include "sysstuff.h"
#include <string.h>
+
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "pbc.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/commandline/pargs.h"
-#include "index.h"
+#include "gromacs/topology/index.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "rmpbc.h"
+#include "gromacs/pbcutil/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));
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+
#include <math.h>
#include <string.h>
-#include "sysstuff.h"
-#include "physics.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/commandline/pargs.h"
-#include "vec.h"
-#include "index.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/topology/index.h"
#include "macros.h"
-#include "xvgr.h"
-#include "rmpbc.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
+#include "gromacs/pbcutil/rmpbc.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "gmx_ana.h"
/* 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. */
#include <string.h>
#include "gromacs/commandline/pargs.h"
-#include "sysstuff.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "vec.h"
-#include "pbc.h"
-#include "gromacs/fileio/futil.h"
-#include "index.h"
-#include "mshift.h"
-#include "xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/topology/index.h"
#include "princ.h"
-#include "rmpbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
#include "txtdump.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
++#include "gromacs/fileio/xvgr.h"
#include "gstat.h"
#include "gmx_ana.h"
matrix axes, box;
output_env_t oenv;
gmx_rmpbc_t gpbc = NULL;
+ char ** legend;
-
- t_filenm fnm[] = {
+ t_filenm fnm[] = {
{ efTRX, "-f", NULL, ffREAD },
{ efTPS, NULL, NULL, ffREAD },
{ efNDX, NULL, NULL, ffOPTRD },
- { efDAT, "-a1", "paxis1", ffWRITE },
- { efDAT, "-a2", "paxis2", ffWRITE },
- { efDAT, "-a3", "paxis3", ffWRITE },
- { efDAT, "-om", "moi", ffWRITE }
+ { efXVG, "-a1", "paxis1", ffWRITE },
+ { efXVG, "-a2", "paxis2", ffWRITE },
+ { efXVG, "-a3", "paxis3", ffWRITE },
+ { efXVG, "-om", "moi", ffWRITE }
};
#define NFILE asize(fnm)
return 0;
}
- axis1 = gmx_ffopen(opt2fn("-a1", NFILE, fnm), "w");
- axis2 = gmx_ffopen(opt2fn("-a2", NFILE, fnm), "w");
- axis3 = gmx_ffopen(opt2fn("-a3", NFILE, fnm), "w");
- fmoi = gmx_ffopen(opt2fn("-om", NFILE, fnm), "w");
+ snew(legend, DIM);
+ for (i = 0; i < DIM; i++)
+ {
+ snew(legend[i], STRLEN);
+ sprintf(legend[i], "%c component", 'X'+i);
+ }
+
+ axis1 = xvgropen(opt2fn("-a1", NFILE, fnm), "Principal axis 1 (major axis)",
+ output_env_get_xvgr_tlabel(oenv), "Component (nm)", oenv);
+ xvgr_legend(axis1, DIM, (const char **)legend, oenv);
+
+ axis2 = xvgropen(opt2fn("-a2", NFILE, fnm), "Principal axis 2 (middle axis)",
+ output_env_get_xvgr_tlabel(oenv), "Component (nm)", oenv);
+ xvgr_legend(axis2, DIM, (const char **)legend, oenv);
+
+ axis3 = xvgropen(opt2fn("-a3", NFILE, fnm), "Principal axis 3 (minor axis)",
+ output_env_get_xvgr_tlabel(oenv), "Component (nm)", oenv);
+ xvgr_legend(axis3, DIM, (const char **)legend, oenv);
+
+ sprintf(legend[XX], "Axis 1 (major)");
+ sprintf(legend[YY], "Axis 2 (middle)");
+ sprintf(legend[ZZ], "Axis 3 (minor)");
+
+ fmoi = xvgropen(opt2fn("-om", NFILE, fnm), "Moments of inertia around inertial axes",
+ output_env_get_xvgr_tlabel(oenv), "I (au nm\\S2\\N)", oenv);
+ xvgr_legend(fmoi, DIM, (const char **)legend, oenv);
+
+ for (i = 0; i < DIM; i++)
+ {
+ sfree(legend[i]);
+ }
+ sfree(legend);
read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box, TRUE);
gmx_rmpbc_done(gpbc);
-
close_trj(status);
- gmx_ffclose(axis1);
- gmx_ffclose(axis2);
- gmx_ffclose(axis3);
- gmx_ffclose(fmoi);
+
+ xvgrclose(axis1);
+ xvgrclose(axis2);
+ xvgrclose(axis3);
+ xvgrclose(fmoi);
return 0;
}
#endif
#include <math.h>
-#include "sysstuff.h"
#include <string.h>
+
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "physics.h"
-#include "pbc.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
+#include "gromacs/math/units.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/commandline/pargs.h"
-#include "index.h"
#include "nrama.h"
#include "gmx_ana.h"
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
{
#include <config.h>
#endif
-#include "gromacs/utility/smalloc.h"
#include <math.h>
+#include <stdlib.h>
+
+#include "gromacs/utility/smalloc.h"
#include "macros.h"
#include "typedefs.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
#include "copyrite.h"
#include "gromacs/commandline/pargs.h"
-#include "vec.h"
-#include "index.h"
-#include "gmx_fatal.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
#include "princ.h"
-#include "rmpbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
#include "gromacs/fileio/matio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
int ePBC;
t_iatom *iatom = NULL;
- matrix box;
+ matrix box = {{0}};
rvec *x, *xp, *xm = NULL, **mat_x = NULL, **mat_x2, *mat_x2_j = NULL, vec1,
vec2;
t_trxstatus *status;
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++)
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
-#include <stdio.h>
#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
#include "gromacs/fileio/confio.h"
-#include "gmx_fatal.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
#include "gstat.h"
#include "macros.h"
#include "gromacs/math/utilities.h"
-#include "physics.h"
-#include "index.h"
+#include "gromacs/math/units.h"
+#include "gromacs/topology/index.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/commandline/pargs.h"
-#include <string.h>
-#include "sysstuff.h"
#include "txtdump.h"
#include "typedefs.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "pbc.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gmx_ana.h"
#include "gromacs/fileio/trxio.h"
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");
}
#endif
#include <math.h>
+#include <stdlib.h>
#include <string.h>
+
#include "gromacs/commandline/pargs.h"
-#include "sysstuff.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "vec.h"
-#include "pbc.h"
-#include "gromacs/fileio/futil.h"
-#include "index.h"
-#include "mshift.h"
-#include "xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "rmpbc.h"
-#include "physics.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/math/units.h"
#include "gromacs/fileio/confio.h"
#include "gmx_ana.h"
#include "gromacs/linearalgebra/nrjac.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
static void low_print_data(FILE *fp, real time, rvec x[], int n, atom_id *index,
gmx_bool bDim[], const char *sffmt)
read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
+
+ if ((bOV || bOF) && fn2ftp(ftp2fn(efTRX, NFILE, fnm)) == efXTC)
+ {
+ gmx_fatal(FARGS, "Cannot extract velocities or forces since your input XTC file does not contain them.");
+ }
+
if (bCV || bCF)
{
snew(sumx, fr.natoms);
#include <config.h>
#endif
-#include <string.h>
#include <math.h>
+#include <stdlib.h>
+#include <string.h>
#include "copyrite.h"
#include "macros.h"
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
#include "typedefs.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/trnio.h"
#include "gromacs/fileio/tngio_for_tools.h"
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/fileio/pdbio.h"
#include "gromacs/fileio/confio.h"
#include "names.h"
-#include "index.h"
-#include "vec.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/math/vec.h"
#include "gromacs/fileio/xtcio.h"
-#include "rmpbc.h"
-#include "pbc.h"
#include "viewit.h"
-#include "xvgr.h"
#include "gmx_ana.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/fileio/xvgr.h"
#include "gromacs/math/do_fit.h"
-#include "gmx_fatal.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
#ifdef HAVE_UNISTD_H
#include <unistd.h>
/* Make atoms struct for output in GRO or PDB files */
if ((ftp == efGRO) || ((ftp == efG96) && bTPS) || (ftp == efPDB))
{
- /* get memory for stuff to go in .pdb file */
- init_t_atoms(&useatoms, atoms->nr, FALSE);
+ /* get memory for stuff to go in .pdb file, and initialize
+ * the pdbinfo structure part if the input has it.
+ */
+ init_t_atoms(&useatoms, atoms->nr, (atoms->pdbinfo != NULL));
sfree(useatoms.resinfo);
useatoms.resinfo = atoms->resinfo;
for (i = 0; (i < nout); i++)
{
useatoms.atomname[i] = atoms->atomname[index[i]];
useatoms.atom[i] = atoms->atom[index[i]];
+ if (atoms->pdbinfo != NULL)
+ {
+ useatoms.pdbinfo[i] = atoms->pdbinfo[index[i]];
+ }
useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
}
useatoms.nr = nout;
#include <config.h>
#endif
+#include <assert.h>
#include <math.h>
+#include <stdlib.h>
#include <string.h>
-#include "sysstuff.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
#include "gromacs/commandline/pargs.h"
#include "gromacs/math/utilities.h"
-#include "gromacs/fileio/futil.h"
-#include "index.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/topology/index.h"
#include "typedefs.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
+#include "viewit.h"
#include "gstat.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "gromacs/fileio/matio.h"
#include "gmx_ana.h"
srenew(sbox, nalloc);
srenew(sx, nalloc);
}
+ assert(time != NULL); assert(sbox != NULL);
time[nfr] = t;
copy_mat(box, sbox[nfr]);
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));
#include "gromacs/commandline/pargs.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "copyrite.h"
#include "gromacs/fileio/tpxio.h"
#include "names.h"
#include "gmx_ana.h"
#include "macros.h"
#include "gromacs/utility/cstringutil.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
//! longest file names allowed in input files
#define WHAM_MAXFILELEN 2048
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++)
#include <config.h>
#endif
-#include <math.h>
#include <assert.h>
-#include "physics.h"
-#include "vec.h"
+#include <math.h>
+
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
#include "gromacs/math/utilities.h"
#include "txtdump.h"
#include "bondf.h"
-#include "gromacs/utility/smalloc.h"
-#include "pbc.h"
#include "ns.h"
#include "macros.h"
#include "names.h"
-#include "gmx_fatal.h"
-#include "mshift.h"
-#include "main.h"
#include "disre.h"
#include "orires.h"
#include "force.h"
#include "nonbonded.h"
#include "restcbt.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/simd/simd.h"
#include "gromacs/simd/simd_math.h"
#include "gromacs/simd/vector_operations.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
/* Find a better place for this? */
const int cmap_coeff_matrix[] = {
* 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)
/* The source code in this file should be thread-safe.
Please keep it that way. */
-
+#include "checkpoint.h"
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+#include <errno.h>
+#include <stdlib.h>
#include <string.h>
#include <time.h>
+#include <fcntl.h>
#ifdef HAVE_SYS_TIME_H
#include <sys/time.h>
#endif
#include "names.h"
#include "typedefs.h"
#include "types/commrec.h"
-#include "gromacs/utility/smalloc.h"
#include "txtdump.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "network.h"
-#include "checkpoint.h"
-#include "main.h"
-#include "gromacs/utility/cstringutil.h"
-#include <fcntl.h>
#include "gromacs/fileio/filenm.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xdrf.h"
#include "gromacs/fileio/xdr_datatype.h"
+#include "gromacs/utility/basenetwork.h"
#include "gromacs/utility/baseversion.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
#include "buildinfo.h"
return 0;
#else
FILE *fp;
- int rc;
fp = fopen(filename, "rb+");
{
if (dtc == xdr_datatype_double)
{
+ /* cppcheck-suppress invalidPointerCast
+ * Only executed if real is anyhow double */
vd = (double *)vp;
}
else
static int do_cpte_real(XDR *xd, int cptp, int ecpt, int sflags,
real *r, FILE *list)
{
- int n;
-
return do_cpte_reals_low(xd, cptp, ecpt, sflags, 1, NULL, &r, list, ecprREAL);
}
bool_t res = 0;
int dtc = xdr_datatype_int;
int *vp, *va = NULL;
- int nf, dt, i;
+ int nf, dt;
nf = n;
res = xdr_int(xd, &nf);
bool_t res = 0;
int dtc = xdr_datatype_double;
double *vp, *va = NULL;
- int nf, dt, i;
+ int nf, dt;
nf = n;
res = xdr_int(xd, &nf);
static int do_cpte_rvecs(XDR *xd, int cptp, int ecpt, int sflags,
int n, rvec **v, FILE *list)
{
- int n3;
-
return do_cpte_reals_low(xd, cptp, ecpt, sflags,
n*DIM, NULL, (real **)v, list, ecprRVEC);
}
matrix v, FILE *list)
{
real *vr;
- real ret;
+ int ret;
vr = (real *)&(v[0][0]);
ret = do_cpte_reals_low(xd, cptp, ecpt, sflags,
int n, real **v, FILE *list)
{
int i;
- real *vr;
- real ret, reti;
+ int ret, reti;
char name[CPTSTRLEN];
ret = 0;
}
for (i = 0; i < n; i++)
{
- reti = 0;
- vr = v[i];
reti = do_cpte_reals_low(xd, cptp, ecpt, sflags, n, NULL, &(v[i]), NULL, ecprREAL);
if (list && reti == 0)
{
sprintf(name, "%s[%d]", st_names(cptp, ecpt), i);
pr_reals(list, 0, name, v[i], n);
}
- if (reti == 0)
+ if (reti != 0)
{
- ret = 0;
+ ret = reti;
}
}
return ret;
bool_t res = 0;
int magic;
int idum = 0;
- int i;
char *fhost;
if (bRead)
{
enerhist->ener_sum_sim[i] = enerhist->ener_sum[i];
}
- fflags |= (1<<eenhENERGY_SUM_SIM);
}
if ( (fflags & (1<<eenhENERGY_NSUM)) &&
{
/* Assume we have an old file format and copy nsum to nsteps */
enerhist->nsteps = enerhist->nsum;
- fflags |= (1<<eenhENERGY_NSTEPS);
}
if ( (fflags & (1<<eenhENERGY_NSUM_SIM)) &&
!(fflags & (1<<eenhENERGY_NSTEPS_SIM)))
{
/* Assume we have an old file format and copy nsum to nsteps */
enerhist->nsteps_sim = enerhist->nsum_sim;
- fflags |= (1<<eenhENERGY_NSTEPS_SIM);
}
return ret;
static int do_cpt_EDstate(XDR *xd, gmx_bool bRead,
edsamstate_t *EDstate, FILE *list)
{
- int i, j;
+ int i;
int ret = 0;
char buf[STRLEN];
gmx_file_position_t **p_outputfiles, int *nfiles,
FILE *list, int file_version)
{
- int i, j;
+ int i;
gmx_off_t offset;
gmx_off_t mask = 0xFFFFFFFFL;
int offset_high, offset_low;
char *fntemp; /* the temporary checkpoint file name */
time_t now;
char timebuf[STRLEN];
- int nppnodes, npmenodes, flag_64bit;
+ int nppnodes, npmenodes;
char buf[1024], suffix[5+STEPSTRSIZE], sbuf[STEPSTRSIZE];
gmx_file_position_t *outputfiles;
int noutputfiles;
char *ftime;
- int flags_eks, flags_enh, flags_dfh, i;
+ int flags_eks, flags_enh, flags_dfh;
t_fileio *ret;
if (DOMAINDECOMP(cr))
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);
- sscanf(gmx_version(), "VERSION %d.%d", &gmx_major, &gmx_minor);
- sscanf(version, "VERSION %d.%d", &cpt_major, &cpt_minor);
+ 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 %5d.%5d", &gmx_major, &gmx_minor);
++ sscanf(version, "VERSION %5d.%5d", &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);
+ }
}
}
}
int file_version;
char *version, *btime, *buser, *bhost, *fprog, *ftime;
int double_prec;
- char filename[STRLEN], buf[STEPSTRSIZE];
- int nppnodes, eIntegrator_f, nppnodes_f, npmenodes_f;
+ char buf[STEPSTRSIZE];
+ int eIntegrator_f, nppnodes_f, npmenodes_f;
ivec dd_nc_f;
int natoms, ngtc, nnhpres, nhchainlength, nlambda, fflags, flags_eks, flags_enh, flags_dfh;
int d;
if (!PAR(cr))
{
- nppnodes = 1;
cr->npmenodes = 0;
}
else if (cr->nnodes == nppnodes_f + npmenodes_f)
{
cr->npmenodes = npmenodes_f;
}
- nppnodes = cr->nnodes - cr->npmenodes;
+ int nppnodes = cr->nnodes - cr->npmenodes;
if (nppnodes == nppnodes_f)
{
for (d = 0; d < DIM; d++)
}
}
}
- else
- {
- /* The number of PP nodes has not been set yet */
- nppnodes = -1;
- }
if (fflags != state->flags)
{
ivec dd_nc;
t_state state;
int flags_eks, flags_enh, flags_dfh;
- int indent;
- int i, j;
int ret;
gmx_file_position_t *outputfiles;
int nfiles;
}
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 <math.h>
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "typedefs.h"
#include "nonbonded.h"
#include "nb_kernel.h"
#include "macros.h"
#include "nb_free_energy.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
void
gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist,
real Vvdw6, Vvdw12, vvtot;
real ix, iy, iz, fix, fiy, fiz;
real dx, dy, dz, rsq, rinv;
- real c6[NSTATES], c12[NSTATES], c6grid[NSTATES];
+ real c6[NSTATES], c12[NSTATES], c6grid;
real LFC[NSTATES], LFV[NSTATES], DLF[NSTATES];
double dvdl_coul, dvdl_vdw;
real lfac_coul[NSTATES], dlfac_coul[NSTATES], lfac_vdw[NSTATES], dlfac_vdw[NSTATES];
real * Vv;
real * Vc;
gmx_bool bDoForces, bDoShiftForces, bDoPotential;
- real rcoulomb, sh_ewald;
- real rvdw, sh_invrc6;
- gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoffAll, bEwald;
+ real rcoulomb, rvdw, sh_invrc6;
+ gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoffAll;
+ gmx_bool bEwald, bEwaldLJ;
real rcutoff_max2;
- real rcutoff, rcutoff2, rswitch, d, d2, swV3, swV4, swV5, swF2, swF3, swF4, sw, dsw, rinvcorr;
- const real * tab_ewald_F;
- const real * tab_ewald_V;
const real * tab_ewald_F_lj;
const real * tab_ewald_V_lj;
- real tab_ewald_scale, tab_ewald_halfsp;
+ real d, d2, sw, dsw, rinvcorr;
+ real elec_swV3, elec_swV4, elec_swV5, elec_swF2, elec_swF3, elec_swF4;
+ real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
+ gmx_bool bConvertEwaldToCoulomb, bConvertLJEwaldToLJ6;
+ gmx_bool bComputeVdwInteraction, bComputeElecInteraction;
+ const real * ewtab;
+ int ewitab;
+ real ewrt, eweps, ewtabscale, ewtabhalfspace, sh_ewald;
+
+ sh_ewald = fr->ic->sh_ewald;
+ ewtab = fr->ic->tabq_coul_FDV0;
+ ewtabscale = fr->ic->tabq_scale;
+ ewtabhalfspace = 0.5/ewtabscale;
+ tab_ewald_F_lj = fr->ic->tabq_vdw_F;
+ tab_ewald_V_lj = fr->ic->tabq_vdw_V;
x = xx[0];
f = ff[0];
bDoPotential = kernel_data->flags & GMX_NONBONDED_DO_POTENTIAL;
rcoulomb = fr->rcoulomb;
-- sh_ewald = fr->ic->sh_ewald;
rvdw = fr->rvdw;
sh_invrc6 = fr->ic->sh_invrc6;
- /* Ewald (PME) reciprocal force and energy quadratic spline tables */
- tab_ewald_F = fr->ic->tabq_coul_F;
- tab_ewald_V = fr->ic->tabq_coul_V;
- tab_ewald_scale = fr->ic->tabq_scale;
- tab_ewald_F_lj = fr->ic->tabq_vdw_F;
- tab_ewald_V_lj = fr->ic->tabq_vdw_V;
- tab_ewald_halfsp = 0.5/tab_ewald_scale;
+ if (fr->coulomb_modifier == eintmodPOTSWITCH)
+ {
+ d = fr->rcoulomb-fr->rcoulomb_switch;
+ elec_swV3 = -10.0/(d*d*d);
+ elec_swV4 = 15.0/(d*d*d*d);
+ elec_swV5 = -6.0/(d*d*d*d*d);
+ elec_swF2 = -30.0/(d*d*d);
+ elec_swF3 = 60.0/(d*d*d*d);
+ elec_swF4 = -30.0/(d*d*d*d*d);
+ }
+ else
+ {
+ /* Avoid warnings from stupid compilers (looking at you, Clang!) */
+ elec_swV3 = elec_swV4 = elec_swV5 = elec_swF2 = elec_swF3 = elec_swF4 = 0.0;
+ }
- if (fr->coulomb_modifier == eintmodPOTSWITCH || fr->vdw_modifier == eintmodPOTSWITCH)
+ if (fr->vdw_modifier == eintmodPOTSWITCH)
{
- rcutoff = (fr->coulomb_modifier == eintmodPOTSWITCH) ? fr->rcoulomb : fr->rvdw;
- rcutoff2 = rcutoff*rcutoff;
- rswitch = (fr->coulomb_modifier == eintmodPOTSWITCH) ? fr->rcoulomb_switch : fr->rvdw_switch;
- d = rcutoff-rswitch;
- swV3 = -10.0/(d*d*d);
- swV4 = 15.0/(d*d*d*d);
- swV5 = -6.0/(d*d*d*d*d);
- swF2 = -30.0/(d*d*d);
- swF3 = 60.0/(d*d*d*d);
- swF4 = -30.0/(d*d*d*d*d);
+ d = fr->rvdw-fr->rvdw_switch;
+ vdw_swV3 = -10.0/(d*d*d);
+ vdw_swV4 = 15.0/(d*d*d*d);
+ vdw_swV5 = -6.0/(d*d*d*d*d);
+ vdw_swF2 = -30.0/(d*d*d);
+ vdw_swF3 = 60.0/(d*d*d*d);
+ vdw_swF4 = -30.0/(d*d*d*d*d);
}
else
{
- /* Stupid compilers dont realize these variables will not be used */
- rswitch = 0.0;
- swV3 = 0.0;
- swV4 = 0.0;
- swV5 = 0.0;
- swF2 = 0.0;
- swF3 = 0.0;
- swF4 = 0.0;
+ /* Avoid warnings from stupid compilers (looking at you, Clang!) */
+ vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
}
if (fr->cutoff_scheme == ecutsVERLET)
rcutoff_max2 = rcutoff_max2*rcutoff_max2;
bEwald = (icoul == GMX_NBKERNEL_ELEC_EWALD);
+ bEwaldLJ = (ivdw == GMX_NBKERNEL_VDW_LJEWALD);
+
+ /* For Ewald/PME interactions we cannot easily apply the soft-core component to
+ * reciprocal space. When we use vanilla (not switch/shift) Ewald interactions, we
+ * can apply the small trick of subtracting the _reciprocal_ space contribution
+ * in this kernel, and instead apply the free energy interaction to the 1/r
+ * (standard coulomb) interaction.
+ *
+ * However, we cannot use this approach for switch-modified since we would then
+ * effectively end up evaluating a significantly different interaction here compared to the
+ * normal (non-free-energy) kernels, either by applying a cutoff at a different
+ * position than what the user requested, or by switching different
+ * things (1/r rather than short-range Ewald). For these settings, we just
+ * use the traditional short-range Ewald interaction in that case.
+ */
+ bConvertEwaldToCoulomb = (bEwald && (fr->coulomb_modifier != eintmodPOTSWITCH));
+ /* For now the below will always be true (since LJ-PME only works with Shift in Gromacs-5.0),
+ * but writing it this way means we stay in sync with coulomb, and it avoids future bugs.
+ */
+ bConvertLJEwaldToLJ6 = (bEwaldLJ && (fr->vdw_modifier != eintmodPOTSWITCH));
/* fix compiler warnings */
nj1 = 0;
tj[STATE_A] = ntiA+2*typeA[jnr];
tj[STATE_B] = ntiB+2*typeB[jnr];
- if (ivdw == GMX_NBKERNEL_VDW_LJEWALD)
- {
- c6grid[STATE_A] = nbfp_grid[tj[STATE_A]];
- c6grid[STATE_B] = nbfp_grid[tj[STATE_B]];
- }
-
if (nlist->excl_fep == NULL || nlist->excl_fep[k])
{
c6[STATE_A] = nbfp[tj[STATE_A]];
n1V = tab_elemsize*n0;
}
- /* With Ewald and soft-core we should put the cut-off on r,
- * not on the soft-cored rC, as the real-space and
- * reciprocal space contributions should (almost) cancel.
+ /* Only process the coulomb interactions if we have charges,
+ * and if we either include all entries in the list (no cutoff
+ * used in the kernel), or if we are within the cutoff.
*/
- if (qq[i] != 0 &&
- !(bExactElecCutoff &&
- ((!bEwald && rC >= rcoulomb) ||
- (bEwald && r >= rcoulomb))))
+ bComputeElecInteraction = !bExactElecCutoff ||
+ ( bConvertEwaldToCoulomb && r < rcoulomb) ||
+ (!bConvertEwaldToCoulomb && rC < rcoulomb);
+
+ if ( (qq[i] != 0) && bComputeElecInteraction)
{
switch (icoul)
{
/* simple cutoff */
Vcoul[i] = qq[i]*rinvC;
FscalC[i] = Vcoul[i];
- break;
-
- case GMX_NBKERNEL_ELEC_EWALD:
- /* Ewald FEP is done only on the 1/r part */
- Vcoul[i] = qq[i]*(rinvC - sh_ewald);
- FscalC[i] = Vcoul[i];
+ /* The shift for the Coulomb potential is stored in
+ * the RF parameter c_rf, which is 0 without shift
+ */
+ Vcoul[i] -= qq[i]*fr->ic->c_rf;
break;
case GMX_NBKERNEL_ELEC_REACTIONFIELD:
gmx_fatal(FARGS, "Free energy and GB not implemented.\n");
break;
+ case GMX_NBKERNEL_ELEC_EWALD:
+ if (bConvertEwaldToCoulomb)
+ {
+ /* Ewald FEP is done only on the 1/r part */
+ Vcoul[i] = qq[i]*(rinvC-sh_ewald);
+ FscalC[i] = qq[i]*rinvC;
+ }
+ else
+ {
+ ewrt = rC*ewtabscale;
+ ewitab = (int) ewrt;
+ eweps = ewrt-ewitab;
+ ewitab = 4*ewitab;
+ FscalC[i] = ewtab[ewitab]+eweps*ewtab[ewitab+1];
+ rinvcorr = rinvC-sh_ewald;
+ Vcoul[i] = qq[i]*(rinvcorr-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+FscalC[i])));
+ FscalC[i] = qq[i]*(rinvC-rC*FscalC[i]);
+ }
+ break;
+
case GMX_NBKERNEL_ELEC_NONE:
FscalC[i] = 0.0;
Vcoul[i] = 0.0;
if (fr->coulomb_modifier == eintmodPOTSWITCH)
{
- d = rC-rswitch;
+ d = rC-fr->rcoulomb_switch;
d = (d > 0.0) ? d : 0.0;
d2 = d*d;
- sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
- dsw = d2*(swF2+d*(swF3+d*swF4));
+ sw = 1.0+d2*d*(elec_swV3+d*(elec_swV4+d*elec_swV5));
+ dsw = d2*(elec_swF2+d*(elec_swF3+d*elec_swF4));
- Vcoul[i] *= sw;
- FscalC[i] = FscalC[i]*sw + Vcoul[i]*dsw;
+ FscalC[i] = FscalC[i]*sw - rC*Vcoul[i]*dsw;
+ Vcoul[i] *= sw;
+
+ FscalC[i] = (rC < rcoulomb) ? FscalC[i] : 0.0;
+ Vcoul[i] = (rC < rcoulomb) ? Vcoul[i] : 0.0;
}
}
- if ((c6[i] != 0 || c12[i] != 0) &&
- !(bExactVdwCutoff &&
- ((ivdw != GMX_NBKERNEL_VDW_LJEWALD && rV >= rvdw) ||
- (ivdw == GMX_NBKERNEL_VDW_LJEWALD && r >= rvdw))))
+ /* Only process the VDW interactions if we have
+ * some non-zero parameters, and if we either
+ * include all entries in the list (no cutoff used
+ * in the kernel), or if we are within the cutoff.
+ */
+ bComputeVdwInteraction = !bExactVdwCutoff ||
+ ( bConvertLJEwaldToLJ6 && r < rvdw) ||
+ (!bConvertLJEwaldToLJ6 && rV < rvdw);
+ if ((c6[i] != 0 || c12[i] != 0) && bComputeVdwInteraction)
{
switch (ivdw)
{
if (fr->vdw_modifier == eintmodPOTSWITCH)
{
- d = rV-rswitch;
- d = (d > 0.0) ? d : 0.0;
- d2 = d*d;
- sw = 1.0+d2*d*(swV3+d*(swV4+d*swV5));
- dsw = d2*(swF2+d*(swF3+d*swF4));
+ d = rV-fr->rvdw_switch;
+ d = (d > 0.0) ? d : 0.0;
+ d2 = d*d;
+ sw = 1.0+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
+ dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
- Vvdw[i] *= sw;
- FscalV[i] = FscalV[i]*sw + Vvdw[i]*dsw;
+ FscalV[i] = FscalV[i]*sw - rV*Vvdw[i]*dsw;
+ Vvdw[i] *= sw;
FscalV[i] = (rV < rvdw) ? FscalV[i] : 0.0;
Vvdw[i] = (rV < rvdw) ? Vvdw[i] : 0.0;
}
}
- if (icoul == GMX_NBKERNEL_ELEC_EWALD &&
- !(bExactElecCutoff && r >= rcoulomb))
+ if (bConvertEwaldToCoulomb && ( !bExactElecCutoff || r < rcoulomb ) )
{
- /* Because we compute the soft-core normally,
- * we have to remove the Ewald short range portion.
- * Done outside of the states loop because this part
- * doesn't depend on the scaled R.
+ /* See comment in the preamble. When using Ewald interactions
+ * (unless we use a switch modifier) we subtract the reciprocal-space
+ * Ewald component here which made it possible to apply the free
+ * energy interaction to 1/r (vanilla coulomb short-range part)
+ * above. This gets us closer to the ideal case of applying
+ * the softcore to the entire electrostatic interaction,
+ * including the reciprocal-space component.
*/
- real rs, frac, f_lr;
- int ri;
+ real v_lr, f_lr;
- rs = rsq*rinv*tab_ewald_scale;
- ri = (int)rs;
- frac = rs - ri;
- f_lr = (1 - frac)*tab_ewald_F[ri] + frac*tab_ewald_F[ri+1];
- FF = f_lr*rinv;
- VV = tab_ewald_V[ri] - tab_ewald_halfsp*frac*(tab_ewald_F[ri] + f_lr);
+ ewrt = r*ewtabscale;
+ ewitab = (int) ewrt;
+ eweps = ewrt-ewitab;
+ ewitab = 4*ewitab;
+ f_lr = ewtab[ewitab]+eweps*ewtab[ewitab+1];
+ v_lr = (ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+f_lr));
+ f_lr *= rinv;
if (ii == jnr)
{
- VV *= 0.5;
+ /* If we get here, the i particle (ii) has itself (jnr)
+ * in its neighborlist. This can only happen with the Verlet
+ * scheme, and corresponds to a self-interaction that will
+ * occur twice. Scale it down by 50% to only include it once.
+ */
+ v_lr *= 0.5;
}
for (i = 0; i < NSTATES; i++)
{
- vctot -= LFC[i]*qq[i]*VV;
- Fscal -= LFC[i]*qq[i]*FF;
- dvdl_coul -= (DLF[i]*qq[i])*VV;
+ vctot -= LFC[i]*qq[i]*v_lr;
+ Fscal -= LFC[i]*qq[i]*f_lr;
+ dvdl_coul -= (DLF[i]*qq[i])*v_lr;
}
}
- if (ivdw == GMX_NBKERNEL_VDW_LJEWALD &&
- !(bExactVdwCutoff && r >= rvdw))
+ if (bConvertLJEwaldToLJ6 && (!bExactVdwCutoff || r < rvdw))
{
+ /* See comment in the preamble. When using LJ-Ewald interactions
+ * (unless we use a switch modifier) we subtract the reciprocal-space
+ * Ewald component here which made it possible to apply the free
+ * energy interaction to r^-6 (vanilla LJ6 short-range part)
+ * above. This gets us closer to the ideal case of applying
+ * the softcore to the entire VdW interaction,
+ * including the reciprocal-space component.
+ */
real rs, frac, f_lr;
int ri;
- rs = rsq*rinv*tab_ewald_scale;
+ rs = rsq*rinv*ewtabscale;
ri = (int)rs;
frac = rs - ri;
f_lr = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1];
FF = f_lr*rinv;
- VV = tab_ewald_V_lj[ri] - tab_ewald_halfsp*frac*(tab_ewald_F_lj[ri] + f_lr);
+ VV = tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr);
+
+ if (ii == jnr)
+ {
+ /* If we get here, the i particle (ii) has itself (jnr)
+ * in its neighborlist. This can only happen with the Verlet
+ * scheme, and corresponds to a self-interaction that will
+ * occur twice. Scale it down by 50% to only include it once.
+ */
+ VV *= 0.5;
+ }
+
for (i = 0; i < NSTATES; i++)
{
- vvtot += LFV[i]*c6grid[i]*VV*(1.0/6.0);
- Fscal += LFV[i]*c6grid[i]*FF*(1.0/6.0);
- dvdl_vdw += (DLF[i]*c6grid[i])*VV*(1.0/6.0);
+ c6grid = nbfp_grid[tj[i]];
+ vvtot += LFV[i]*c6grid*VV*(1.0/6.0);
+ Fscal += LFV[i]*c6grid*FF*(1.0/6.0);
+ dvdl_vdw += (DLF[i]*c6grid)*VV*(1.0/6.0);
}
}
#include <math.h>
#include "types/simple.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "typedefs.h"
#include "nb_generic.h"
#include "nrnb.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
#include "nonbonded.h"
#include "nb_kernel.h"
fvdw = (vvdw_rep - vvdw_disp - c6grid*(1.0/6.0)*exponent*ewclj6)*rinvsq;
if (fr->vdw_modifier == eintmodPOTSHIFT)
{
- vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion + c6grid*sh_lj_ewald)/6.0;
+ vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion - c6grid*sh_lj_ewald)/6.0;
}
else
{
#include "typedefs.h"
#include "txtdump.h"
-#include "gromacs/utility/smalloc.h"
#include "ns.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "gromacs/math/utilities.h"
#include "macros.h"
-#include "gromacs/utility/cstringutil.h"
#include "force.h"
#include "names.h"
-#include "main.h"
-#include "xvgr.h"
-#include "gmx_fatal.h"
-#include "physics.h"
#include "force.h"
#include "bondf.h"
#include "nrnb.h"
#include "nonbonded.h"
#include "gromacs/simd/simd.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+
#include "nb_kernel.h"
#include "nb_free_energy.h"
#include "nb_generic.h"
void
- gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl)
+ gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl, gmx_bool bElecAndVdwSwitchDiffers)
{
const char * elec;
const char * elec_mod;
}
}
- /* Give up. If this was a water kernel, leave the pointer as NULL, which
- * will disable water optimization in NS. If it is a particle kernel, set
- * the pointer to the generic NB kernel.
+ /* For now, the accelerated kernels cannot handle the combination of switch functions for both
+ * electrostatics and VdW that use different switch radius or switch cutoff distances
+ * (both of them enter in the switch function calculation). This would require
+ * us to evaluate two completely separate switch functions for every interaction.
+ * Instead, we disable such kernels by setting the pointer to NULL.
+ * This will cause the generic kernel (which can handle it) to be called instead.
+ *
+ * Note that we typically already enable tabulated coulomb interactions for this case,
+ * so this is mostly a safe-guard to make sure we call the generic kernel if the
+ * tables are disabled.
+ */
+ if ((nl->ielec != GMX_NBKERNEL_ELEC_NONE) && (nl->ielecmod == eintmodPOTSWITCH) &&
+ (nl->ivdw != GMX_NBKERNEL_VDW_NONE) && (nl->ivdwmod == eintmodPOTSWITCH) &&
+ bElecAndVdwSwitchDiffers)
+ {
+ nl->kernelptr_vf = NULL;
+ nl->kernelptr_f = NULL;
+ }
+
+ /* Give up, pick a generic one instead.
+ * We only do this for particle-particle kernels; by leaving the water-optimized kernel
+ * pointers to NULL, the water optimization will automatically be disabled for this interaction.
*/
if (nl->kernelptr_vf == NULL && !gmx_strcasecmp_min(geom, "Particle-Particle"))
{
fprintf(debug,
"WARNING - Slow generic NB kernel used for neighborlist with\n"
" Elec: '%s', Modifier: '%s'\n"
- " Vdw: '%s', Modifier: '%s'\n"
- " Geom: '%s', Other: '%s'\n\n",
- elec, elec_mod, vdw, vdw_mod, geom, other);
+ " Vdw: '%s', Modifier: '%s'\n",
+ elec, elec_mod, vdw, vdw_mod);
}
}
}
-
return;
}
nlist = nblists->nlist_sr;
f = f_shortrange;
}
- else if (range == 1)
+ else
{
/* Long-range */
if (!(flags & GMX_NONBONDED_DO_LR))
{
(*kernelptr)(&(nlist[i]), x, f, fr, mdatoms, &kernel_data, nrnb);
}
+ else
+ {
+ gmx_fatal(FARGS, "Non-empty neighborlist does not have any kernel pointer assigned.");
+ }
}
}
}
/* This file is completely threadsafe - please keep it that way! */
#include <stdio.h>
+#include <stdlib.h>
+
#include "typedefs.h"
#include "types/commrec.h"
#include "names.h"
#include "txtdump.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "macros.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
int pr_indent(FILE *fp, int n)
pr_ivec(fp, indent, "QMbasis", opts->QMbasis, opts->ngQM, FALSE);
pr_ivec(fp, indent, "QMcharge", opts->QMcharge, opts->ngQM, FALSE);
pr_ivec(fp, indent, "QMmult", opts->QMmult, opts->ngQM, FALSE);
- pr_bvec(fp, indent, "bSH", opts->bSH, opts->ngQM, FALSE);
+ pr_bvec(fp, indent, "SH", opts->bSH, opts->ngQM, FALSE);
pr_ivec(fp, indent, "CASorbitals", opts->CASorbitals, opts->ngQM, FALSE);
pr_ivec(fp, indent, "CASelectrons", opts->CASelectrons, opts->ngQM, FALSE);
pr_rvec(fp, indent, "SAon", opts->SAon, opts->ngQM, FALSE);
- pr_rvec(fp, indent, "SAon", opts->SAon, opts->ngQM, FALSE);
+ pr_rvec(fp, indent, "SAoff", opts->SAoff, opts->ngQM, FALSE);
pr_ivec(fp, indent, "SAsteps", opts->SAsteps, opts->ngQM, FALSE);
pr_bvec(fp, indent, "bOPT", opts->bOPT, opts->ngQM, FALSE);
pr_bvec(fp, indent, "bTS", opts->bTS, opts->ngQM, FALSE);
fprintf(out, "\n");
/* Pretty-print the simulated annealing info */
- fprintf(out, "anneal%s", bMDPformat ? " = " : ":");
+ fprintf(out, "annealing%s", bMDPformat ? " = " : ":");
for (i = 0; (i < opts->ngtc); i++)
{
fprintf(out, " %10s", EANNEAL(opts->annealing[i]));
}
fprintf(out, "\n");
- fprintf(out, "ann-npoints%s", bMDPformat ? " = " : ":");
+ fprintf(out, "annealing-npoints%s", bMDPformat ? " = " : ":");
for (i = 0; (i < opts->ngtc); i++)
{
fprintf(out, " %10d", opts->anneal_npoints[i]);
{
if (opts->anneal_npoints[i] > 0)
{
- fprintf(out, "ann. times [%d]:\t", i);
+ fprintf(out, "annealing-time [%d]:\t", i);
for (j = 0; (j < opts->anneal_npoints[i]); j++)
{
fprintf(out, " %10.1f", opts->anneal_time[i][j]);
}
fprintf(out, "\n");
- fprintf(out, "ann. temps [%d]:\t", i);
+ fprintf(out, "annealing-temp [%d]:\t", i);
for (j = 0; (j < opts->anneal_npoints[i]); j++)
{
fprintf(out, " %10.1f", opts->anneal_temp[i][j]);
static void pr_simtempvals(FILE *fp, int indent, t_simtemp *simtemp, int n_lambda)
{
- PR("simtemp_low", simtemp->simtemp_low);
- PR("simtemp_high", simtemp->simtemp_high);
PS("simulated-tempering-scaling", ESIMTEMP(simtemp->eSimTempScale));
+ PR("sim-temp-low", simtemp->simtemp_low);
+ PR("sim-temp-high", simtemp->simtemp_high);
pr_rvec(fp, indent, "simulated tempering temperatures", simtemp->temperatures, n_lambda, TRUE);
}
{
PI("nstexpanded", expand->nstexpanded);
- PS("lambda-stats", elamstats_names[expand->elamstats]);
- PS("lambda-mc-move", elmcmove_names[expand->elmcmove]);
- PI("lmc-repeats", expand->lmc_repeats);
- PI("lmc-gibbsdelta", expand->gibbsdeltalam);
- PI("lmc-nstart", expand->lmc_forced_nstart);
- PS("symmetrized-transition-matrix", EBOOL(expand->bSymmetrizedTMatrix));
- PI("nst-transition-matrix", expand->nstTij);
- PI("mininum-var-min", expand->minvarmin); /*default is reasonable */
- PI("weight-c-range", expand->c_range); /* default is just C=0 */
- PR("wl-scale", expand->wl_scale);
- PR("init-wl-delta", expand->init_wl_delta);
- PR("wl-ratio", expand->wl_ratio);
- PS("bWLoneovert", EBOOL(expand->bWLoneovert));
- PI("lmc-seed", expand->lmc_seed);
- PR("mc-temperature", expand->mc_temp);
+ PS("lmc-stats", elamstats_names[expand->elamstats]);
+ PS("lmc-move", elmcmove_names[expand->elmcmove]);
PS("lmc-weights-equil", elmceq_names[expand->elmceq]);
if (expand->elmceq == elmceqNUMATLAM)
{
{
PR("weight-equil-count-ratio", expand->equil_ratio);
}
+ PI("lmc-seed", expand->lmc_seed);
+ PR("mc-temperature", expand->mc_temp);
+ PI("lmc-repeats", expand->lmc_repeats);
+ PI("lmc-gibbsdelta", expand->gibbsdeltalam);
+ PI("lmc-forced-nstart", expand->lmc_forced_nstart);
+ PS("symmetrized-transition-matrix", EBOOL(expand->bSymmetrizedTMatrix));
+ PI("nst-transition-matrix", expand->nstTij);
+ PI("mininum-var-min", expand->minvarmin); /*default is reasonable */
+ PI("weight-c-range", expand->c_range); /* default is just C=0 */
+ PR("wl-scale", expand->wl_scale);
+ PR("wl-ratio", expand->wl_ratio);
+ PR("init-wl-delta", expand->init_wl_delta);
+ PS("wl-oneovert", EBOOL(expand->bWLoneovert));
pr_indent(fp, indent);
pr_rvec(fp, indent, "init-lambda-weights", expand->init_lambda_weights, n_lambda, TRUE);
{
int i, j;
- PI("nstdhdl", fep->nstdhdl);
- PI("init-lambda-state", fep->init_fep_state);
PR("init-lambda", fep->init_lambda);
+ PI("init-lambda-state", fep->init_fep_state);
PR("delta-lambda", fep->delta_lambda);
+ PI("nstdhdl", fep->nstdhdl);
+
if (!bMDPformat)
{
PI("n-lambdas", fep->n_lambda);
}
}
PI("calc-lambda-neighbors", fep->lambda_neighbors);
-
+ PS("dhdl-print-energy", EBOOL(fep->bPrintEnergy));
PR("sc-alpha", fep->sc_alpha);
- PS("bScCoul", EBOOL(fep->bScCoul));
- PS("bScPrintEnergy", EBOOL(fep->bPrintEnergy));
PI("sc-power", fep->sc_power);
PR("sc-r-power", fep->sc_r_power);
PR("sc-sigma", fep->sc_sigma);
PR("sc-sigma-min", fep->sc_sigma_min);
- PS("separate-dhdl-file", SEPDHDLFILETYPE(fep->separate_dhdl_file));
- PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
+ PS("sc-coul", EBOOL(fep->bScCoul));
PI("dh-hist-size", fep->dh_hist_size);
PD("dh-hist-spacing", fep->dh_hist_spacing);
+ PS("separate-dhdl-file", SEPDHDLFILETYPE(fep->separate_dhdl_file));
+ PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
};
static void pr_pull(FILE *fp, int indent, t_pull *pull)
PR("pull-r1", pull->cyl_r1);
PR("pull-r0", pull->cyl_r0);
PR("pull-constr-tol", pull->constr_tol);
- PS("pull-bPrintRef", EBOOL(pull->bPrintRef));
+ PS("pull-print-reference", EBOOL(pull->bPrintRef));
PI("pull-nstxout", pull->nstxout);
PI("pull-nstfout", pull->nstfout);
- PI("pull-ngroup", pull->ngroup);
+ PI("pull-ngroups", pull->ngroup);
for (g = 0; g < pull->ngroup; g++)
{
pr_pull_group(fp, indent, g, &pull->group[g]);
}
- PI("pull-ncoord", pull->ncoord);
+ PI("pull-ncoords", pull->ncoord);
for (g = 0; g < pull->ncoord; g++)
{
pr_pull_coord(fp, indent, g, &pull->coord[g]);
static void pr_rotgrp(FILE *fp, int indent, int g, t_rotgrp *rotg)
{
pr_indent(fp, indent);
- fprintf(fp, "rotation_group %d:\n", g);
+ fprintf(fp, "rot-group %d:\n", g);
indent += 2;
- PS("type", EROTGEOM(rotg->eType));
- PS("massw", EBOOL(rotg->bMassW));
+ PS("rot-type", EROTGEOM(rotg->eType));
+ PS("rot-massw", EBOOL(rotg->bMassW));
pr_ivec_block(fp, indent, "atom", rotg->ind, rotg->nat, TRUE);
- pr_rvecs(fp, indent, "x_ref", rotg->x_ref, rotg->nat);
- pr_rvec(fp, indent, "vec", rotg->vec, DIM, TRUE);
- pr_rvec(fp, indent, "pivot", rotg->pivot, DIM, TRUE);
- PR("rate", rotg->rate);
- PR("k", rotg->k);
- PR("slab_dist", rotg->slab_dist);
- PR("min_gaussian", rotg->min_gaussian);
- PR("epsilon", rotg->eps);
- PS("fit_method", EROTFIT(rotg->eFittype));
- PI("potfitangle_nstep", rotg->PotAngle_nstep);
- PR("potfitangle_step", rotg->PotAngle_step);
+ pr_rvecs(fp, indent, "x-ref", rotg->x_ref, rotg->nat);
+ pr_rvec(fp, indent, "rot-vec", rotg->vec, DIM, TRUE);
+ pr_rvec(fp, indent, "rot-pivot", rotg->pivot, DIM, TRUE);
+ PR("rot-rate", rotg->rate);
+ PR("rot-k", rotg->k);
+ PR("rot-slab-dist", rotg->slab_dist);
+ PR("rot-min-gauss", rotg->min_gaussian);
+ PR("rot-eps", rotg->eps);
+ PS("rot-fit-method", EROTFIT(rotg->eFittype));
+ PI("rot_potfit_nstep", rotg->PotAngle_nstep);
+ PR("rot_potfit_step", rotg->PotAngle_step);
}
static void pr_rot(FILE *fp, int indent, t_rot *rot)
{
int g;
- PI("rot_nstrout", rot->nstrout);
- PI("rot_nstsout", rot->nstsout);
- PI("rot_ngrp", rot->ngrp);
+ PI("rot-nstrout", rot->nstrout);
+ PI("rot-nstsout", rot->nstsout);
+ PI("rot-ngroups", rot->ngrp);
for (g = 0; g < rot->ngrp; g++)
{
pr_rotgrp(fp, indent, g, &rot->grp[g]);
char str[STRLEN];
- PI("frequency", swap->nstswap);
- for (j = 0; j < 2; j++)
- {
- sprintf(str, "nanions%c", j+'A');
- PI(str, swap->nanions[j]);
- sprintf(str, "ncations%c", j+'A');
- PI(str, swap->ncations[j]);
- }
- PI("coupling_steps", swap->nAverage);
- PR("threshold", swap->threshold);
+ PI("swap-frequency", swap->nstswap);
for (j = 0; j < 2; j++)
{
- sprintf(str, "splitgroup%d_massw", j);
+ sprintf(str, "massw_split%d", j);
PS(str, EBOOL(swap->massw_split[j]));
sprintf(str, "split atoms group %d", j);
pr_ivec_block(fp, indent, str, swap->ind_split[j], swap->nat_split[j], TRUE);
}
pr_ivec_block(fp, indent, "swap atoms", swap->ind, swap->nat, TRUE);
pr_ivec_block(fp, indent, "solvent atoms", swap->ind_sol, swap->nat_sol, TRUE);
- PR("cyl0_radius", swap->cyl0r);
- PR("cyl0_upper", swap->cyl0u);
- PR("cyl0_lower", swap->cyl0l);
- PR("cyl1_radius", swap->cyl1r);
- PR("cyl1_upper", swap->cyl1u);
- PR("cyl1_lower", swap->cyl1l);
+ PR("cyl0-r", swap->cyl0r);
+ PR("cyl0-up", swap->cyl0u);
+ PR("cyl0-down", swap->cyl0l);
+ PR("cyl1-r", swap->cyl1r);
+ PR("cyl1-up", swap->cyl1u);
+ PR("cyl1-down", swap->cyl1l);
+ PI("coupl-steps", swap->nAverage);
+ for (j = 0; j < 2; j++)
+ {
+ sprintf(str, "anions%c", j+'A');
+ PI(str, swap->nanions[j]);
+ sprintf(str, "cations%c", j+'A');
+ PI(str, swap->ncations[j]);
+ }
+ PR("threshold", swap->threshold);
}
static void pr_imd(FILE *fp, int indent, t_IMD *imd)
{
- PI("IMD_atoms", imd->nat);
+ PI("IMD-atoms", imd->nat);
pr_ivec_block(fp, indent, "atom", imd->ind, imd->nat, TRUE);
}
{
indent = pr_title(fp, indent, title);
}
- /* This strings do not all have a direct correspondence to
- .mdp entries, but we should follow the same convention of
- using hyphens in the names users read and write. */
+ /* Try to make this list appear in the same order as the
+ * options are written in the default mdout.mdp, and with
+ * the same user-exposed names to facilitate debugging.
+ */
PS("integrator", EI(ir->eI));
+ PR("tinit", ir->init_t);
+ PR("dt", ir->delta_t);
PSTEP("nsteps", ir->nsteps);
PSTEP("init-step", ir->init_step);
- PS("cutoff-scheme", ECUTSCHEME(ir->cutoff_scheme));
- PS("ns-type", ENS(ir->ns_type));
- PI("nstlist", ir->nstlist);
- PI("ndelta", ir->ndelta);
- PI("nstcomm", ir->nstcomm);
+ PI("simulation-part", ir->simulation_part);
PS("comm-mode", ECOM(ir->comm_mode));
- PI("nstlog", ir->nstlog);
+ PI("nstcomm", ir->nstcomm);
+
+ /* Langevin dynamics */
+ PR("bd-fric", ir->bd_fric);
+ PSTEP("ld-seed", ir->ld_seed);
+
+ /* Energy minimization */
+ PR("emtol", ir->em_tol);
+ PR("emstep", ir->em_stepsize);
+ PI("niter", ir->niter);
+ PR("fcstep", ir->fc_stepsize);
+ PI("nstcgsteep", ir->nstcgsteep);
+ PI("nbfgscorr", ir->nbfgscorr);
+
+ /* Test particle insertion */
+ PR("rtpi", ir->rtpi);
+
+ /* Output control */
PI("nstxout", ir->nstxout);
PI("nstvout", ir->nstvout);
PI("nstfout", ir->nstfout);
+ PI("nstlog", ir->nstlog);
PI("nstcalcenergy", ir->nstcalcenergy);
PI("nstenergy", ir->nstenergy);
PI("nstxout-compressed", ir->nstxout_compressed);
- PR("init-t", ir->init_t);
- PR("delta-t", ir->delta_t);
+ PR("compressed-x-precision", ir->x_compression_precision);
- PR("x-compression-precision", ir->x_compression_precision);
- PR("fourierspacing", ir->fourier_spacing);
- PI("nkx", ir->nkx);
- PI("nky", ir->nky);
- PI("nkz", ir->nkz);
- PI("pme-order", ir->pme_order);
- PR("ewald-rtol", ir->ewald_rtol);
- PR("ewald-rtol-lj", ir->ewald_rtol_lj);
- PR("ewald-geometry", ir->ewald_geometry);
- PR("epsilon-surface", ir->epsilon_surface);
- PS("optimize-fft", EBOOL(ir->bOptFFT));
- PS("lj-pme-comb-rule", ELJPMECOMBNAMES(ir->ljpme_combination_rule));
- PS("ePBC", EPBC(ir->ePBC));
- PS("bPeriodicMols", EBOOL(ir->bPeriodicMols));
- PS("bContinuation", EBOOL(ir->bContinuation));
- PS("bShakeSOR", EBOOL(ir->bShakeSOR));
- PS("etc", ETCOUPLTYPE(ir->etc));
- PS("bPrintNHChains", EBOOL(ir->bPrintNHChains));
- PI("nsttcouple", ir->nsttcouple);
- PS("epc", EPCOUPLTYPE(ir->epc));
- PS("epctype", EPCOUPLTYPETYPE(ir->epct));
- PI("nstpcouple", ir->nstpcouple);
- PR("tau-p", ir->tau_p);
- pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
- pr_matrix(fp, indent, "compress", ir->compress, bMDPformat);
- PS("refcoord-scaling", EREFSCALINGTYPE(ir->refcoord_scaling));
- if (bMDPformat)
- {
- fprintf(fp, "posres-com = %g %g %g\n", ir->posres_com[XX],
- ir->posres_com[YY], ir->posres_com[ZZ]);
- }
- else
- {
- pr_rvec(fp, indent, "posres-com", ir->posres_com, DIM, TRUE);
- }
- if (bMDPformat)
- {
- fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX],
- ir->posres_comB[YY], ir->posres_comB[ZZ]);
- }
- else
- {
- pr_rvec(fp, indent, "posres-comB", ir->posres_comB, DIM, TRUE);
- }
+ /* Neighborsearching parameters */
+ PS("cutoff-scheme", ECUTSCHEME(ir->cutoff_scheme));
+ PI("nstlist", ir->nstlist);
+ PS("ns-type", ENS(ir->ns_type));
+ PS("pbc", EPBC(ir->ePBC));
+ PS("periodic-molecules", EBOOL(ir->bPeriodicMols));
PR("verlet-buffer-tolerance", ir->verletbuf_tol);
PR("rlist", ir->rlist);
PR("rlistlong", ir->rlistlong);
PR("nstcalclr", ir->nstcalclr);
- PR("rtpi", ir->rtpi);
+
+ /* Options for electrostatics and VdW */
PS("coulombtype", EELTYPE(ir->coulombtype));
PS("coulomb-modifier", INTMODIFIER(ir->coulomb_modifier));
PR("rcoulomb-switch", ir->rcoulomb_switch);
PR("rcoulomb", ir->rcoulomb);
- PS("vdwtype", EVDWTYPE(ir->vdwtype));
- PS("vdw-modifier", INTMODIFIER(ir->vdw_modifier));
- PR("rvdw-switch", ir->rvdw_switch);
- PR("rvdw", ir->rvdw);
if (ir->epsilon_r != 0)
{
PR("epsilon-r", ir->epsilon_r);
{
PS("epsilon-rf", infbuf);
}
- PR("tabext", ir->tabext);
+ PS("vdw-type", EVDWTYPE(ir->vdwtype));
+ PS("vdw-modifier", INTMODIFIER(ir->vdw_modifier));
+ PR("rvdw-switch", ir->rvdw_switch);
+ PR("rvdw", ir->rvdw);
+ PS("DispCorr", EDISPCORR(ir->eDispCorr));
+ PR("table-extension", ir->tabext);
+
+ PR("fourierspacing", ir->fourier_spacing);
+ PI("fourier-nx", ir->nkx);
+ PI("fourier-ny", ir->nky);
+ PI("fourier-nz", ir->nkz);
+ PI("pme-order", ir->pme_order);
+ PR("ewald-rtol", ir->ewald_rtol);
+ PR("ewald-rtol-lj", ir->ewald_rtol_lj);
+ PS("lj-pme-comb-rule", ELJPMECOMBNAMES(ir->ljpme_combination_rule));
+ PR("ewald-geometry", ir->ewald_geometry);
+ PR("epsilon-surface", ir->epsilon_surface);
+
+ /* Implicit solvent */
PS("implicit-solvent", EIMPLICITSOL(ir->implicit_solvent));
+
+ /* Generalized born electrostatics */
PS("gb-algorithm", EGBALGORITHM(ir->gb_algorithm));
- PR("gb-epsilon-solvent", ir->gb_epsilon_solvent);
PI("nstgbradii", ir->nstgbradii);
PR("rgbradii", ir->rgbradii);
+ PR("gb-epsilon-solvent", ir->gb_epsilon_solvent);
PR("gb-saltconc", ir->gb_saltconc);
PR("gb-obc-alpha", ir->gb_obc_alpha);
PR("gb-obc-beta", ir->gb_obc_beta);
PR("gb-dielectric-offset", ir->gb_dielectric_offset);
PS("sa-algorithm", ESAALGORITHM(ir->gb_algorithm));
PR("sa-surface-tension", ir->sa_surface_tension);
- PS("DispCorr", EDISPCORR(ir->eDispCorr));
- PS("bSimTemp", EBOOL(ir->bSimTemp));
- if (ir->bSimTemp)
- {
- pr_simtempvals(fp, indent, ir->simtempvals, ir->fepvals->n_lambda);
- }
- PS("free-energy", EFEPTYPE(ir->efep));
- if (ir->efep != efepNO || ir->bSimTemp)
+
+ /* Options for weak coupling algorithms */
+ PS("tcoupl", ETCOUPLTYPE(ir->etc));
+ PI("nsttcouple", ir->nsttcouple);
+ PI("nh-chain-length", ir->opts.nhchainlength);
+ PS("print-nose-hoover-chain-variables", EBOOL(ir->bPrintNHChains));
+
+ PS("pcoupl", EPCOUPLTYPE(ir->epc));
+ PS("pcoupltype", EPCOUPLTYPETYPE(ir->epct));
+ PI("nstpcouple", ir->nstpcouple);
+ PR("tau-p", ir->tau_p);
+ pr_matrix(fp, indent, "compressibility", ir->compress, bMDPformat);
+ pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
+ PS("refcoord-scaling", EREFSCALINGTYPE(ir->refcoord_scaling));
+
+ if (bMDPformat)
{
- pr_fepvals(fp, indent, ir->fepvals, bMDPformat);
+ fprintf(fp, "posres-com = %g %g %g\n", ir->posres_com[XX],
+ ir->posres_com[YY], ir->posres_com[ZZ]);
+ fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX],
+ ir->posres_comB[YY], ir->posres_comB[ZZ]);
}
- if (ir->bExpanded)
+ else
{
- pr_expandedvals(fp, indent, ir->expandedvals, ir->fepvals->n_lambda);
+ pr_rvec(fp, indent, "posres-com", ir->posres_com, DIM, TRUE);
+ pr_rvec(fp, indent, "posres-comB", ir->posres_comB, DIM, TRUE);
}
+ /* QMMM */
+ PS("QMMM", EBOOL(ir->bQMMM));
+ PI("QMconstraints", ir->QMconstraints);
+ PI("QMMMscheme", ir->QMMMscheme);
+ PR("MMChargeScaleFactor", ir->scalefactor);
+ pr_qm_opts(fp, indent, "qm-opts", &(ir->opts));
+
+ /* CONSTRAINT OPTIONS */
+ PS("constraint-algorithm", ECONSTRTYPE(ir->eConstrAlg));
+ PS("continuation", EBOOL(ir->bContinuation));
+
+ PS("Shake-SOR", EBOOL(ir->bShakeSOR));
+ PR("shake-tol", ir->shake_tol);
+ PI("lincs-order", ir->nProjOrder);
+ PI("lincs-iter", ir->nLincsIter);
+ PR("lincs-warnangle", ir->LincsWarnAngle);
+
+ /* Walls */
PI("nwall", ir->nwall);
PS("wall-type", EWALLTYPE(ir->wall_type));
+ PR("wall-r-linpot", ir->wall_r_linpot);
+ /* wall-atomtype */
PI("wall-atomtype[0]", ir->wall_atomtype[0]);
PI("wall-atomtype[1]", ir->wall_atomtype[1]);
+ /* wall-density */
PR("wall-density[0]", ir->wall_density[0]);
PR("wall-density[1]", ir->wall_density[1]);
PR("wall-ewald-zfac", ir->wall_ewald_zfac);
+ /* COM PULLING */
PS("pull", EPULLTYPE(ir->ePull));
if (ir->ePull != epullNO)
{
pr_pull(fp, indent, ir->pull);
}
+ /* ENFORCED ROTATION */
PS("rotation", EBOOL(ir->bRot));
if (ir->bRot)
{
pr_rot(fp, indent, ir->rot);
}
+ /* INTERACTIVE MD */
PS("interactiveMD", EBOOL(ir->bIMD));
if (ir->bIMD)
{
pr_imd(fp, indent, ir->imd);
}
+ /* NMR refinement stuff */
PS("disre", EDISRETYPE(ir->eDisre));
PS("disre-weighting", EDISREWEIGHTING(ir->eDisreWeighting));
PS("disre-mixed", EBOOL(ir->bDisreMixed));
PR("dr-fc", ir->dr_fc);
PR("dr-tau", ir->dr_tau);
PR("nstdisreout", ir->nstdisreout);
- PR("orires-fc", ir->orires_fc);
- PR("orires-tau", ir->orires_tau);
- PR("nstorireout", ir->nstorireout);
- PR("dihre-fc", ir->dihre_fc);
+ PR("orire-fc", ir->orires_fc);
+ PR("orire-tau", ir->orires_tau);
+ PR("nstorireout", ir->nstorireout);
- PR("em-stepsize", ir->em_stepsize);
- PR("em-tol", ir->em_tol);
- PI("niter", ir->niter);
- PR("fc-stepsize", ir->fc_stepsize);
- PI("nstcgsteep", ir->nstcgsteep);
- PI("nbfgscorr", ir->nbfgscorr);
+ /* FREE ENERGY VARIABLES */
+ PS("free-energy", EFEPTYPE(ir->efep));
+ if (ir->efep != efepNO || ir->bSimTemp)
+ {
+ pr_fepvals(fp, indent, ir->fepvals, bMDPformat);
+ }
+ if (ir->bExpanded)
+ {
+ pr_expandedvals(fp, indent, ir->expandedvals, ir->fepvals->n_lambda);
+ }
- PS("ConstAlg", ECONSTRTYPE(ir->eConstrAlg));
- PR("shake-tol", ir->shake_tol);
- PI("lincs-order", ir->nProjOrder);
- PR("lincs-warnangle", ir->LincsWarnAngle);
- PI("lincs-iter", ir->nLincsIter);
- PR("bd-fric", ir->bd_fric);
- PSTEP("ld-seed", ir->ld_seed);
- PR("cos-accel", ir->cos_accel);
+ /* NON-equilibrium MD stuff */
+ PR("cos-acceleration", ir->cos_accel);
pr_matrix(fp, indent, "deform", ir->deform, bMDPformat);
+ /* SIMULATED TEMPERING */
+ PS("simulated-tempering", EBOOL(ir->bSimTemp));
+ if (ir->bSimTemp)
+ {
+ pr_simtempvals(fp, indent, ir->simtempvals, ir->fepvals->n_lambda);
+ }
+
+ /* ELECTRIC FIELDS */
+ pr_cosine(fp, indent, "E-x", &(ir->ex[XX]), bMDPformat);
+ pr_cosine(fp, indent, "E-xt", &(ir->et[XX]), bMDPformat);
+ pr_cosine(fp, indent, "E-y", &(ir->ex[YY]), bMDPformat);
+ pr_cosine(fp, indent, "E-yt", &(ir->et[YY]), bMDPformat);
+ pr_cosine(fp, indent, "E-z", &(ir->ex[ZZ]), bMDPformat);
+ pr_cosine(fp, indent, "E-zt", &(ir->et[ZZ]), bMDPformat);
+
+ /* ION/WATER SWAPPING FOR COMPUTATIONAL ELECTROPHYSIOLOGY */
+ PS("swapcoords", ESWAPTYPE(ir->eSwapCoords));
+ if (ir->eSwapCoords != eswapNO)
+ {
+ pr_swap(fp, indent, ir->swap);
+ }
+
+ /* AdResS PARAMETERS */
PS("adress", EBOOL(ir->bAdress));
if (ir->bAdress)
{
PR("adress-const-wf", ir->adress->const_wf);
PR("adress-ex-width", ir->adress->ex_width);
PR("adress-hy-width", ir->adress->hy_width);
+ PR("adress-ex-forcecap", ir->adress->ex_forcecap);
PS("adress-interface-correction", EADRESSICTYPE(ir->adress->icor));
PS("adress-site", EADRESSSITETYPE(ir->adress->site));
- PR("adress-ex-force-cap", ir->adress->ex_forcecap);
- PS("adress-do-hybridpairs", EBOOL(ir->adress->do_hybridpairs));
-
pr_rvec(fp, indent, "adress-reference-coords", ir->adress->refs, DIM, TRUE);
+ PS("adress-do-hybridpairs", EBOOL(ir->adress->do_hybridpairs));
}
+
+ /* USER-DEFINED THINGIES */
PI("userint1", ir->userint1);
PI("userint2", ir->userint2);
PI("userint3", ir->userint3);
PR("userreal2", ir->userreal2);
PR("userreal3", ir->userreal3);
PR("userreal4", ir->userreal4);
+
pr_grp_opts(fp, indent, "grpopts", &(ir->opts), bMDPformat);
- pr_cosine(fp, indent, "efield-x", &(ir->ex[XX]), bMDPformat);
- pr_cosine(fp, indent, "efield-xt", &(ir->et[XX]), bMDPformat);
- pr_cosine(fp, indent, "efield-y", &(ir->ex[YY]), bMDPformat);
- pr_cosine(fp, indent, "efield-yt", &(ir->et[YY]), bMDPformat);
- pr_cosine(fp, indent, "efield-z", &(ir->ex[ZZ]), bMDPformat);
- pr_cosine(fp, indent, "efield-zt", &(ir->et[ZZ]), bMDPformat);
- PS("eSwapCoords", ESWAPTYPE(ir->eSwapCoords));
- if (ir->eSwapCoords != eswapNO)
- {
- pr_swap(fp, indent, ir->swap);
- }
- PS("bQMMM", EBOOL(ir->bQMMM));
- PI("QMconstraints", ir->QMconstraints);
- PI("QMMMscheme", ir->QMMMscheme);
- PR("scalefactor", ir->scalefactor);
- pr_qm_opts(fp, indent, "qm-opts", &(ir->opts));
}
}
#undef PS
for (i = 0; i <= block->nr; i++)
{
(void) pr_indent(fp, indent+INDENT);
- (void) fprintf(fp, "%s->index[%d]=%u\n",
+ (void) fprintf(fp, "%s->index[%d]=%d\n",
title, bShowNumbers ? i : -1, block->index[i]);
}
for (i = 0; i < block->nra; i++)
{
(void) pr_indent(fp, indent+INDENT);
- (void) fprintf(fp, "%s->a[%d]=%u\n",
+ (void) fprintf(fp, "%s->a[%d]=%d\n",
title, bShowNumbers ? i : -1, block->a[i]);
}
}
(void) fprintf(fp, "\n");
size = pr_indent(fp, indent+INDENT);
}
- size += fprintf(fp, "%u", block->a[j]);
+ size += fprintf(fp, "%d", block->a[j]);
}
(void) fprintf(fp, "}\n");
start = end;
#include <stdlib.h>
#include <string.h>
-#include "vec.h"
+
+#include "gromacs/math/vec.h"
#include "macros.h"
-#include "gromacs/utility/smalloc.h"
#include "types/commrec.h"
#include "force.h"
#include "names.h"
#include "mdatoms.h"
#include "nrnb.h"
#include "ns.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "chargegroup.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/smalloc.h"
+
static real box_margin;
static real max_dist(rvec *x, real *r, int start, int end)
ir->vdw_modifier = eintmodNONE;
ir->coulombtype = eelCUT;
ir->vdwtype = evdwCUT;
- ir->ndelta = 2;
ir->ns_type = ensGRID;
snew(ir->opts.egp_flags, 1);
#include <limits.h>
#include <assert.h>
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
#include "macros.h"
#include "readir.h"
#include "toputil.h"
#include "topio.h"
#include "gromacs/fileio/confio.h"
#include "readir.h"
-#include "symtab.h"
#include "names.h"
#include "grompp-impl.h"
-#include "gromacs/random/random.h"
#include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
-#include "vec.h"
-#include "gromacs/fileio/futil.h"
-#include "gromacs/commandline/pargs.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
#include "splitter.h"
#include "gromacs/gmxpreprocess/sortwater.h"
#include "convparm.h"
-#include "gmx_fatal.h"
#include "warninp.h"
-#include "index.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/trnio.h"
#include "gromacs/fileio/tpxio.h"
#include "perf_est.h"
#include "compute_io.h"
#include "gpp_atomtype.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "genborn.h"
#include "calc_verletbuf.h"
#include "tomorse.h"
#include "gromacs/imd/imd.h"
+ #include "gromacs/utility/cstringutil.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/random/random.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.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;
int i, j;
read_posres (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
- if (strcmp(fnA, fnB) != 0)
- {
- read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
- }
+ /* It is safer to simply read the b-state posres rather than trying
+ * to be smart and copy the positions.
+ */
+ read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
}
static void set_wall_atomtype(gpp_atomtype_t at, t_gromppopts *opts,
}
/* If we are using CMAP, setup the pre-interpolation grid */
- if (plist->ncmap > 0)
+ if (plist[F_CMAP].ncmap > 0)
{
- init_cmap_grid(&sys->ffparams.cmap_grid, plist->nc, plist->grid_spacing);
- setup_cmap(plist->grid_spacing, plist->nc, plist->cmap, &sys->ffparams.cmap_grid);
+ init_cmap_grid(&sys->ffparams.cmap_grid, plist[F_CMAP].nc, plist[F_CMAP].grid_spacing);
+ setup_cmap(plist[F_CMAP].grid_spacing, plist[F_CMAP].nc, plist[F_CMAP].cmap, &sys->ffparams.cmap_grid);
}
set_wall_atomtype(atype, opts, ir, wi);
check_vel(sys, state.v);
}
+ /* check for shells and inpurecs */
+ check_shells_inputrec(sys, ir, wi);
+
/* check masses */
check_mol(sys, wi);
#include <string.h>
#include "hackblock.h"
#include "gromacs/utility/smalloc.h"
-#include "vec.h"
-#include "macros.h"
+#include "gromacs/math/vec.h"
/* these MUST correspond to the enum in hackblock.h */
const char *btsNames[ebtsNR] = { "bonds", "angles", "dihedrals", "impropers", "exclusions", "cmap" };
free_t_bondeds(&(*rtp)[i].rb[j]);
}
}
- free(*rtp);
+ sfree(*rtp);
}
void free_t_hack(int nh, t_hack **h)
return bRet;
}
+ gmx_bool rbonded_atoms_exist_in_list(t_rbonded *b, t_rbonded blist[], int nlist, int natoms)
+ {
+ int i, k;
+ gmx_bool matchFound = FALSE;
+ gmx_bool atomsMatch;
+
+ for (i = 0; i < nlist && !matchFound; i++)
+ {
+ atomsMatch = TRUE;
+ for (k = 0; k < natoms && atomsMatch; k++)
+ {
+ atomsMatch = atomsMatch && !strcmp(b->a[k], blist[i].a[k]);
+ }
+ /* Try reverse if forward match did not work */
+ if (!atomsMatch)
+ {
+ atomsMatch = TRUE;
+ for (k = 0; k < natoms && atomsMatch; k++)
+ {
+ atomsMatch = atomsMatch && !strcmp(b->a[k], blist[i].a[natoms-1-k]);
+ }
+ }
+ matchFound = atomsMatch;
+ }
+ return matchFound;
+ }
+
gmx_bool merge_t_bondeds(t_rbondeds s[], t_rbondeds d[], gmx_bool bMin, gmx_bool bPlus)
{
int i, j;
srenew(d[i].b, d[i].nb + s[i].nb);
for (j = 0; j < s[i].nb; j++)
{
- if (!(bMin && contains_char(&s[i].b[j], '-'))
- && !(bPlus && contains_char(&s[i].b[j], '+')))
+ /* Check if this bonded string already exists before adding.
+ * We are merging from the main rtp to the hackblocks, so this
+ * will mean the hackblocks overwrite the man rtp, as intended.
+ */
+ if (!rbonded_atoms_exist_in_list(&s[i].b[j], d[i].b, d[i].nb, btsNiatoms[i]))
{
- copy_t_rbonded(&s[i].b[j], &d[i].b[ d[i].nb ]);
- d[i].nb++;
- }
- else if (i == ebtsBONDS)
- {
- bBondsRemoved = TRUE;
+ if (!(bMin && contains_char(&s[i].b[j], '-'))
+ && !(bPlus && contains_char(&s[i].b[j], '+')))
+ {
+ copy_t_rbonded(&s[i].b[j], &d[i].b[ d[i].nb ]);
+ d[i].nb++;
+ }
+ else if (i == ebtsBONDS)
+ {
+ bBondsRemoved = TRUE;
+ }
}
}
}
}
-
return bBondsRemoved;
}
#endif
#include <ctype.h>
+#include <stdlib.h>
#include <string.h>
#include <time.h>
-#include "sysstuff.h"
#include "typedefs.h"
#include "gromacs/fileio/gmxfio.h"
-#include "gromacs/utility/smalloc.h"
#include "copyrite.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/fileio/confio.h"
-#include "symtab.h"
-#include "vec.h"
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/fileio/futil.h"
-#include "gmx_fatal.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/fileio/pdbio.h"
#include "toputil.h"
#include "h_db.h"
-#include "physics.h"
#include "pgutil.h"
-#include "calch.h"
#include "resall.h"
#include "pdb2top.h"
#include "ter_db.h"
#include "gromacs/gmxlib/conformation-utilities.h"
#include "genhydro.h"
#include "readinp.h"
-#include "atomprop.h"
-#include "index.h"
+#include "gromacs/topology/index.h"
#include "fflibutil.h"
#include "macros.h"
+#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/strdb.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/residuetypes.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/fatalerror.h"
#include "hizzie.h"
#include "specbond.h"
static int read_pdball(const char *inf, const char *outf, char *title,
t_atoms *atoms, rvec **x,
int *ePBC, matrix box, gmx_bool bRemoveH,
- t_symtab *symtab, gmx_residuetype_t rt, const char *watres,
+ t_symtab *symtab, gmx_residuetype_t *rt, const char *watres,
gmx_atomprop_t aps, gmx_bool bVerbose)
/* Read a pdb file. (containing proteins) */
{
return pdba->nr;
}
-void find_nc_ter(t_atoms *pdba, int r0, int r1, int *r_start, int *r_end, gmx_residuetype_t rt)
+void find_nc_ter(t_atoms *pdba, int r0, int r1, int *r_start, int *r_end,
+ gmx_residuetype_t *rt)
{
int i;
const char *p_startrestype;
t_hackblock *ah;
t_symtab symtab;
gpp_atomtype_t atype;
- gmx_residuetype_t rt;
+ gmx_residuetype_t*rt;
const char *top_fn;
char fn[256], itp_fn[STRLEN], posre_fn[STRLEN], buf_fn[STRLEN];
char molname[STRLEN], title[STRLEN], quote[STRLEN];
printf("There are %d chains and %d blocks of water and "
"%d residues with %d atoms\n",
nch-nwaterchain, nwaterchain,
- pdba_all.resinfo[pdba_all.atom[natom-1].resind].nr, natom);
+ pdba_all.nres, natom);
printf("\n %5s %4s %6s\n", "chain", "#res", "#atoms");
for (i = 0; (i < nch); i++)
#include <math.h>
#include <ctype.h>
-#include "vec.h"
-#include "gromacs/utility/smalloc.h"
+#include "gromacs/math/vec.h"
#include "macros.h"
-#include "symtab.h"
-#include "gromacs/fileio/futil.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/futil.h"
#include "pdb2top.h"
#include "gpp_nextnb.h"
#include "topdirs.h"
#include "resall.h"
#include "topio.h"
#include "gromacs/utility/cstringutil.h"
-#include "physics.h"
#include "gromacs/fileio/pdbio.h"
#include "gen_ad.h"
#include "gromacs/fileio/filenm.h"
-#include "index.h"
#include "gen_vsite.h"
#include "add_par.h"
#include "toputil.h"
#include "copyrite.h"
#include "gromacs/fileio/strdb.h"
+#include "gromacs/topology/residuetypes.h"
+#include "gromacs/topology/symtab.h"
#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/programcontext.h"
+#include "gromacs/utility/smalloc.h"
/* this must correspond to enum in pdb2top.h */
const char *hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
sfree(model);
}
- static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
+ static int name2type(t_atoms *at, int **cgnr,
- t_restp restp[], gmx_residuetype_t rt)
+ t_restp restp[], gmx_residuetype_t *rt)
{
int i, j, prevresind, resind, i0, prevcg, cg, curcg;
char *name;
j = search_jtype(&restp[resind], name, bNterm);
at->atom[i].type = restp[resind].atom[j].type;
at->atom[i].q = restp[resind].atom[j].q;
- at->atom[i].m = get_atomtype_massA(restp[resind].atom[j].type,
- atype);
- cg = restp[resind].cgnr[j];
+ at->atom[i].m = restp[resind].atom[j].m;
+ cg = restp[resind].cgnr[j];
/* A charge group number -1 signals a separate charge group
* for this atom.
*/
}
#define NUM_CMAP_ATOMS 5
-static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t rt)
+static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t *rt)
{
int residx, i, j, k;
const char *ptr;
int *vsite_type;
int i, nmissat;
int bts[ebtsNR];
- gmx_residuetype_t rt;
+ gmx_residuetype_t*rt;
init_plist(plist);
gmx_residuetype_init(&rt);
atoms, nssbonds, ssbonds,
bAllowMissing);
- nmissat = name2type(atoms, &cgnr, atype, restp, rt);
+ nmissat = name2type(atoms, &cgnr, restp, rt);
if (nmissat)
{
if (bAllowMissing)
#endif
#include <ctype.h>
-#include <stdlib.h>
#include <limits.h>
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
+#include <stdlib.h>
+
#include "typedefs.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "names.h"
-#include "gmx_fatal.h"
#include "macros.h"
-#include "index.h"
-#include "symtab.h"
+#include "gromacs/topology/index.h"
#include "gromacs/utility/cstringutil.h"
#include "readinp.h"
#include "warninp.h"
#include "readir.h"
#include "toputil.h"
-#include "index.h"
#include "network.h"
-#include "vec.h"
-#include "pbc.h"
-#include "mtop_util.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/mtop_util.h"
#include "chargegroup.h"
#include "inputrec.h"
#include "calc_verletbuf.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+
#define MAXPTR 254
#define NOGID 255
}
}
+ if (ir->nsteps == 0 && !ir->bContinuation)
+ {
+ warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
+ }
+
/* LD STUFF */
if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
ir->bContinuation && ir->ld_seed != -1)
}
}
+ if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
+ {
+ sprintf(err_buf,
+ "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
+ CHECK( ir->coulomb_modifier != eintmodNONE);
+ }
+ if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
+ {
+ sprintf(err_buf,
+ "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
+ CHECK( ir->vdw_modifier != eintmodNONE);
+ }
+
if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
{
if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
{
real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
- sprintf(warn_buf, "The switching range for should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
+ sprintf(warn_buf, "The switching range should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
warning(wi, warn_buf);
}
warning_note(wi, warn_buf);
}
- /* remove the following deprecated commands */
+ /* ignore the following deprecated commands */
REM_TYPE("title");
REM_TYPE("cpp");
REM_TYPE("domain-decomposition");
REM_TYPE("dihre-tau");
REM_TYPE("nstdihreout");
REM_TYPE("nstcheckpoint");
+ REM_TYPE("optimize-fft");
/* replace the following commands with the clearer new versions*/
REPL_TYPE("unconstrained-start", "continuation");
ITYPE ("nstxout", ir->nstxout, 0);
ITYPE ("nstvout", ir->nstvout, 0);
ITYPE ("nstfout", ir->nstfout, 0);
- ir->nstcheckpoint = 1000;
CTYPE ("Output frequency for energies to log file and energy file");
ITYPE ("nstlog", ir->nstlog, 1000);
ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
ITYPE ("nstlist", ir->nstlist, 10);
CTYPE ("ns algorithm (simple or grid)");
EETYPE("ns-type", ir->ns_type, ens_names);
- /* set ndelta to the optimal value of 2 */
- ir->ndelta = 2;
CTYPE ("Periodic boundary conditions: xyz, no, xy");
EETYPE("pbc", ir->ePBC, epbc_names);
EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
- EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
CCTYPE("IMPLICIT SOLVENT ALGORITHM");
EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
}
else
{
- warning(wi, "Can not couple a molecule with free_energy = no");
+ warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
}
}
/* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
{
- int ig = -1, i;
+ int ig, i;
ig = search_string(IMDgname, grps->nr, gnames);
#include <stdlib.h>
#include "gromacs/utility/cstringutil.h"
-#include "sysstuff.h"
-#include "princ.h"
-#include "gromacs/fileio/futil.h"
-#include "vec.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/math/vec.h"
#include "gromacs/utility/smalloc.h"
#include "typedefs.h"
#include "names.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
#include "macros.h"
-#include "index.h"
-#include "symtab.h"
#include "readinp.h"
#include "readir.h"
#include "mdatoms.h"
-#include "pbc.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
RTYPE("pull-r1", pull->cyl_r1, 1.0);
CTYPE("Switch from r1 to r0 in case of dynamic reaction force");
RTYPE("pull-r0", pull->cyl_r0, 1.5);
- RTYPE("pull_constr_tol", pull->constr_tol, 1E-6);
+ RTYPE("pull-constr-tol", pull->constr_tol, 1E-6);
EETYPE("pull-start", *bStart, yesno_names);
EETYPE("pull-print-reference", pull->bPrintRef, yesno_names);
ITYPE("pull-nstxout", pull->nstxout, 10);
{
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];
#include <config.h>
#endif
+#include <assert.h>
+#include <ctype.h>
+#include <errno.h>
#include <math.h>
-#include <sys/types.h>
#include <stdio.h>
+#include <stdlib.h>
#include <string.h>
-#include <errno.h>
-#include <ctype.h>
-#include <assert.h>
-#include "gromacs/fileio/futil.h"
-#include "sysstuff.h"
+#include <sys/types.h>
+
+#include "gromacs/utility/futil.h"
#include "typedefs.h"
#include "gromacs/utility/smalloc.h"
#include "macros.h"
#include "gromacs/fileio/gmxfio.h"
#include "txtdump.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "macros.h"
#include "names.h"
#include "gromacs/utility/cstringutil.h"
-#include "symtab.h"
-#include "gmx_fatal.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/fatalerror.h"
#include "warninp.h"
#include "vsite_parm.h"
comb = 0;
/* Init the number of CMAP torsion angles and grid spacing */
- plist->grid_spacing = 0;
- plist->nc = 0;
+ plist[F_CMAP].grid_spacing = 0;
+ plist[F_CMAP].nc = 0;
bWarn_copy_A_B = bFEP;
j = 0;
while (j < molt->ilist[i].nr)
{
- bexcl = FALSE;
a1 = molt->ilist[i].iatoms[j+1];
a2 = molt->ilist[i].iatoms[j+2];
bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
#include <config.h>
#endif
+#include <assert.h>
#include <ctype.h>
#include <math.h>
-#include <assert.h>
+#include <stdlib.h>
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "gromacs/utility/cstringutil.h"
#include "names.h"
#include "toputil.h"
#include "toppush.h"
#include "topdirs.h"
#include "readir.h"
-#include "symtab.h"
-#include "gmx_fatal.h"
#include "warninp.h"
#include "gpp_atomtype.h"
#include "gpp_bond_atomtype.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+
void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype,
warninp_t wi)
{
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);
}
}
nrfp = nrfpA+nrfpB;
/* Allocate memory for the CMAP grid */
- bt->ncmap += nrfp;
- srenew(bt->cmap, bt->ncmap);
+ bt[F_CMAP].ncmap += nrfp;
+ srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap);
/* Read in CMAP parameters */
sl = 0;
}
nn = sscanf(line+start+sl, " %s ", s);
sl += strlen(s);
- bt->cmap[i+(bt->ncmap)-nrfp] = strtod(s, NULL);
+ bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, NULL);
if (nn == 1)
{
{
for (i = 0; i < ncmap; i++)
{
- bt->cmap[i+ncmap] = bt->cmap[i];
+ bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i];
}
}
else
/* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
* so we can safely assign them each time
*/
- bt->grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
- bt->nc = bt->nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
- nct = (nral+1) * bt->nc;
+ bt[F_CMAP].grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
+ bt[F_CMAP].nc = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
+ nct = (nral+1) * bt[F_CMAP].nc;
/* Allocate memory for the cmap_types information */
- srenew(bt->cmap_types, nct);
+ srenew(bt[F_CMAP].cmap_types, nct);
for (i = 0; (i < nral); i++)
{
}
/* Assign a grid number to each cmap_type */
- bt->cmap_types[bt->nct++] = get_bond_atomtype_type(alc[i], bat);
+ bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat);
}
/* Assign a type number to this cmap */
- bt->cmap_types[bt->nct++] = bt->nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
+ bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = bt[F_CMAP].nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
/* Check for the correct number of atoms (again) */
- if (bt->nct != nct)
+ if (bt[F_CMAP].nct != nct)
{
- gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt->nc);
+ gmx_fatal(FARGS, "Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
}
/* Is this correct?? */
ct = 0;
/* Match the current cmap angle against the list of cmap_types */
- for (i = 0; i < bondtype->nct && !bFound; i += 6)
+ for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6)
{
if (bB)
{
else
{
if (
- (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype->cmap_types[i]) &&
- (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype->cmap_types[i+1]) &&
- (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype->cmap_types[i+2]) &&
- (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype->cmap_types[i+3]) &&
- (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype->cmap_types[i+4]))
+ (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i]) &&
+ (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) &&
+ (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) &&
+ (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) &&
+ (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[i+4]))
{
/* Found cmap torsion */
bFound = TRUE;
- ct = bondtype->cmap_types[i+5];
+ ct = bondtype[F_CMAP].cmap_types[i+5];
nparam_found = 1;
}
}
#include <ctype.h>
#include <string.h>
+
#include "typedefs.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/smalloc.h"
-#include "symtab.h"
-#include "index.h"
-#include "gromacs/fileio/futil.h"
#include "fflibutil.h"
#include "hackblock.h"
-#include "gmx_fatal.h"
#include "xlate.h"
#include "gromacs/fileio/strdb.h"
+#include "gromacs/topology/residuetypes.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
typedef struct {
char *filebase;
void rename_atoms(const char *xlfile, const char *ffdir,
t_atoms *atoms, t_symtab *symtab, const t_restp *restp,
- gmx_bool bResname, gmx_residuetype_t rt, gmx_bool bReorderNum,
+ gmx_bool bResname, gmx_residuetype_t *rt, gmx_bool bReorderNum,
gmx_bool bVerbose)
{
FILE *fp;
char **f;
char c, *rnm, atombuf[32], *ptr0, *ptr1;
gmx_bool bReorderedNum, bRenamed, bMatch;
+ gmx_bool bStartTerm, bEndTerm;
nxlate = 0;
xlatom = NULL;
for (a = 0; (a < atoms->nr); a++)
{
resind = atoms->atom[a].resind;
+
+ bStartTerm = (resind == 0) || atoms->resinfo[resind].chainnum != atoms->resinfo[resind-1].chainnum;
+ bEndTerm = (resind >= atoms->nres-1) || atoms->resinfo[resind].chainnum != atoms->resinfo[resind+1].chainnum;
+
if (bResname)
{
rnm = *(atoms->resinfo[resind].name);
{
/* Match the residue name */
bMatch = (xlatom[i].res == NULL ||
+ (gmx_strcasecmp("protein-nterm", xlatom[i].res) == 0 &&
+ gmx_residuetype_is_protein(rt, rnm) && bStartTerm) ||
+ (gmx_strcasecmp("protein-cterm", xlatom[i].res) == 0 &&
+ gmx_residuetype_is_protein(rt, rnm) && bEndTerm) ||
(gmx_strcasecmp("protein", xlatom[i].res) == 0 &&
gmx_residuetype_is_protein(rt, rnm)) ||
(gmx_strcasecmp("DNA", xlatom[i].res) == 0 &&
#define _nonbonded_h
#include "typedefs.h"
-#include "pbc.h"
#include "network.h"
#include "tgroup.h"
#include "genborn.h"
} /* fixes auto-indentation problems */
#endif
+struct t_graph;
+struct t_pbc;
void
void
gmx_nonbonded_set_kernel_pointers(FILE * fplog,
- t_nblist * nl);
+ t_nblist * nl,
+ gmx_bool bElecAndVdwSwitchDiffers);
*/
real
do_nonbonded_listed(int ftype, int nbonds, const t_iatom iatoms[], const t_iparams iparams[],
- const rvec x[], rvec f[], rvec fshift[], const t_pbc *pbc, const t_graph *g,
+ const rvec x[], rvec f[], rvec fshift[],
+ const struct t_pbc *pbc, const struct t_graph *g,
real *lambda, real *dvdl, const t_mdatoms *md, const t_forcerec *fr,
gmx_grppairener_t *grppener, int *global_atom_index);
* the research papers on the package. Check out http://www.gromacs.org.
*/
+#include "../timing/wallcycle.h"
+
#include "typedefs.h"
#include "vsite.h"
extern "C" {
#endif
+struct t_graph;
+
/* Initialization function, also predicts the initial shell postions.
* 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);
tensor force_vir,
t_mdatoms *md,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- t_graph *graph,
+ struct t_graph *graph,
gmx_groups_t *groups,
gmx_shellfc_t shfc,
t_forcerec *fr,
#ifndef _inputrec_h_
#define _inputrec_h_
+#include <stdio.h>
#include "simple.h"
#include "enums.h"
-#include "../sysstuff.h"
#include "../../swap/enums.h"
#ifdef __cplusplus
int nstlist; /* number of steps before pairlist is generated */
int ndelta; /* number of cells per rlong */
int nstcomm; /* number of steps after which center of mass */
- /* motion is removed */
+ /* motion is removed */
int comm_mode; /* Center of mass motion removal algorithm */
- int nstcheckpoint; /* checkpointing frequency */
int nstlog; /* number of steps after which print to logfile */
int nstxout; /* number of steps after which X is output */
int nstvout; /* id. for V */
real ewald_rtol_lj; /* Real space tolerance for LJ-Ewald */
int ewald_geometry; /* normal/3d ewald, or pseudo-2d LR corrections */
real epsilon_surface; /* Epsilon for PME dipole correction */
- gmx_bool bOptFFT; /* optimize the fft plan at start */
int ljpme_combination_rule; /* Type of combination rule in LJ-PME */
int ePBC; /* Type of periodic boundary conditions */
int bPeriodicMols; /* Periodic molecules */
real orires_fc; /* force constant for orientational restraints */
real orires_tau; /* time constant for memory function in orires */
int nstorireout; /* frequency of writing tr(SD) to enx */
- real dihre_fc; /* force constant for dihedral restraints (obsolete) */
real em_stepsize; /* The stepsize for updating */
real em_tol; /* The tolerance */
int niter; /* Number of iterations for convergence of */
#include <config.h>
#endif
+#include <stdlib.h>
+
#include "gromacs/fileio/confio.h"
#include "types/commrec.h"
#include "constr.h"
#include "copyrite.h"
-#include "invblock.h"
-#include "main.h"
#include "mdrun.h"
#include "nrnb.h"
-#include "gromacs/utility/smalloc.h"
-#include "vec.h"
-#include "physics.h"
+#include "gromacs/math/vec.h"
#include "names.h"
#include "txtdump.h"
#include "domdec.h"
#include "gromacs/fileio/pdbio.h"
#include "splitter.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "gromacs/fileio/gmxfio.h"
#include "macros.h"
#include "gmx_omp_nthreads.h"
#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/pulling/pull.h"
-#include "gmx_fatal.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/invblock.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
typedef struct gmx_constr {
int ncon_tot; /* The total number of constraints */
int start, int homenr, t_commrec *cr,
rvec x[], matrix box)
{
- char fname[STRLEN], format[STRLEN];
+ char fname[STRLEN];
FILE *out;
int dd_ac0 = 0, dd_ac1 = 0, i, ii, resnr;
gmx_domdec_t *dd;
{
sprintf(fname, "%s.pdb", fn);
}
- sprintf(format, "%s\n", get_pdbformat());
out = gmx_fio_fopen(fname, "w");
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]);
+ gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, anm, ' ', resnm, ' ', resnr, ' ',
+ 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 0.0, "");
}
fprintf(out, "TER\n");
j, constr->nblocks, ncons);
for (i = 0; (i < ncons); i++)
{
- fprintf(stderr, "i: %5d sb[i].blocknr: %5u\n", i, sb[i].blocknr);
+ fprintf(stderr, "i: %5d sb[i].blocknr: %5d\n", i, sb[i].blocknr);
}
for (j = 0; (j <= constr->nblocks); j++)
{
#include <assert.h>
#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "gmx_fatal.h"
-#include "gmx_fatal_collective.h"
-#include "vec.h"
+#include "network.h"
+#include "gromacs/math/vec.h"
#include "domdec.h"
#include "domdec_network.h"
#include "nrnb.h"
-#include "pbc.h"
#include "chargegroup.h"
#include "constr.h"
#include "mdatoms.h"
#include "mdrun.h"
#include "nsgrid.h"
#include "shellfc.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "gmx_ga2la.h"
#include "macros.h"
#include "nbnxn_search.h"
#include "gmx_omp_nthreads.h"
#include "gpu_utils.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/pdbio.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/pulling/pull_rotation.h"
+#include "gromacs/swap/swapcoords.h"
#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/basenetwork.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxmpi.h"
-#include "gromacs/swap/swapcoords.h"
#include "gromacs/utility/qsort_threadsafe.h"
-#include "gromacs/pulling/pull.h"
-#include "gromacs/pulling/pull_rotation.h"
-#include "gromacs/imd/imd.h"
+#include "gromacs/utility/smalloc.h"
#define DDRANK(dd, rank) (rank)
#define DDMASTERRANK(dd) (dd->masterrank)
gmx_domdec_t *dd, matrix box, gmx_ddbox_t *ddbox)
{
rvec grid_s[2], *grid_r = NULL, cx, r;
- char fname[STRLEN], format[STRLEN], buf[22];
+ char fname[STRLEN], buf[22];
FILE *out;
int a, i, d, z, y, x;
matrix tric;
snew(grid_r, 2*dd->nnodes);
}
- dd_gather(dd, 2*sizeof(rvec), grid_s[0], DDMASTER(dd) ? grid_r[0] : NULL);
+ dd_gather(dd, 2*sizeof(rvec), grid_s, DDMASTER(dd) ? grid_r : NULL);
if (DDMASTER(dd))
{
}
}
sprintf(fname, "%s_%s.pdb", fn, gmx_step_str(step, buf));
- sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
out = gmx_fio_fopen(fname, "w");
gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
a = 1;
cx[YY] = grid_r[i*2+y][YY];
cx[ZZ] = grid_r[i*2+z][ZZ];
mvmul(tric, cx, r);
- fprintf(out, format, "ATOM", a++, "CA", "GLY", ' ', 1+i,
- ' ', 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol);
+ gmx_fprintf_pdb_atomline(out, epdbATOM, a++, "CA", ' ', "GLY", ' ', i+1, ' ',
+ 10*r[XX], 10*r[YY], 10*r[ZZ], 1.0, vol, "");
}
}
}
gmx_mtop_t *mtop, t_commrec *cr,
int natoms, rvec x[], matrix box)
{
- char fname[STRLEN], format[STRLEN], format4[STRLEN], buf[22];
+ char fname[STRLEN], buf[22];
FILE *out;
int i, ii, resnr, c;
char *atomname, *resname;
sprintf(fname, "%s_%s_n%d.pdb", fn, gmx_step_str(step, buf), cr->sim_nodeid);
- sprintf(format, "%s%s\n", get_pdbformat(), "%6.2f%6.2f");
- sprintf(format4, "%s%s\n", get_pdbformat4(), "%6.2f%6.2f");
-
out = gmx_fio_fopen(fname, "w");
fprintf(out, "TITLE %s\n", title);
{
b = dd->comm->zones.n + 1;
}
- fprintf(out, strlen(atomname) < 4 ? format : format4,
- "ATOM", (ii+1)%100000,
- atomname, resname, ' ', resnr%10000, ' ',
- 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b);
+ gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, atomname, ' ', resname, ' ', resnr, ' ',
+ 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, b, "");
}
fprintf(out, "TER\n");
#include <config.h>
#endif
+#include <assert.h>
#include <math.h>
#include <string.h>
-#include <assert.h>
-#include "sysstuff.h"
+
#include "typedefs.h"
#include "types/commrec.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "gromacs/math/utilities.h"
#include "macros.h"
-#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "gmx_fatal.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "force.h"
#include "tables.h"
#include "nonbonded.h"
-#include "invblock.h"
#include "names.h"
#include "network.h"
-#include "pbc.h"
#include "ns.h"
-#include "mshift.h"
#include "txtdump.h"
#include "coulomb.h"
#include "md_support.h"
#include "domdec.h"
#include "qmmm.h"
#include "copyrite.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "nbnxn_simd.h"
#include "nbnxn_search.h"
#include "nbnxn_atomdata.h"
#include "gmx_detect_hardware.h"
#include "inputrec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+
#include "types/nbnxn_cuda_types_ext.h"
#include "gpu_utils.h"
#include "nbnxn_cuda_data_mgmt.h"
sigmaj = pow(c12j / c6j, 1.0/6.0);
epsi = c6i * c6i / c12i;
epsj = c6j * c6j / c12j;
- c6 = epsi * epsj * pow(0.5*(sigmai+sigmaj), 6);
- c12 = epsi * epsj * pow(0.5*(sigmai+sigmaj), 12);
+ c6 = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 6);
+ c12 = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 12);
}
C6(nbfp, atnr, i, j) = c6*6.0;
C12(nbfp, atnr, i, j) = c12*12.0;
fr->nbkernel_elec_modifier = fr->coulomb_modifier;
fr->nbkernel_vdw_modifier = fr->vdw_modifier;
+ fr->rvdw = cutoff_inf(ir->rvdw);
+ fr->rvdw_switch = ir->rvdw_switch;
+ fr->rcoulomb = cutoff_inf(ir->rcoulomb);
+ fr->rcoulomb_switch = ir->rcoulomb_switch;
+
fr->bTwinRange = fr->rlistlong > fr->rlist;
fr->bEwald = (EEL_PME(fr->eeltype) || fr->eeltype == eelEWALD);
/* If the user absolutely wants different switch/shift settings for coul/vdw, it is likely
* going to be faster to tabulate the interaction than calling the generic kernel.
+ * However, if generic kernels have been requested we keep things analytically.
*/
- if (fr->nbkernel_elec_modifier == eintmodPOTSWITCH && fr->nbkernel_vdw_modifier == eintmodPOTSWITCH)
+ if (fr->nbkernel_elec_modifier == eintmodPOTSWITCH &&
+ fr->nbkernel_vdw_modifier == eintmodPOTSWITCH &&
+ bGenericKernelOnly == FALSE)
{
if ((fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw))
{
fr->bcoultab = TRUE;
+ /* Once we tabulate electrostatics, we can use the switch function for LJ,
+ * which would otherwise need two tables.
+ */
}
}
else if ((fr->nbkernel_elec_modifier == eintmodPOTSHIFT && fr->nbkernel_vdw_modifier == eintmodPOTSHIFT) ||
fr->nbkernel_elec_modifier == eintmodEXACTCUTOFF &&
(fr->nbkernel_vdw_modifier == eintmodPOTSWITCH || fr->nbkernel_vdw_modifier == eintmodPOTSHIFT))))
{
- if (fr->rcoulomb != fr->rvdw)
+ if ((fr->rcoulomb != fr->rvdw) && (bGenericKernelOnly == FALSE))
{
fr->bcoultab = TRUE;
}
}
+ if (fr->nbkernel_elec_modifier == eintmodFORCESWITCH)
+ {
+ fr->bcoultab = TRUE;
+ }
+ if (fr->nbkernel_vdw_modifier == eintmodFORCESWITCH)
+ {
+ fr->bvdwtab = TRUE;
+ }
+
if (getenv("GMX_REQUIRE_TABLES"))
{
fr->bvdwtab = TRUE;
fr->epsilon_r = ir->epsilon_r;
fr->epsilon_rf = ir->epsilon_rf;
fr->fudgeQQ = mtop->ffparams.fudgeQQ;
- fr->rcoulomb_switch = ir->rcoulomb_switch;
- fr->rcoulomb = cutoff_inf(ir->rcoulomb);
/* Parameters for generalized RF */
fr->zsquare = 0.0;
fr->egp_flags = ir->opts.egp_flags;
/* Van der Waals stuff */
- fr->rvdw = cutoff_inf(ir->rvdw);
- fr->rvdw_switch = ir->rvdw_switch;
if ((fr->vdwtype != evdwCUT) && (fr->vdwtype != evdwUSER) && !fr->bBHAM)
{
if (fr->rvdw_switch >= fr->rvdw)
(ir->eDispCorr != edispcNO && ir_vdw_switched(ir));
bMakeSeparate14Table = ((!bMakeTables || fr->eeltype != eelCUT || fr->vdwtype != evdwCUT ||
+ fr->coulomb_modifier != eintmodNONE ||
+ fr->vdw_modifier != eintmodNONE ||
fr->bBHAM || fr->bEwald) &&
(gmx_mtop_ftype_count(mtop, F_LJ14) > 0 ||
gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0 ||
}
}
}
+ else if ((fr->eDispCorr != edispcNO) &&
+ ((fr->vdw_modifier == eintmodPOTSWITCH) ||
+ (fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdw_modifier == eintmodPOTSHIFT)))
+ {
+ /* Tables might not be used for the potential modifier interactions per se, but
+ * we still need them to evaluate switch/shift dispersion corrections in this case.
+ */
+ make_nbf_tables(fp, oenv, fr, rtab, cr, tabfn, NULL, NULL, &fr->nblists[0]);
+ }
+
if (bMakeSeparate14Table)
{
/* generate extra tables with plain Coulomb for 1-4 interactions only */
#endif
#include <math.h>
+#include <stdlib.h>
#include <string.h>
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
+
#include "macros.h"
#include "gromacs/math/utilities.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "types/commrec.h"
#include "network.h"
#include "nsgrid.h"
#include "force.h"
#include "nonbonded.h"
#include "ns.h"
-#include "pbc.h"
#include "names.h"
-#include "gmx_fatal.h"
#include "nrnb.h"
#include "txtdump.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "domdec.h"
#include "adress.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
/*
* E X C L U S I O N H A N D L I N G
int maxsr, int maxlr,
int ivdw, int ivdwmod,
int ielec, int ielecmod,
- int igeometry, int type)
+ int igeometry, int type,
+ gmx_bool bElecAndVdwSwitchDiffers)
{
t_nblist *nl;
int homenr;
}
/* This will also set the simd_padding_width field */
- gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl);
+ gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl, bElecAndVdwSwitchDiffers);
/* maxnri is influenced by the number of shifts (maximum is 8)
* and the number of energy groups.
* cache trashing.
*/
int maxsr, maxsr_wat, maxlr, maxlr_wat;
- int ielec, ielecf, ivdw, ielecmod, ielecmodf, ivdwmod, type;
+ int ielec, ivdw, ielecmod, ivdwmod, type;
int solvent;
int igeometry_def, igeometry_w, igeometry_ww;
int i;
+ gmx_bool bElecAndVdwSwitchDiffers;
t_nblists *nbl;
/* maxsr = homenr-fr->nWatMol*3; */
}
/* Determine the values for ielec/ivdw. */
- ielec = fr->nbkernel_elec_interaction;
- ivdw = fr->nbkernel_vdw_interaction;
- ielecmod = fr->nbkernel_elec_modifier;
- ivdwmod = fr->nbkernel_vdw_modifier;
- type = GMX_NBLIST_INTERACTION_STANDARD;
+ ielec = fr->nbkernel_elec_interaction;
+ ivdw = fr->nbkernel_vdw_interaction;
+ ielecmod = fr->nbkernel_elec_modifier;
+ ivdwmod = fr->nbkernel_vdw_modifier;
+ type = GMX_NBLIST_INTERACTION_STANDARD;
+ bElecAndVdwSwitchDiffers = ( (fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw));
fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
if (!fr->ns.bCGlist)
type = GMX_NBLIST_INTERACTION_ADRESS;
}
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ], &nbl->nlist_lr[eNL_VDWQQ],
- maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type);
+ maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDW], &nbl->nlist_lr[eNL_VDW],
- maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type);
+ maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ], &nbl->nlist_lr[eNL_QQ],
- maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type);
+ maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATER], &nbl->nlist_lr[eNL_VDWQQ_WATER],
- maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type);
+ maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATER], &nbl->nlist_lr[eNL_QQ_WATER],
- maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type);
+ maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATERWATER], &nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
- maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type);
+ maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ_WATERWATER], &nbl->nlist_lr[eNL_QQ_WATERWATER],
- maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type);
+ maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
/* Did we get the solvent loops so we can use optimized water kernels? */
if (nbl->nlist_sr[eNL_VDWQQ_WATER].kernelptr_vf == NULL
if (fr->efep != efepNO)
{
- if ((fr->bEwald) && (fr->sc_alphacoul > 0)) /* need to handle long range differently if using softcore */
- {
- ielecf = GMX_NBKERNEL_ELEC_EWALD;
- ielecmodf = eintmodNONE;
- }
- else
- {
- ielecf = ielec;
- ielecmodf = ielecmod;
- }
-
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_FREE], &nbl->nlist_lr[eNL_VDWQQ_FREE],
- maxsr, maxlr, ivdw, ivdwmod, ielecf, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+ maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_VDW_FREE], &nbl->nlist_lr[eNL_VDW_FREE],
- maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+ maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
init_nblist(log, &nbl->nlist_sr[eNL_QQ_FREE], &nbl->nlist_lr[eNL_QQ_FREE],
- maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielecf, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+ maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
}
}
/* QMMM MM list */
if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
{
init_nblist(log, &fr->QMMMlist, NULL,
- maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD);
+ maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD, bElecAndVdwSwitchDiffers);
}
if (log != NULL)
gmx_ns_t *ns;
atom_id **nl_lr_ljc, **nl_lr_one, **nl_sr;
int *nlr_ljc, *nlr_one, *nsr;
- gmx_domdec_t *dd = NULL;
+ gmx_domdec_t *dd;
t_block *cgs = &(top->cgs);
int *cginfo = fr->cginfo;
/* atom_id *i_atoms,*cgsindex=cgs->index; */
ns = &fr->ns;
bDomDec = DOMAINDECOMP(cr);
- if (bDomDec)
- {
- dd = cr->dd;
- }
+ dd = cr->dd;
bTriclinicX = ((YY < grid->npbcdim &&
(!bDomDec || dd->nc[YY] == 1) && box[YY][XX] != 0) ||
#include <config.h>
#endif
+#include <assert.h>
+#include <math.h>
#include <stdio.h>
+#include <stdlib.h>
#include <string.h>
-#include <math.h>
-#include <assert.h>
+
#include "typedefs.h"
#include "txtdump.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "gromacs/utility/smalloc.h"
#include "coulomb.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
#include "pme.h"
#include "network.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "nrnb.h"
#include "macros.h"
+#include "gromacs/legacyheaders/types/commrec.h"
#include "gromacs/fft/parallel_3dfft.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/fileio/pdbio.h"
#include "gromacs/math/gmxcomplex.h"
#include "gromacs/timing/cyclecounter.h"
{
#ifdef DEBUG_PME
FILE *fp, *fp2;
- char fn[STRLEN], format[STRLEN];
+ char fn[STRLEN];
real val;
sprintf(fn, "pmegrid%d.pdb", pme->nodeid);
fp = gmx_ffopen(fn, "w");
sprintf(fn, "pmegrid%d.txt", pme->nodeid);
fp2 = gmx_ffopen(fn, "w");
- sprintf(format, "%s%s\n", pdbformat, "%6.2f%6.2f");
#endif
for (ix = 0; ix < local_fft_ndata[XX]; ix++)
val = 100*pmegrid[pmeidx];
if (pmegrid[pmeidx] != 0)
{
- fprintf(fp, format, "ATOM", pmeidx, "CA", "GLY", ' ', pmeidx, ' ',
- 5.0*ix, 5.0*iy, 5.0*iz, 1.0, val);
+ gmx_fprintf_pdb_atomline(fp, epdbATOM, pmeidx, "CA", ' ', "GLY", ' ', pmeidx, ' ',
+ 5.0*ix, 5.0*iy, 5.0*iz, 1.0, val, "");
}
if (pmegrid[pmeidx] != 0)
{
/* 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--)
{
#include <config.h>
#endif
+#include <stdlib.h>
#include <string.h>
+
#include "typedefs.h"
#include "types/commrec.h"
-#include "gromacs/utility/smalloc.h"
-#include "gmx_fatal.h"
-#include "vec.h"
#include "txtdump.h"
#include "force.h"
#include "mdrun.h"
#include "names.h"
#include "constr.h"
#include "domdec.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "shellfc.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "chargegroup.h"
#include "macros.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/mshift.h"
++#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
typedef struct {
int nnucl;
}
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);
#include "typedefs.h"
#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/smalloc.h"
#include "names.h"
#include "txtdump.h"
-#include "pbc.h"
+#include "gromacs/pbcutil/pbc.h"
#include "chargegroup.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "nrnb.h"
-#include "mshift.h"
#include "mdrun.h"
#include "sim_util.h"
#include "update.h"
-#include "physics.h"
-#include "main.h"
+#include "gromacs/math/units.h"
#include "mdatoms.h"
#include "force.h"
#include "bondf.h"
#include "network.h"
#include "calcmu.h"
#include "constr.h"
-#include "xvgr.h"
#include "copyrite.h"
#include "domdec.h"
#include "genborn.h"
#include "../gmxlib/nonbonded/nb_kernel.h"
#include "../gmxlib/nonbonded/nb_free_energy.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/mshift.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/timing/walltime_accounting.h"
#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/utility/smalloc.h"
#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/pulling/pull_rotation.h"
void calc_enervirdiff(FILE *fplog, int eDispCorr, t_forcerec *fr)
{
- double eners[2], virs[2], enersum, virsum, y0, f, g, h;
- double r0, r1, r, rc3, rc9, ea, eb, ec, pa, pb, pc, pd;
- double invscale, invscale2, invscale3;
- int ri0, ri1, ri, i, offstart, offset;
- real scale, *vdwtab, tabfactor, tmp;
+ double eners[2], virs[2], enersum, virsum, y0, f, g, h;
+ double r0, r1, r, rc3, rc9, ea, eb, ec, pa, pb, pc, pd;
+ double invscale, invscale2, invscale3;
+ int ri0, ri1, ri, i, offstart, offset;
+ real scale, *vdwtab, tabfactor, tmp;
fr->enershiftsix = 0;
fr->enershifttwelve = 0;
eners[i] = 0;
virs[i] = 0;
}
- if (fr->vdwtype == evdwSWITCH || fr->vdwtype == evdwSHIFT ||
- fr->vdw_modifier == eintmodPOTSWITCH ||
- fr->vdw_modifier == eintmodFORCESWITCH)
+ if ((fr->vdw_modifier == eintmodPOTSHIFT) ||
+ (fr->vdw_modifier == eintmodPOTSWITCH) ||
+ (fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdwtype == evdwSHIFT) ||
+ (fr->vdwtype == evdwSWITCH))
{
- if (fr->rvdw_switch == 0)
+ if (((fr->vdw_modifier == eintmodPOTSWITCH) ||
+ (fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdwtype == evdwSWITCH)) && fr->rvdw_switch == 0)
{
gmx_fatal(FARGS,
"With dispersion correction rvdw-switch can not be zero "
"for vdw-type = %s", evdw_names[fr->vdwtype]);
}
- scale = fr->nblists[0].table_elec_vdw.scale;
+ scale = fr->nblists[0].table_vdw.scale;
vdwtab = fr->nblists[0].table_vdw.data;
/* Round the cut-offs to exact table values for precision */
ri0 = floor(fr->rvdw_switch*scale);
ri1 = ceil(fr->rvdw*scale);
+
+ /* The code below has some support for handling force-switching, i.e.
+ * when the force (instead of potential) is switched over a limited
+ * region. This leads to a constant shift in the potential inside the
+ * switching region, which we can handle by adding a constant energy
+ * term in the force-switch case just like when we do potential-shift.
+ *
+ * For now this is not enabled, but to keep the functionality in the
+ * code we check separately for switch and shift. When we do force-switch
+ * the shifting point is rvdw_switch, while it is the cutoff when we
+ * have a classical potential-shift.
+ *
+ * For a pure potential-shift the potential has a constant shift
+ * all the way out to the cutoff, and that is it. For other forms
+ * we need to calculate the constant shift up to the point where we
+ * start modifying the potential.
+ */
+ ri0 = (fr->vdw_modifier == eintmodPOTSHIFT) ? ri1 : ri0;
+
r0 = ri0/scale;
r1 = ri1/scale;
rc3 = r0*r0*r0;
rc9 = rc3*rc3*rc3;
- if (fr->vdwtype == evdwSHIFT ||
- fr->vdw_modifier == eintmodFORCESWITCH)
+ if ((fr->vdw_modifier == eintmodFORCESWITCH) ||
+ (fr->vdwtype == evdwSHIFT))
{
/* Determine the constant energy shift below rvdw_switch.
* Table has a scale factor since we have scaled it down to compensate
fr->enershiftsix = (real)(-1.0/(rc3*rc3)) - 6.0*vdwtab[8*ri0];
fr->enershifttwelve = (real)( 1.0/(rc9*rc3)) - 12.0*vdwtab[8*ri0 + 4];
}
+ else if (fr->vdw_modifier == eintmodPOTSHIFT)
+ {
+ fr->enershiftsix = (real)(-1.0/(rc3*rc3));
+ fr->enershifttwelve = (real)( 1.0/(rc9*rc3));
+ }
+
/* Add the constant part from 0 to rvdw_switch.
* This integration from 0 to rvdw_switch overcounts the number
* of interactions by 1, as it also counts the self interaction.
*/
eners[0] += 4.0*M_PI*fr->enershiftsix*rc3/3.0;
eners[1] += 4.0*M_PI*fr->enershifttwelve*rc3/3.0;
+
+ /* Calculate the contribution in the range [r0,r1] where we
+ * modify the potential. For a pure potential-shift modifier we will
+ * have ri0==ri1, and there will not be any contribution here.
+ */
for (i = 0; i < 2; i++)
{
enersum = 0;
virs[i] -= virsum;
}
- /* now add the correction for rvdw_switch to infinity */
+ /* Alright: Above we compensated by REMOVING the parts outside r0
+ * corresponding to the ideal VdW 1/r6 and /r12 potentials.
+ *
+ * Regardless of whether r0 is the point where we start switching,
+ * or the cutoff where we calculated the constant shift, we include
+ * all the parts we are missing out to infinity from r0 by
+ * calculating the analytical dispersion correction.
+ */
eners[0] += -4.0*M_PI/(3.0*rc3);
eners[1] += 4.0*M_PI/(9.0*rc9);
virs[0] += 8.0*M_PI/rc3;
evdw_names[fr->vdwtype]);
}
- /* TODO: remove this code once we have group LJ-PME kernels
- * that calculate the exact, full LJ param C6/r^6 within the cut-off,
- * as the current nbnxn kernels do.
- */
+ /* When we deprecate the group kernels the code below can go too */
if (fr->vdwtype == evdwPME && fr->cutoff_scheme == ecutsGROUP)
{
/* Calculate self-interaction coefficient (assuming that
#include "typedefs.h"
#include "names.h"
#include "gromacs/utility/smalloc.h"
-#include "gmx_fatal.h"
-#include "gromacs/fileio/futil.h"
-#include "xvgr.h"
-#include "vec.h"
-#include "main.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/math/vec.h"
#include "network.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "force.h"
#include "gromacs/fileio/gmxfio.h"
#include "macros.h"
real *table_v,
real *table_fdv0,
int ntab,
- real dx,
+ double dx,
real beta,
real_space_grid_contribution_computer v_lr)
{
sfree(td->f);
}
- static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr)
+ static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
+ gmx_bool b14only)
{
/* Fill the table according to the formulas in the manual.
* In principle, we only need the potential and the second
int i;
double reppow, p;
double r1, rc, r12, r13;
- double r, r2, r6, rc6;
+ double r, r2, r6, rc2, rc6, rc12;
double expr, Vtab, Ftab;
/* Parameters for David's function */
double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
/* Parameters for the switching function */
double ksw, swi, swi1;
/* Temporary parameters */
- gmx_bool bSwitch, bShift;
+ gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
double ewc = fr->ewaldcoeff_q;
double ewclj = fr->ewaldcoeff_lj;
+ double Vcut = 0;
- bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
- (tp == etabCOULSwitch) ||
- (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch));
-
- bShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
- (tp == etabShift));
+ if (b14only)
+ {
+ bPotentialSwitch = FALSE;
+ bForceSwitch = FALSE;
+ bPotentialShift = FALSE;
+ }
+ else
+ {
+ bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
+ (tp == etabCOULSwitch) ||
+ (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
+ bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
+ (tp == etabShift) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodFORCESWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodFORCESWITCH)));
+ bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
+ }
reppow = fr->reppow;
r1 = fr->rvdw_switch;
rc = fr->rvdw;
}
- if (bSwitch)
+ if (bPotentialSwitch)
{
ksw = 1.0/(pow5(rc-r1));
}
{
ksw = 0.0;
}
- if (bShift)
+ if (bForceSwitch)
{
if (tp == etabShift)
{
fp = xvgropen("switch.xvg", "switch", "r", "s");
#endif
+ if (bPotentialShift)
+ {
+ rc2 = rc*rc;
+ rc6 = 1.0/(rc2*rc2*rc2);
+ if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ {
+ rc12 = rc6*rc6;
+ }
+ else
+ {
+ rc12 = pow(rc, -reppow);
+ }
+
+ switch (tp)
+ {
+ case etabLJ6:
+ /* Dispersion */
+ Vcut = -rc6;
+ break;
+ case etabLJ6Ewald:
+ Vcut = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + pow4(ewclj)*rc2*rc2/2);
+ break;
+ case etabLJ12:
+ /* Repulsion */
+ Vcut = rc12;
+ break;
+ case etabCOUL:
+ Vcut = 1.0/rc;
+ break;
+ case etabEwald:
+ case etabEwaldSwitch:
+ Vtab = gmx_erfc(ewc*rc)/rc;
+ break;
+ case etabEwaldUser:
+ /* Only calculate minus the reciprocal space contribution */
+ Vtab = -gmx_erf(ewc*rc)/rc;
+ break;
+ case etabRF:
+ case etabRF_ZERO:
+ /* No need for preventing the usage of modifiers with RF */
+ Vcut = 0.0;
+ break;
+ case etabEXPMIN:
+ Vcut = exp(-rc);
+ break;
+ default:
+ gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
+ tprops[tp].name, __FILE__, __LINE__);
+ }
+ }
+
for (i = td->nx0; (i < td->nx); i++)
{
r = td->x[i];
}
Vtab = 0.0;
Ftab = 0.0;
- if (bSwitch)
+ if (bPotentialSwitch)
{
/* swi is function, swi1 1st derivative and swi2 2nd derivative */
/* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
tp, __FILE__, __LINE__);
}
- if (bShift)
+ if (bForceSwitch)
{
/* Normal coulomb with cut-off correction for potential */
if (r < rc)
Ftab += A*r12 + B*r13;
}
}
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ }
+ if (bPotentialShift)
+ {
+ if (r < rc)
+ {
+ Vtab -= Vcut;
+ }
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
}
if (ETAB_USER(tp))
Ftab += td->f[i];
}
- if ((r > r1) && bSwitch)
+ if (bPotentialSwitch)
{
- Ftab = Ftab*swi - Vtab*swi1;
- Vtab = Vtab*swi;
+ if (r >= rc)
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ else if (r > r1)
+ {
+ Ftab = Ftab*swi - Vtab*swi1;
+ Vtab = Vtab*swi;
+ }
}
-
/* Convert to single precision when we store to mem */
td->v[i] = Vtab;
td->f[i] = Ftab;
gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
}
- switch (fr->vdw_modifier)
+ /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
+ * to the original interaction forms when we fill the table, so we only check cutoffs here.
+ */
+ if (fr->vdwtype == evdwCUT)
{
- case eintmodNONE:
- case eintmodPOTSHIFT:
- case eintmodEXACTCUTOFF:
- /* No modification */
- break;
- case eintmodPOTSWITCH:
- tabsel[etiLJ6] = etabLJ6Switch;
- tabsel[etiLJ12] = etabLJ12Switch;
- break;
- case eintmodFORCESWITCH:
- tabsel[etiLJ6] = etabLJ6Shift;
- tabsel[etiLJ12] = etabLJ12Shift;
- break;
- default:
- gmx_incons("Unsupported vdw_modifier");
+ switch (fr->vdw_modifier)
+ {
+ case eintmodNONE:
+ case eintmodPOTSHIFT:
+ case eintmodEXACTCUTOFF:
+ /* No modification */
+ break;
+ case eintmodPOTSWITCH:
+ tabsel[etiLJ6] = etabLJ6Switch;
+ tabsel[etiLJ12] = etabLJ12Switch;
+ break;
+ case eintmodFORCESWITCH:
+ tabsel[etiLJ6] = etabLJ6Shift;
+ tabsel[etiLJ12] = etabLJ12Shift;
+ break;
+ default:
+ gmx_incons("Unsupported vdw_modifier");
+ }
}
}
}
init_table(nx, nx0,
(tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
&(td[k]), !bReadTab);
- fill_table(&(td[k]), tabsel[k], fr);
+ fill_table(&(td[k]), tabsel[k], fr, b14only);
if (out)
{
fprintf(out, "%s table with %d data points for %s%s.\n"
#include <math.h>
#include "types/commrec.h"
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
#include "typedefs.h"
#include "nrnb.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "macros.h"
-#include "vec.h"
-#include "main.h"
+#include "gromacs/math/vec.h"
#include "update.h"
#include "gromacs/random/random.h"
-#include "mshift.h"
#include "tgroup.h"
#include "force.h"
#include "names.h"
#include "gmx_omp_nthreads.h"
#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
-#include "gromacs/pulling/pull.h"
+#include "gromacs/utility/smalloc.h"
/*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
/*#define STARTFROMDT2*/
}
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;
#include <math.h>
#include <string.h>
-#include "gromacs/math/utilities.h"
-#include "sysstuff.h"
+
#include "typedefs.h"
#include "macros.h"
-#include "gromacs/utility/smalloc.h"
#include "force.h"
-#include "gromacs/fileio/filenm.h"
#include "nrnb.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
+
+#include "gromacs/fileio/filenm.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/smalloc.h"
void make_wall_tables(FILE *fplog, const output_env_t oenv,
const t_inputrec *ir, const char *tabfn,
}
else
{
- lamfac = 0;
+ lamfac = lambda;
type = md->typeB;
}
}
#include <config.h>
#endif
-
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
-#include "gromacs/fileio/futil.h"
-#include "index.h"
-#include "gromacs/fileio/gmxfio.h"
-#include "vec.h"
+#include <string.h>
+
+#include "gromacs/utility/futil.h"
#include "typedefs.h"
#include "types/commrec.h"
#include "network.h"
-#include "gromacs/fileio/filenm.h"
-#include <string.h>
-#include "gromacs/utility/smalloc.h"
#include "pull.h"
-#include "xvgr.h"
#include "names.h"
-#include "pbc.h"
-#include "mtop_util.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/mtop_util.h"
#include "mdrun.h"
#include "gmx_ga2la.h"
#include "copyrite.h"
#include "macros.h"
-#include "vec.h"
+
+#include "gromacs/fileio/filenm.h"
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/smalloc.h"
static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
{
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]],
#include <math.h>
#include <stdio.h>
#include <string.h>
-#include "main.h"
#include "macros.h"
#include "gromacs/utility/smalloc.h"
-#include "gromacs/fileio/futil.h"
-#include "sysstuff.h"
+#include "gromacs/utility/futil.h"
#include "txtdump.h"
-#include "gmx_fatal.h"
+#include "gromacs/utility/fatalerror.h"
#include "names.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/enxio.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/cstringutil.h"
static void cmp_int(FILE *fp, const char *s, int index, int i1, int i2)
cmp_int(fp, "inputrec->cutoff_scheme", -1, ir1->cutoff_scheme, ir2->cutoff_scheme);
cmp_int(fp, "inputrec->ns_type", -1, ir1->ns_type, ir2->ns_type);
cmp_int(fp, "inputrec->nstlist", -1, ir1->nstlist, ir2->nstlist);
- cmp_int(fp, "inputrec->ndelta", -1, ir1->ndelta, ir2->ndelta);
cmp_int(fp, "inputrec->nstcomm", -1, ir1->nstcomm, ir2->nstcomm);
cmp_int(fp, "inputrec->comm_mode", -1, ir1->comm_mode, ir2->comm_mode);
- cmp_int(fp, "inputrec->nstcheckpoint", -1, ir1->nstcheckpoint, ir2->nstcheckpoint);
cmp_int(fp, "inputrec->nstlog", -1, ir1->nstlog, ir2->nstlog);
cmp_int(fp, "inputrec->nstxout", -1, ir1->nstxout, ir2->nstxout);
cmp_int(fp, "inputrec->nstvout", -1, ir1->nstvout, ir2->nstvout);
cmp_real(fp, "inputrec->ewald_rtol", -1, ir1->ewald_rtol, ir2->ewald_rtol, ftol, abstol);
cmp_int(fp, "inputrec->ewald_geometry", -1, ir1->ewald_geometry, ir2->ewald_geometry);
cmp_real(fp, "inputrec->epsilon_surface", -1, ir1->epsilon_surface, ir2->epsilon_surface, ftol, abstol);
- cmp_int(fp, "inputrec->bOptFFT", -1, ir1->bOptFFT, ir2->bOptFFT);
cmp_int(fp, "inputrec->bContinuation", -1, ir1->bContinuation, ir2->bContinuation);
cmp_int(fp, "inputrec->bShakeSOR", -1, ir1->bShakeSOR, ir2->bShakeSOR);
cmp_int(fp, "inputrec->etc", -1, ir1->etc, ir2->etc);
cmp_real(fp, "inputrec->orires_fc", -1, ir1->orires_fc, ir2->orires_fc, ftol, abstol);
cmp_real(fp, "inputrec->orires_tau", -1, ir1->orires_tau, ir2->orires_tau, ftol, abstol);
cmp_int(fp, "inputrec->nstorireout", -1, ir1->nstorireout, ir2->nstorireout);
- cmp_real(fp, "inputrec->dihre_fc", -1, ir1->dihre_fc, ir2->dihre_fc, ftol, abstol);
cmp_real(fp, "inputrec->em_stepsize", -1, ir1->em_stepsize, ir2->em_stepsize, ftol, abstol);
cmp_real(fp, "inputrec->em_tol", -1, ir1->em_tol, ir2->em_tol, ftol, abstol);
cmp_int(fp, "inputrec->niter", -1, ir1->niter, ir2->niter);
}
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))
{
#include <config.h>
#endif
+#include <stdlib.h>
+
#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "sysstuff.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "vcm.h"
#include "mdebin.h"
#include "nrnb.h"
#include "calcmu.h"
-#include "index.h"
#include "vsite.h"
#include "update.h"
#include "ns.h"
#include "md_support.h"
#include "md_logging.h"
#include "network.h"
-#include "xvgr.h"
-#include "physics.h"
#include "names.h"
#include "force.h"
#include "disre.h"
#include "qmmm.h"
#include "domdec.h"
#include "domdec_network.h"
-#include "gromacs/gmxlib/topsort.h"
#include "coulomb.h"
#include "constr.h"
#include "shellfc.h"
#include "gromacs/gmxpreprocess/compute_io.h"
#include "checkpoint.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "sighandler.h"
#include "txtdump.h"
#include "gromacs/utility/cstringutil.h"
#include "types/iteratedconstraints.h"
#include "nbnxn_cuda_data_mgmt.h"
-#include "gromacs/utility/gmxmpi.h"
#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/mdoutf.h"
#include "gromacs/fileio/trajectory_writing.h"
#include "gromacs/fileio/trnio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/xtcio.h"
-#include "gromacs/timing/wallcycle.h"
-#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/swap/swapcoords.h"
-#include "gromacs/imd/imd.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/utility/smalloc.h"
#ifdef GMX_FAHCORE
#include "corewrap.h"
rvec mu_tot;
t_vcm *vcm;
t_state *bufstate = NULL;
- matrix *scale_tot, pcoupl_mu, M, ebox;
+ matrix pcoupl_mu, M;
gmx_nlheur_t nlh;
t_trxframe rerun_fr;
gmx_repl_ex_t repl_ex = NULL;
gmx_ekindata_t *ekind, *ekind_save;
gmx_shellfc_t shellfc;
int count, nconverged = 0;
- real timestep = 0;
- double tcount = 0;
- gmx_bool bConverged = TRUE, bOK, bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
- gmx_bool bAppend;
+ double tcount = 0;
+ gmx_bool bConverged = TRUE, bOK, bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
gmx_bool bResetCountersHalfMaxH = FALSE;
gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
gmx_bool bUpdateDoLR;
double cycles;
real saved_conserved_quantity = 0;
real last_ekin = 0;
- int iter_i;
t_extmass MassQ;
int **trotter_seq;
char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
/* Check for special mdrun options */
bRerunMD = (Flags & MD_RERUN);
- bAppend = (Flags & MD_APPENDFILES);
if (Flags & MD_RESETCOUNTERSHALFWAY)
{
if (ir->nsteps > 0)
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 */
bStateFromTPX = !bStateFromCP;
bInitStep = bFirstStep && (bStateFromTPX || bVV);
bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
- bLastStep = FALSE;
bSumEkinhOld = FALSE;
- bDoReplEx = FALSE;
bExchanged = FALSE;
bNeedRepartition = FALSE;
step = ir->init_step;
step_rel = 0;
- if (ir->nstlist == -1)
- {
- init_nlistheuristics(&nlh, bGStatEveryStep, step);
- }
+ init_nlistheuristics(&nlh, bGStatEveryStep, step);
if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
{
#include "gromacs/legacyheaders/checkpoint.h"
#include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/legacyheaders/main.h"
#include "gromacs/legacyheaders/mdrun.h"
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/filenm.h"
+#include "gromacs/utility/fatalerror.h"
int gmx_mdrun(int argc, char *argv[])
{
"[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",
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
+
+#include "repl_ex.h"
+
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>
-#include "repl_ex.h"
+
#include "network.h"
#include "gromacs/random/random.h"
#include "gromacs/utility/smalloc.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "copyrite.h"
-#include "macros.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "names.h"
#include "domdec.h"
+#include "main.h"
#include "gromacs/random/random.h"
#define PROBABILITYCUTOFF 100
{
real *qall;
gmx_bool bDiff;
- int i, s;
+ int s;
snew(qall, ms->nsim);
qall[re->repl] = q;
const t_inputrec *ir,
int nst, int nex, int init_seed)
{
- real temp, pres;
+ real pres;
int i, j, k;
struct gmx_repl_ex *re;
gmx_bool bTemp;
{
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;
}
real time)
{
int m, i, j, a, b, ap, bp, i0, i1, tmp;
- real ediff = 0, delta = 0, dpV = 0;
+ real delta = 0;
gmx_bool bPrint, bMultiEx;
gmx_bool *bEx = re->bEx;
real *prob = re->prob;
gmx_bool bEpot = FALSE;
gmx_bool bDLambda = FALSE;
gmx_bool bVol = FALSE;
- 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);
++ 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;
t_state *state, gmx_enerdata_t *enerd,
t_state *state_local, gmx_int64_t step, real time)
{
- int i, j;
+ int j;
int replica_id = 0;
int exchange_partner;
int maxswap = 0;
{
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
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
+
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+
+#include <algorithm>
+
+#include <assert.h>
#include <signal.h>
#include <stdlib.h>
+#include <string.h>
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
-#include <string.h>
-#include <assert.h>
#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "sysstuff.h"
#include "copyrite.h"
#include "force.h"
#include "mdrun.h"
#include "constr.h"
#include "mvdata.h"
#include "checkpoint.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "sighandler.h"
#include "txtdump.h"
#include "gmx_detect_hardware.h"
#include "gmx_omp_nthreads.h"
#include "gromacs/gmxpreprocess/calc_verletbuf.h"
-#include "gmx_fatal_collective.h"
#include "membed.h"
-#include "macros.h"
#include "gmx_thread_affinity.h"
#include "inputrec.h"
+#include "main.h"
+#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/fileio/tpxio.h"
+#include "gromacs/math/vec.h"
#include "gromacs/mdlib/nbnxn_search.h"
#include "gromacs/mdlib/nbnxn_consts.h"
-#include "gromacs/timing/wallcycle.h"
-#include "gromacs/utility/gmxmpi.h"
-#include "gromacs/utility/gmxomp.h"
-#include "gromacs/swap/swapcoords.h"
-#include "gromacs/essentialdynamics/edsam.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/pulling/pull_rotation.h"
+#include "gromacs/swap/swapcoords.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/utility/smalloc.h"
#ifdef GMX_FAHCORE
#include "corewrap.h"
else if (hw_opt->nthreads_omp > 0)
{
/* Here we could oversubscribe, when we do, we issue a warning later */
- nthreads_tmpi = max(1, nthreads_tot/hw_opt->nthreads_omp);
+ nthreads_tmpi = std::max(1, nthreads_tot/hw_opt->nthreads_omp);
}
else
{
const int nthreads_omp_always_faster = 4;
const int nthreads_omp_always_faster_Nehalem = 12;
const int nthreads_omp_always_faster_SandyBridge = 16;
- const int first_model_Nehalem = 0x1A;
- const int first_model_SandyBridge = 0x2A;
gmx_bool bIntel_Family6;
bIntel_Family6 =
{
int nthreads_hw, nthreads_tot_max, nthreads_tmpi, nthreads_new, ngpu;
int min_atoms_per_mpi_thread;
- char *env;
- char sbuf[STRLEN];
gmx_bool bCanUseGPU;
if (hw_opt->nthreads_tmpi > 0)
{
/* the thread number was chosen automatically, but there are too many
threads (too few atoms per thread) */
- nthreads_new = max(1, mtop->natoms/min_atoms_per_mpi_thread);
+ nthreads_new = std::max(1, mtop->natoms/min_atoms_per_mpi_thread);
/* Avoid partial use of Hyper-Threading */
if (gmx_cpuid_x86_smt(hwinfo->cpuid_info) == GMX_CPUID_X86_SMT_ENABLED &&
/* We determine the extra cost of the non-bonded kernels compared to
* a reference nstlist value of 10 (which is the default in grompp).
*/
-static const int nbnxn_reference_nstlist = 10;
+static const int nbnxnReferenceNstlist = 10;
/* The values to try when switching */
const int nstlist_try[] = { 20, 25, 40 };
#define NNSTL sizeof(nstlist_try)/sizeof(nstlist_try[0])
float listfac_ok, listfac_max;
int nstlist_orig, nstlist_prev;
verletbuf_list_setup_t ls;
- real rlist_nstlist10, rlist_inc, rlist_ok, rlist_max;
+ real rlistWithReferenceNstlist, rlist_inc, rlist_ok, rlist_max;
real rlist_new, rlist_prev;
- int nstlist_ind = 0;
+ size_t nstlist_ind = 0;
t_state state_tmp;
gmx_bool bBox, bDD, bCont;
const char *nstl_gpu = "\nFor optimal performance with a GPU nstlist (now %d) should be larger.\nThe optimum depends on your CPU and GPU resources.\nYou might want to try several nstlist values.\n";
const char *box_err = "Can not increase nstlist because the box is too small";
const char *dd_err = "Can not increase nstlist because of domain decomposition limitations";
char buf[STRLEN];
+ const float oneThird = 1.0f / 3.0f;
if (nstlist_cmdline <= 0)
{
verletbuf_get_list_setup(bGPU, &ls);
/* Allow rlist to make the list a given factor larger than the list
- * would be with nstlist=10.
+ * would be with the reference value for nstlist (10).
*/
nstlist_prev = ir->nstlist;
- ir->nstlist = 10;
+ ir->nstlist = nbnxnReferenceNstlist;
calc_verlet_buffer_size(mtop, det(box), ir, -1, &ls, NULL,
- &rlist_nstlist10);
+ &rlistWithReferenceNstlist);
ir->nstlist = nstlist_prev;
/* Determine the pair list size increase due to zero interactions */
rlist_inc = nbnxn_get_rlist_effective_inc(ls.cluster_size_j,
mtop->natoms/det(box));
- rlist_ok = (rlist_nstlist10 + rlist_inc)*pow(listfac_ok, 1.0/3.0) - rlist_inc;
- rlist_max = (rlist_nstlist10 + rlist_inc)*pow(listfac_max, 1.0/3.0) - rlist_inc;
+ rlist_ok = (rlistWithReferenceNstlist + rlist_inc)*pow(listfac_ok, oneThird) - rlist_inc;
+ rlist_max = (rlistWithReferenceNstlist + rlist_inc)*pow(listfac_max, oneThird) - rlist_inc;
if (debug)
{
fprintf(debug, "nstlist tuning: rlist_inc %.3f rlist_ok %.3f rlist_max %.3f\n",
t_inputrec *ir,
gmx_mtop_t *mtop, real box_vol)
{
- char *conv_mesg = "Converting input file with group cut-off scheme to the Verlet cut-off scheme";
+ const char *conv_mesg = "Converting input file with group cut-off scheme to the Verlet cut-off scheme";
md_print_warn(NULL, fplog, "%s\n", conv_mesg);
rlist_fac = 1 + verlet_buffer_ratio_nodynamics;
}
ir->verletbuf_tol = -1;
- ir->rlist = rlist_fac*max(ir->rvdw, ir->rcoulomb);
+ ir->rlist = rlist_fac*std::max(ir->rvdw, ir->rcoulomb);
}
gmx_mtop_remove_chargegroups(mtop);
gmx_bool bIsPPrankUsingGPU;
char gpu_err_str[STRLEN];
- bIsPPrankUsingGPU = (cr->duty & DUTY_PP) && fr->nbv != NULL && fr->nbv->bUseGPU;
+ bIsPPrankUsingGPU = (cr->duty & DUTY_PP) && fr != NULL && fr->nbv != NULL && fr->nbv->bUseGPU;
if (bIsPPrankUsingGPU)
{
const char *deviceOptions, int imdport, unsigned long Flags)
{
gmx_bool bForceUseGPU, bTryUseGPU;
- double nodetime = 0, realtime;
t_inputrec *inputrec;
t_state *state = NULL;
matrix box;
gmx_ddbox_t ddbox = {0};
int npme_major, npme_minor;
- real tmpr1, tmpr2;
t_nrnb *nrnb;
gmx_mtop_t *mtop = NULL;
t_mdatoms *mdatoms = NULL;
gmx_pme_t *pmedata = NULL;
gmx_vsite_t *vsite = NULL;
gmx_constr_t constr;
- int i, m, nChargePerturbed = -1, nTypePerturbed = 0, status, nalloc;
- char *gro;
+ int nChargePerturbed = -1, nTypePerturbed = 0, status;
gmx_wallcycle_t wcycle;
gmx_bool bReadEkin;
- int list;
gmx_walltime_accounting_t walltime_accounting = NULL;
int rc;
gmx_int64_t reset_counters;
gmx_edsam_t ed = NULL;
- t_commrec *cr_old = cr;
int nthreads_pme = 1;
int nthreads_pp = 1;
gmx_membed_t membed = NULL;
}
}
- /* Check for externally set OpenMP affinity and turn off internal
- * pinning if any is found. We need to do this check early to tell
- * thread-MPI whether it should do pinning when spawning threads.
- * TODO: the above no longer holds, we should move these checks down
- */
- gmx_omp_check_thread_affinity(fplog, cr, hw_opt);
-
/* Check and update the hardware options for internal consistency */
check_and_update_hw_opt_1(hw_opt, SIMMASTER(cr));
+ /* Early check for externally set process affinity. */
+ gmx_check_thread_affinity_set(fplog, cr,
+ hw_opt, hwinfo->nthreads_hw_avail, FALSE);
if (SIMMASTER(cr))
{
-#ifdef GMX_THREAD_MPI
- /* Early check for externally set process affinity.
- * With thread-MPI this is needed as pinning might get turned off,
- * which needs to be known before starting thread-MPI.
- * With thread-MPI hw_opt is processed here on the master rank
- * and passed to the other ranks later, so we only do this on master.
- */
- gmx_check_thread_affinity_set(fplog,
- NULL,
- hw_opt, hwinfo->nthreads_hw_avail, FALSE);
-#endif
#ifdef GMX_THREAD_MPI
if (cr->npmenodes > 0 && hw_opt->nthreads_tmpi <= 0)
if (hw_opt->nthreads_tmpi > 1)
{
+ t_commrec *cr_old = cr;
/* now start the threads. */
cr = mdrunner_start_threads(hw_opt, fplog, cr_old, nfile, fnm,
oenv, bVerbose, bCompact, nstglobalcomm,
/* Initialize per-physical-node MPI process/thread ID and counters. */
gmx_init_intranode_counters(cr);
-
#ifdef GMX_MPI
+ if (MULTISIM(cr))
+ {
+ md_print_info(cr, fplog,
+ "This is simulation %d out of %d running as a composite Gromacs\n"
+ "multi-simulation job. Setup for this simulation:\n\n",
+ cr->ms->sim, cr->ms->nsim);
+ }
md_print_info(cr, fplog, "Using %d MPI %s\n",
cr->nnodes,
#ifdef GMX_THREAD_MPI
if (DOMAINDECOMP(cr))
{
+ GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP");
dd_init_bondeds(fplog, cr->dd, mtop, vsite, inputrec,
Flags & MD_DDBONDCHECK, fr->cginfo_mb);
}
else
{
+ GMX_RELEASE_ASSERT(pmedata, "pmedata was NULL while cr->duty was not DUTY_PP");
/* do PME only */
walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntPME));
gmx_pmeonly(*pmedata, cr, nrnb, wcycle, walltime_accounting, ewaldcoeff_q, ewaldcoeff_lj, inputrec);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff