Clashes in set_pull_init signature.
Issues with new pull checks from LocalAtomSet refactoring.
Transfered change to PME-on-GPU availability to new
location.
Ignored changes to pre-submit, because no longer appropriate.
Adopted changes to gmx solvate into master. The new test
coverage finds a memory leak that is fixed in this patch with
a refactoring of sort_molecule and add_solv.
Needed to include domdec/collect.h now used in master because
a function call from it was introduced in release-2018.
Change-Id: I5f7fcdf52a6093b455bd6e31c264696d88ced2ac
# Test oldest supported CUDA
# Test oldest supported Ubuntu
# Test MPI with CUDA
+# Test cmake version from before new FindCUDA support (in 3.8)
# Test MPMD PME with library MPI
# Test recent cmake (3.7+), to cover minor FindCUDA changes from 3.7.0
gcc-4.8 gpu cuda-7.0 cmake-3.8.1 mpi npme=1 nranks=2 openmp
if(NOT GMX_OPENMP)
message(WARNING "To use GPU acceleration efficiently, mdrun requires OpenMP multi-threading. Without OpenMP a single CPU core can be used with a GPU which is not optimal. Note that with MPI multiple processes can be forced to use a single GPU, but this is typically inefficient. You need to set both C and C++ compilers that support OpenMP (CC and CXX environment variables, respectively) when using GPUs.")
endif()
-
- if(NOT GMX_CLANG_CUDA)
- gmx_check_if_changed(GMX_CHECK_NVCC CUDA_NVCC_EXECUTABLE CUDA_HOST_COMPILER CUDA_NVCC_FLAGS)
-
- if(GMX_CHECK_NVCC OR NOT GMX_NVCC_WORKS)
- message(STATUS "Check for working NVCC/C compiler combination")
- execute_process(COMMAND ${CUDA_NVCC_EXECUTABLE} -ccbin ${CUDA_HOST_COMPILER} -c ${CUDA_NVCC_FLAGS} ${CMAKE_SOURCE_DIR}/cmake/TestCUDA.cu
- RESULT_VARIABLE _cuda_test_res
- OUTPUT_VARIABLE _cuda_test_out
- ERROR_VARIABLE _cuda_test_err
- OUTPUT_STRIP_TRAILING_WHITESPACE)
-
- if(${_cuda_test_res})
- message(STATUS "Check for working NVCC/C compiler combination - broken")
- if(${_cuda_test_err} MATCHES "nsupported")
- message(FATAL_ERROR "NVCC/C compiler combination does not seem to be supported. CUDA frequently does not support the latest versions of the host compiler, so you might want to try an earlier C/C++ compiler version and make sure your CUDA compiler and driver are as recent as possible.")
- else()
- message(FATAL_ERROR "CUDA compiler does not seem to be functional.")
- endif()
- elseif(NOT GMX_CUDA_TEST_COMPILER_QUIETLY)
- message(STATUS "Check for working NVCC/C compiler combination - works")
- set(GMX_NVCC_WORKS TRUE CACHE INTERNAL "Nvcc can compile a trivial test program")
- endif()
- endif() # GMX_CHECK_NVCC
- endif() #GMX_CLANG_CUDA
endif() # GMX_GPU
option(GMX_CUDA_NB_SINGLE_COMPILATION_UNIT "Whether to compile the CUDA non-bonded module using a single compilation unit." OFF)
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
+# Copyright (c) 2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
# Manage CUDA nvcc compilation configuration, try to be smart to ease the users'
# pain as much as possible:
# - use the CUDA_HOST_COMPILER if defined by the user, otherwise
-# - auto-detect compatible nvcc host compiler and set nvcc -ccbin (if not MPI wrapper)
-# - set icc compatibility mode to gcc 4.8.1
+# - check if nvcc works with CUDA_HOST_COMPILER and the generated nvcc and C++ flags
+#
# - (advanced) variables set:
-# * CUDA_HOST_COMPILER - the host compiler for nvcc (only with cmake <2.8.10)
# * CUDA_HOST_COMPILER_OPTIONS - the full host-compiler related option list passed to nvcc
#
# Note that from CMake 2.8.10 FindCUDA defines CUDA_HOST_COMPILER internally,
# so we won't set it ourselves, but hope that the module does a good job.
-gmx_check_if_changed(CUDA_HOST_COMPILER_CHANGED CUDA_HOST_COMPILER)
-
-# CUDA_HOST_COMPILER changed hence it is not auto-set anymore
-if (CUDA_HOST_COMPILER_CHANGED AND CUDA_HOST_COMPILER_AUTOSET)
- unset(CUDA_HOST_COMPILER_AUTOSET CACHE)
-endif()
-
# glibc 2.23 changed string.h in a way that breaks CUDA compilation in
# many projects, but which has a trivial workaround. It would be nicer
# to compile with nvcc and see that the workaround is necessary and
endif()
endfunction()
+gmx_check_if_changed(CUDA_HOST_COMPILER_CHANGED CUDA_HOST_COMPILER)
+
# set up host compiler and its options
if(CUDA_HOST_COMPILER_CHANGED)
set(CUDA_HOST_COMPILER_OPTIONS "")
# assemble the CUDA host compiler flags
list(APPEND GMX_CUDA_NVCC_FLAGS "${CUDA_HOST_COMPILER_OPTIONS}")
+string(TOUPPER "${CMAKE_BUILD_TYPE}" _build_type)
+gmx_check_if_changed(_cuda_nvcc_executable_or_flags_changed CUDA_NVCC_EXECUTABLE CUDA_NVCC_FLAGS CUDA_NVCC_FLAGS_${_build_type})
+
+if(_cuda_nvcc_executable_or_flags_changed OR CUDA_HOST_COMPILER_CHANGED OR NOT GMX_NVCC_WORKS)
+ message(STATUS "Check for working NVCC/C compiler combination")
+ execute_process(COMMAND ${CUDA_NVCC_EXECUTABLE} -ccbin ${CUDA_HOST_COMPILER} -c ${CUDA_NVCC_FLAGS} ${CUDA_NVCC_FLAGS_${_build_type}} ${CMAKE_SOURCE_DIR}/cmake/TestCUDA.cu
+ RESULT_VARIABLE _cuda_test_res
+ OUTPUT_VARIABLE _cuda_test_out
+ ERROR_VARIABLE _cuda_test_err
+ OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+ if(${_cuda_test_res})
+ message(${_cuda_test_err})
+ message(STATUS "Check for working NVCC/C compiler combination - broken")
+ if(${_cuda_test_err} MATCHES "nsupported")
+ message(FATAL_ERROR "NVCC/C compiler combination does not seem to be supported. CUDA frequently does not support the latest versions of the host compiler, so you might want to try an earlier C/C++ compiler version and make sure your CUDA compiler and driver are as recent as possible.")
+ else()
+ message(FATAL_ERROR "CUDA compiler does not seem to be functional.")
+ endif()
+ elseif(NOT GMX_CUDA_TEST_COMPILER_QUIETLY)
+ message(STATUS "Check for working NVCC/C compiler combination - works")
+ set(GMX_NVCC_WORKS TRUE CACHE INTERNAL "Nvcc can compile a trivial test program")
+ endif()
+endif() # GMX_CHECK_NVCC
+
+
# The flags are set as local variables which shadow the cache variables. The cache variables
# (can be set by the user) are appended. This is done in a macro to set the flags when all
# host compiler flags are already set.
# build the regressiontests tarball with all the right naming. The
# naming affects the md5sum that has to go here, and if it isn't right
# release workflow will report a failure.
-set(REGRESSIONTEST_MD5SUM "135149af467a37714a92bc29801cafda" CACHE INTERNAL "MD5 sum of the regressiontests tarball for this GROMACS version")
+set(REGRESSIONTEST_MD5SUM "368aabf0eb4e6fbc9abb7961475f66de" CACHE INTERNAL "MD5 sum of the regressiontests tarball for this GROMACS version")
math(EXPR GMX_VERSION_NUMERIC
"${GMX_VERSION_MAJOR}*10000 + ${GMX_VERSION_PATCH}")
release-notes/removed-features.rst
release-notes/portability.rst
release-notes/miscellaneous.rst
+ release-notes/2018/2018.4.rst
release-notes/2018/2018.3.rst
release-notes/2018/2018.2.rst
release-notes/2018/2018.1.rst
GROMACS 2018.3 release notes
----------------------------
-This version was released on TODO, 2018. These release notes document
+This version was released on August 23, 2018. These release notes document
the changes that have taken place in GROMACS since version 2018.2, to fix known
-issues. It also incorporates all fixes made in version TODO and
+issues. It also incorporates all fixes made in version 2016.5 and
earlier, which you can find described in the :ref:`release-notes`.
Fixes where mdrun could behave incorrectly
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Multi-domain GPU runs can no longer miss pair interactions
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+With systems with empty space in the unit cell, GPU runs with domain
+decomposition would not compute LJ and Coulomb interactions between
+domains when there we no interactions between domains on a rank at some
+point in time.
+
+ - This bug only affects simulations running on GPUs with domain decomposition
+ and containing empty regions of space that can lead to domains being empty.
+ - Possible observations of this bug may have been seemingly random failures
+ of calculations that where not reproducible when restarting a simulation
+ from a checkpoint file, as the domain would then again be filled properly
+ if interactions are present at the beginning.
+ - It is unlikely that this bug will have unnoticed effects on all but
+ very short simulations, as the missing interactions will inevitable lead
+ to simulation instability and crashes.
+ - If a simulation that crashed due to this bug is restarted it can contain
+ a small region around the crash that will be unphysical due to some
+ interactions not being calculated just before the crash itself.
+
+**This is a critical fix and users of 2018.x series that run on GPUs should
+update to this point release**
+
+:issue:`2502`
+
+Fix Conjugate Gradient assertion failure at end of minimization
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When the final step coincided with a coordinate output step,
+conjugate gradient minimization would exit with an assertion failure
+instead of writing confout.gro.
+
+:issue:`2554`
+
+Multi-domain Conjugate Gradient minimimization no longer segfaults.
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2554`
+
+Fix pairlist buffer with Brownian Dynamics
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+With Brownian Dynamics and bd-fric > 0, the Verlet pairlist buffer would
+be determined with incorrect masses for constrained atoms and virtual
+sites. This would lead to a too small buffer for typical atomistic
+systems with constraints.
+
+:issue:`2613`
+
+Avoid "atom moved to far" errors
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The introduction of the dual pair list has led to larger nstlist values,
+which leads to larger atom displacements between domain decomposition
+steps. This has made it more likely that the errors
+"An atom moved too far between two domain decomposition steps" and
+"N particles communicated to PME rank M are more than 2/3 times the cut-off
+out of the domain decomposition cell ..." appear for stable systems.
+Now atom displacements are correctly taken into account when determining
+the minimum cell size, so these errors should only appear for unstable systems.
+
+:issue:`2614`
+
+grompp now checks that pull groups are not close to half the box size
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Pull groups that use a reference atom for periodic boundary treatment
+should have all their atoms well within half the box size of this reference.
+When this is not the case, grompp will issue a warning.
+
+:issue:`2397`
+
+Fixed segmentation fault in mdrun with QM/MM ONIOM scheme
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2617`
+
+Correctly specified that PME on GPUs is only supported for dynamical integrators
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Previously PME on GPU support could run (but fail) for energy
+minimization and normal-mode analysis runs.
+
+:issue:`2578`
+
Fixes for ``gmx`` tools
^^^^^^^^^^^^^^^^^^^^^^^
+Fixed syntax error in make_gromos_rtp.py
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2557`
+
+Fix gmx solvate topology updating
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Removed hard coded solvent names to allow updates to topology based on
+solvent molecule information. Also allows updating with multiple solvent
+types.
+
+:issue:`1929`
+
+Fix bfactor output error caused by fix for :issue:`2511`
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+The fix for the PQR file output broke the output of bfactors from other tools.
+
+:issue:`2575`
+
+Made sure that gmx rms can skip values
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+When requested to skip values, gmx rms would still output all values despite
+the option. Now it only outputs values that are requested to be processed.
+
+:issue:`2565`
+
+Fix trjconv when not providing structure file
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+trjconv would fail with a segmentation violation when running without any structure
+file due to incomplete initialization of the topology data structure. This fix adds
+the missing checks that prevents the failure.
+
+:issue:`2619`
+
+Fix enforced rotation energy output
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+
Fixes to improve portability
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Fix nvcc host compiler check triggering
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2583`
+
+
+Report up to date hwloc version information
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2591`
+
+
+Disable single compilation unit with CUDA 9.0
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2561`
+
+
Miscellaneous
^^^^^^^^^^^^^
+
+Avoid aborting mdrun when GPU sanity check detects errors
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2415`
+
+
+Improve OpenCL kernel performance on AMD Vega GPUs
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+The OpenCL kernel optimization flags did not explicitly turn off denorm handling
+which could lead to performance loss. The optimization is now explicitly turned
+on both for consistency with CUDA and performance reasons.
+On AMD Vega GPUs (with ROCm) kernel performance improves by up to 30%.
+
+
--- /dev/null
+GROMACS 2018.4 release notes
+----------------------------
+
+This version was released on TODO, 2018. These release notes document
+the changes that have taken place in GROMACS since version 2018.3, to fix known
+issues. It also incorporates all fixes made in version TODO and
+earlier, which you can find described in the :ref:`release-notes`.
+
+Fixes where mdrun could behave incorrectly
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Fixes for ``gmx`` tools
+^^^^^^^^^^^^^^^^^^^^^^^
+
+Fixes to improve portability
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Miscellaneous
+^^^^^^^^^^^^^
.. toctree::
:maxdepth: 1
+ 2018/2018.4
2018/2018.3
2018/2018.2
2018/2018.1
For more information about GPU tasks, please refer to
:ref:`Types of GPU tasks<gmx-gpu-tasks>`.
+``-pmefft``
+ Allows choosing whether to execute the 3D FFT computation on a CPU or GPU.
+ Can be set to "auto", "cpu", "gpu.".
+ When PME is offloaded to a GPU ``-pmefft gpu`` is the default,
+ and the entire PME calculation is executed on the GPU. However,
+ in some cases, e.g. with a relatively slow or older generation GPU
+ combined with fast CPU cores in a run, moving some work off of the GPU
+ back to the CPU by computing FFTs on the CPU can improve performance.
+
Examples for :ref:`mdrun <gmx mdrun>` on one node
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Examples for :ref:`mdrun <gmx mdrun>` on more than one node
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The examples and explanations for for single-node :ref:`mdrun <gmx mdrun>` are
-still relevant, but ``-nt`` is no longer the way
+still relevant, but ``-ntmpi`` is no longer the way
to choose the number of MPI ranks.
::
- Free energy calculations where charges are perturbed are not supported,
because only single PME grids can be calculated.
+- Only dynamical integrators are supported (ie. leap-frog, Velocity Verlet,
+ stochastic dynamics)
+
- LJ PME is not supported on GPUs.
Assigning tasks to GPUs
for i in range(len(list)):
if list[i] == string:
return i
- print >> sys.stderr "Could not find string",string,"in list of length",len(list)
+ print >> sys.stderr, "Could not find string", string, "in list of length", len(list)
return -1
#--------------------------#
/** Location of GROMACS-specific data files */
#define GMX_INSTALL_GMXDATADIR "@GMX_INSTALL_GMXDATADIR@"
+/** HWLOC version information */
+#define HWLOC_VERSION "@HWLOC_VERSION@"
+
/** CUDA compiler version information */
#define CUDA_COMPILER_INFO "@CUDA_COMPILER_INFO@"
#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/calc_verletbuf.h"
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/constraintrange.h"
#include "gromacs/mdlib/forcerec.h"
}
comm->cutoff_mbody = 0;
+ /* Determine the minimum cell size limit, affected by many factors */
comm->cellsize_limit = 0;
comm->bBondComm = FALSE;
- /* Atoms should be able to move by up to half the list buffer size (if > 0)
- * within nstlist steps. Since boundaries are allowed to displace by half
- * a cell size, DD cells should be at least the size of the list buffer.
+ /* We do not allow home atoms to move beyond the neighboring domain
+ * between domain decomposition steps, which limits the cell size.
+ * Get an estimate of cell size limit due to atom displacement.
+ * In most cases this is a large overestimate, because it assumes
+ * non-interaction atoms.
+ * We set the chance to 1 in a trillion steps.
*/
+ constexpr real c_chanceThatAtomMovesBeyondDomain = 1e-12;
+ const real limitForAtomDisplacement =
+ minCellSizeForAtomDisplacement(*mtop, *ir, c_chanceThatAtomMovesBeyondDomain);
+ if (fplog)
+ {
+ fprintf(fplog,
+ "Minimum cell size due to atom displacement: %.3f nm\n",
+ limitForAtomDisplacement);
+ }
comm->cellsize_limit = std::max(comm->cellsize_limit,
- ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
+ limitForAtomDisplacement);
+
+ /* TODO: PME decomposition currently requires atoms not to be more than
+ * 2/3 of comm->cutoff, which is >=rlist, outside of their domain.
+ * In nearly all cases, limitForAtomDisplacement will be smaller
+ * than 2/3*rlist, so the PME requirement is satisfied.
+ * But it would be better for both correctness and performance
+ * to use limitForAtomDisplacement instead of 2/3*comm->cutoff.
+ * Note that we would need to improve the pairlist buffer case.
+ */
if (comm->bInterCGBondeds)
{
{
errorReasons.emplace_back("group cutoff scheme");
}
- if (EI_TPI(ir->eI))
+ if (!EI_DYNAMICS(ir->eI))
{
- errorReasons.emplace_back("test particle insertion");
+ errorReasons.emplace_back("not a dynamical integrator");
}
return addMessageIfNotSupported(errorReasons, error);
}
case efENT:
case efPQR:
out = gmx_fio_fopen(outfile, "w");
- write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr, TRUE, TRUE);
+ write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr, TRUE, ftp == efPQR);
gmx_fio_fclose(out);
break;
case efESP:
"HIDDENAverage over this distance in the RMSD matrix" }
};
int natoms_trx, natoms_trx2, natoms;
- int i, j, k, m, teller, teller2, tel_mat, tel_mat2;
+ int i, j, k, m;
#define NFRAME 5000
int maxframe = NFRAME, maxframe2 = NFRAME;
real t, *w_rls, *w_rms, *w_rls_m = nullptr, *w_rms_m = nullptr;
}
/* start looping over frames: */
- tel_mat = 0;
- teller = 0;
+ int tel_mat = 0;
+ int teller = 0;
+ int frame = 0;
do
{
if (bPBC)
do_fit(natoms, w_rls, xp, x);
}
- if (teller % freq == 0)
+ if (frame % freq == 0)
{
/* keep frame for matrix calculation */
if (bMat || bBond || bPrev)
}
}
tel_mat++;
- }
- /*calculate energy of root_least_squares*/
- if (bPrev)
- {
- j = tel_mat-prev-1;
- if (j < 0)
- {
- j = 0;
- }
- for (i = 0; i < n_ind_m; i++)
+ /*calculate energy of root_least_squares*/
+ if (bPrev)
{
- copy_rvec(mat_x[j][i], xp[ind_m[i]]);
- }
- if (bReset)
- {
- reset_x(ifit, ind_fit, natoms, nullptr, xp, w_rls);
+ j = tel_mat-prev-1;
+ if (j < 0)
+ {
+ j = 0;
+ }
+ for (i = 0; i < n_ind_m; i++)
+ {
+ copy_rvec(mat_x[j][i], xp[ind_m[i]]);
+ }
+ if (bReset)
+ {
+ reset_x(ifit, ind_fit, natoms, nullptr, xp, w_rls);
+ }
+ if (bFit)
+ {
+ do_fit(natoms, w_rls, x, xp);
+ }
}
- if (bFit)
+ for (j = 0; (j < nrms); j++)
{
- do_fit(natoms, w_rls, x, xp);
+ rls[j][teller] =
+ calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
}
- }
- for (j = 0; (j < nrms); j++)
- {
- rls[j][teller] =
- calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
- }
- if (bNorm)
- {
- for (j = 0; (j < irms[0]); j++)
+ if (bNorm)
{
- rlsnorm[j] +=
- calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp);
+ for (j = 0; (j < irms[0]); j++)
+ {
+ rlsnorm[j] +=
+ calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp);
+ }
}
- }
- if (bMirror)
- {
- if (bFit)
+ if (bMirror)
{
- /*do the least squares fit to mirror of original structure*/
- do_fit(natoms, w_rls, xm, x);
- }
+ if (bFit)
+ {
+ /*do the least squares fit to mirror of original structure*/
+ do_fit(natoms, w_rls, xm, x);
+ }
- for (j = 0; j < nrms; j++)
- {
- rlsm[j][teller] =
- calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
+ for (j = 0; j < nrms; j++)
+ {
+ rlsm[j][teller] =
+ calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
+ }
}
- }
- time[teller] = output_env_conv_time(oenv, t);
+ time[teller] = output_env_conv_time(oenv, t);
- teller++;
+ teller++;
+ }
+ frame++;
if (teller >= maxframe)
{
maxframe += NFRAME;
while (read_next_x(oenv, status, &t, x, box));
close_trx(status);
+ int tel_mat2 = 0;
+ int teller2 = 0;
+ int frame2 = 0;
if (bFile2)
{
snew(time2, maxframe2);
"Second trajectory (%d atoms) does not match the first one"
" (%d atoms)", natoms_trx2, natoms_trx);
}
- tel_mat2 = 0;
- teller2 = 0;
+ frame2 = 0;
do
{
if (bPBC)
do_fit(natoms, w_rls, xp, x);
}
- if (teller2 % freq2 == 0)
+ if (frame2 % freq2 == 0)
{
/* keep frame for matrix calculation */
if (bMat)
}
}
tel_mat2++;
- }
- time2[teller2] = output_env_conv_time(oenv, t);
+ time2[teller2] = output_env_conv_time(oenv, t);
- teller2++;
+ teller2++;
+ }
+ frame2++;
if (teller2 >= maxframe2)
{
maxframe2 += NFRAME;
real *w_rls = nullptr;
int m, i, d, frame, outframe, natoms, nout, ncent, newstep = 0, model_nr;
#define SKIP 10
- t_topology top;
+ t_topology *top = nullptr;
gmx_conect gc = nullptr;
int ePBC = -1;
t_atoms *atoms = nullptr, useatoms;
if (bTPS)
{
- read_tps_conf(top_file, &top, &ePBC, &xp, nullptr, top_box,
+ snew(top, 1);
+ read_tps_conf(top_file, top, &ePBC, &xp, nullptr, top_box,
bReset || bPBCcomRes);
- std::strncpy(top_title, *top.name, 255);
+ std::strncpy(top_title, *top->name, 255);
top_title[255] = '\0';
- atoms = &top.atoms;
+ atoms = &top->atoms;
- if (0 == top.mols.nr && (bCluster || bPBCcomMol))
+ if (0 == top->mols.nr && (bCluster || bPBCcomMol))
{
gmx_fatal(FARGS, "Option -pbc %s requires a .tpr file for the -s option", pbc_opt[pbc_enum]);
}
if (bCONECT)
{
- gc = gmx_conect_generate(&top);
+ gc = gmx_conect_generate(top);
}
if (bRmPBC)
{
- gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
+ gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
}
}
store original location (to put structure back) */
if (bRmPBC)
{
- gmx_rmpbc(gpbc, top.atoms.nr, top_box, xp);
+ gmx_rmpbc(gpbc, top->atoms.nr, top_box, xp);
}
copy_rvec(xp[index[0]], x_shift);
reset_x_ndim(nfitdim, ifit, ind_fit, atoms->nr, nullptr, xp, w_rls);
}
else if (bCluster)
{
- calc_pbc_cluster(ecenter, ifit, &top, ePBC, fr.x, ind_fit, fr.box);
+ calc_pbc_cluster(ecenter, ifit, top, ePBC, fr.x, ind_fit, fr.box);
}
if (bPFit)
if (bPBCcomMol)
{
put_molecule_com_in_box(unitcell_enum, ecenter,
- &top.mols,
+ &top->mols,
natoms, atoms->atom, ePBC, fr.box, fr.x);
}
/* Copy the input trxframe struct to the output trxframe struct */
if (bTPS)
{
- done_top(&top);
+ done_top(top);
+ sfree(top);
}
sfree(xp);
sfree(xmem);
TrjconvWithIndexGroupSubset,
::testing::ValuesIn(trajectoryFileNames));
+class TrjconvWithoutTopologyFile : public gmx::test::CommandLineTestBase,
+ public ::testing::WithParamInterface<const char *>
+{
+ public:
+ void runTest(const char *fileName)
+ {
+ auto &cmdline = commandLine();
+
+ setInputFile("-f", fileName);
+ setInputFile("-n", "spc2.ndx");
+ setOutputFile("-o", "spc-traj.trr", gmx::test::NoTextMatch());
+
+ gmx::test::StdioTestHelper stdioHelper(&fileManager());
+ stdioHelper.redirectStringToStdin("SecondWaterMolecule\n");
+
+ /* As mentioned above, the tests don't check much besides
+ * that trjconv does not crash.
+ */
+ ASSERT_EQ(0, gmx_trjconv(cmdline.argc(), cmdline.argv()));
+ }
+};
+
+TEST_P(TrjconvWithoutTopologyFile, WithDifferentInputFormats)
+{
+ runTest(GetParam());
+}
+
+INSTANTIATE_TEST_CASE_P(NoFatalErrorWhenWritingFrom,
+ TrjconvWithoutTopologyFile,
+ ::testing::ValuesIn(trajectoryFileNames));
} // namespace
if (ir->bPull)
{
- pull = set_pull_init(ir, &sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS]);
+ pull = set_pull_init(ir, &sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS], wi);
}
/* Modules that supply external potential for pull coordinates
/* Process the pull coordinates after reading the pull groups */
pull_t *set_pull_init(t_inputrec *ir, const gmx_mtop_t *mtop,
- rvec *x, matrix box, real lambda);
+ rvec *x, matrix box, real lambda,
+ warninp_t wi);
/* Prints the initial pull group distances in x.
* If requested, adds the current distance to the initial reference location.
* Returns the pull_t pull work struct. This should be passed to finish_pull()
}
pull_t *set_pull_init(t_inputrec *ir, const gmx_mtop_t *mtop,
- rvec *x, matrix box, real lambda)
+ rvec *x, matrix box, real lambda,
+ warninp_t wi)
{
pull_params_t *pull;
pull_t *pull_work;
pull_calc_coms(nullptr, pull_work, md, &pbc, t_start, x, nullptr);
+ int groupThatFailsPbc = pullCheckPbcWithinGroups(*pull_work, x, pbc);
+ if (groupThatFailsPbc >= 0)
+ {
+ char buf[STRLEN];
+ sprintf(buf,
+ "Pull group %d has atoms at a distance larger than %g times half the box size from the PBC atom (%d). If atoms are or will more beyond half the box size from the PBC atom, the COM will be ill defined.",
+ groupThatFailsPbc,
+ c_pullGroupPbcMargin,
+ pull->group[groupThatFailsPbc].pbcatom);
+ set_warning_line(wi, nullptr, -1);
+ warning(wi, buf);
+ }
+
fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
for (c = 0; c < pull->ncoord; c++)
{
int res0;
} t_moltypes;
-static void sort_molecule(t_atoms **atoms_solvt, std::vector<RVec> *x,
+static void sort_molecule(t_atoms **atoms_solvt,
+ t_atoms **newatoms,
+ std::vector<RVec> *x,
std::vector<RVec> *v)
{
int atnr, i, j, moltp = 0, nrmoltypes, resi_o, resi_n, resnr;
t_moltypes *moltypes;
- t_atoms *atoms, *newatoms;
+ t_atoms *atoms;
fprintf(stderr, "Sorting configuration\n");
}
/* now put them there: */
- snew(newatoms, 1);
- init_t_atoms(newatoms, atoms->nr, FALSE);
- newatoms->nres = atoms->nres;
- snew(newatoms->resinfo, atoms->nres);
+ snew(*newatoms, 1);
+ init_t_atoms(*newatoms, atoms->nr, FALSE);
+ (*newatoms)->nres = atoms->nres;
+ srenew((*newatoms)->resinfo, atoms->nres);
std::vector<RVec> newx(x->size());
std::vector<RVec> newv(v->size());
if (strcmp(*atoms->resinfo[resi_o].name, moltypes[moltp].name) == 0)
{
/* Copy the residue info */
- newatoms->resinfo[resi_n] = atoms->resinfo[resi_o];
- newatoms->resinfo[resi_n].nr = resnr;
+ (*newatoms)->resinfo[resi_n] = atoms->resinfo[resi_o];
+ (*newatoms)->resinfo[resi_n].nr = resnr;
/* Copy the atom info */
do
{
- newatoms->atom[j] = atoms->atom[i];
- newatoms->atomname[j] = atoms->atomname[i];
- newatoms->atom[j].resind = resi_n;
+ (*newatoms)->atom[j] = atoms->atom[i];
+ (*newatoms)->atomname[j] = atoms->atomname[i];
+ (*newatoms)->atom[j].resind = resi_n;
copy_rvec((*x)[i], newx[j]);
if (!v->empty())
{
/* put them back into the original arrays and throw away temporary arrays */
done_atom(atoms);
- *atoms_solvt = newatoms;
+ *atoms_solvt = (*newatoms);
std::swap(*x, newx);
std::swap(*v, newv);
}
}
/* Sort the solvent mixture, not the protein... */
- sort_molecule(&atoms_solvt, &x_solvt, &v_solvt);
+ t_atoms *newatoms = nullptr;
+ sort_molecule(&atoms_solvt, &newatoms, &x_solvt, &v_solvt);
// Merge the two configurations.
x->insert(x->end(), x_solvt.begin(), x_solvt.end());
done_top(top_solvt);
sfree(top_solvt);
+ if (newatoms)
+ {
+ done_atom(newatoms);
+ sfree(newatoms);
+ }
}
-static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
+static void update_top(t_atoms *atoms, int firstSolventResidueIndex, matrix box, int NFILE, t_filenm fnm[],
gmx_atomprop_t aps)
{
- FILE *fpin, *fpout;
- char buf[STRLEN], buf2[STRLEN], *temp;
- const char *topinout;
- int line;
- bool bSystem, bMolecules, bSkip;
- int i, nsol = 0;
- double mtot;
- real vol, mm;
-
- for (i = 0; (i < atoms->nres); i++)
- {
- /* calculate number of SOLvent molecules */
- if ( (strcmp(*atoms->resinfo[i].name, "SOL") == 0) ||
- (strcmp(*atoms->resinfo[i].name, "WAT") == 0) ||
- (strcmp(*atoms->resinfo[i].name, "HOH") == 0) )
- {
- nsol++;
- }
- }
+ FILE *fpin, *fpout;
+ char buf[STRLEN], buf2[STRLEN], *temp;
+ const char *topinout;
+ int line;
+ bool bSystem;
+ int i;
+ double mtot;
+ real vol, mm;
+
+ int nsol = atoms->nres - firstSolventResidueIndex;
+
mtot = 0;
for (i = 0; (i < atoms->nr); i++)
{
fprintf(stderr, "Volume : %10g (nm^3)\n", vol);
fprintf(stderr, "Density : %10g (g/l)\n",
(mtot*1e24)/(AVOGADRO*vol));
- fprintf(stderr, "Number of SOL molecules: %5d \n\n", nsol);
+ fprintf(stderr, "Number of solvent molecules: %5d \n\n", nsol);
/* open topology file and append sol molecules */
topinout = ftp2fn(efTOP, NFILE, fnm);
fpin = gmx_ffopen(topinout, "r");
fpout = gmx_fopen_temporary(temporary_filename);
line = 0;
- bSystem = bMolecules = false;
+ bSystem = false;
while (fgets(buf, STRLEN, fpin))
{
- bSkip = false;
line++;
strcpy(buf2, buf);
if ((temp = strchr(buf2, '\n')) != nullptr)
ltrim(buf2);
rtrim(buf2);
bSystem = (gmx_strcasecmp(buf2, "system") == 0);
- bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
}
}
else if (bSystem && nsol && (buf[0] != ';') )
bSystem = false;
}
}
- else if (bMolecules)
+ fprintf(fpout, "%s", buf);
+ }
+ gmx_ffclose(fpin);
+
+ // Add new solvent molecules to the topology
+ if (nsol > 0)
+ {
+ std::string currRes = *atoms->resinfo[firstSolventResidueIndex].name;
+ int resCount = 0;
+
+ // Iterate through solvent molecules and increment a count until new resname found
+ for (int i = firstSolventResidueIndex; i < atoms->nres; i++)
{
- /* check if this is a line with solvent molecules */
- sscanf(buf, "%4095s", buf2);
- if (strcmp(buf2, "SOL") == 0)
+ if ((currRes == *atoms->resinfo[i].name))
{
- sscanf(buf, "%*4095s %20d", &i);
- nsol -= i;
- if (nsol < 0)
- {
- bSkip = true;
- nsol += i;
- }
+ resCount += 1;
}
- }
- if (bSkip)
- {
- if ((temp = strchr(buf, '\n')) != nullptr)
+ else
{
- temp[0] = '\0';
+ // Change topology and restart count
+ fprintf(stdout, "Adding line for %d solvent molecules with resname (%s) to "
+ "topology file (%s)\n", resCount, currRes.c_str(), topinout);
+ fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
+ currRes = *atoms->resinfo[i].name;
+ resCount = 1;
}
- fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
- line, buf, topinout);
- }
- else
- {
- fprintf(fpout, "%s", buf);
}
- }
- gmx_ffclose(fpin);
- if (nsol)
- {
- fprintf(stdout, "Adding line for %d solvent molecules to "
- "topology file (%s)\n", nsol, topinout);
- fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
+ // One more print needed for last residue type
+ fprintf(stdout, "Adding line for %d solvent molecules with resname (%s) to "
+ "topology file (%s)\n", resCount, currRes.c_str(), topinout);
+ fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
}
gmx_ffclose(fpout);
make_backup(topinout);
};
#define NFILE asize(fnm)
- real defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
- rvec new_box = {0.0, 0.0, 0.0};
- gmx_bool bReadV = FALSE;
- int max_sol = 0;
+ real defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
+ rvec new_box = {0.0, 0.0, 0.0};
+ gmx_bool bReadV = FALSE;
+ int max_sol = 0;
+ int firstSolventResidueIndex = 0;
gmx_output_env_t *oenv;
t_pargs pa[] = {
{ "-box", FALSE, etRVEC, {new_box},
fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
bProt = FALSE;
}
+ else
+ {
+ firstSolventResidueIndex = top->atoms.nres;
+ }
}
if (bBox)
{
/* print size of generated configuration */
fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
top->atoms.nr, top->atoms.nres);
- update_top(&top->atoms, box, NFILE, fnm, aps);
+ update_top(&top->atoms, firstSolventResidueIndex, box, NFILE, fnm, aps);
gmx_atomprop_destroy(aps);
done_top(top);
--- /dev/null
+216H2O,WATJP01,SPC216,SPC-MODEL,300K,BOX(M)=1.86206NM,WFVG,MAR. 1984
+ 648
+ 1HOH OW 1 .230 .628 .113
+ 1HOH HW1 2 .137 .626 .150
+ 1HOH HW2 3 .231 .589 .021
+ 2HOH OW 4 .225 .275 -.866
+ 2HOH HW1 5 .260 .258 -.774
+ 2HOH HW2 6 .137 .230 -.878
+ 3HOH OW 7 .019 .368 .647
+ 3HOH HW1 8 -.063 .411 .686
+ 3HOH HW2 9 -.009 .295 .584
+ 4HOH OW 10 .569 -.587 -.697
+ 4HOH HW1 11 .476 -.594 -.734
+ 4HOH HW2 12 .580 -.498 -.653
+ 5HOH OW 13 -.307 -.351 .703
+ 5HOH HW1 14 -.364 -.367 .784
+ 5HOH HW2 15 -.366 -.341 .623
+ 6HOH OW 16 -.119 .618 .856
+ 6HOH HW1 17 -.086 .712 .856
+ 6HOH HW2 18 -.068 .564 .922
+ 7HOH OW 19 -.727 .703 .717
+ 7HOH HW1 20 -.670 .781 .692
+ 7HOH HW2 21 -.787 .729 .793
+ 8HOH OW 22 -.107 .607 .231
+ 8HOH HW1 23 -.119 .594 .132
+ 8HOH HW2 24 -.137 .526 .280
+ 9HOH OW 25 .768 -.718 -.839
+ 9HOH HW1 26 .690 -.701 -.779
+ 9HOH HW2 27 .802 -.631 -.875
+ 10HOH OW 28 .850 .798 -.039
+ 10HOH HW1 29 .846 .874 .026
+ 10HOH HW2 30 .872 .834 -.130
+ 11HOH OW 31 .685 -.850 .665
+ 11HOH HW1 32 .754 -.866 .735
+ 11HOH HW2 33 .612 -.793 .703
+ 12HOH OW 34 .686 -.701 -.059
+ 12HOH HW1 35 .746 -.622 -.045
+ 12HOH HW2 36 .600 -.670 -.100
+ 13HOH OW 37 .335 -.427 -.801
+ 13HOH HW1 38 .257 -.458 -.854
+ 13HOH HW2 39 .393 -.369 -.858
+ 14HOH OW 40 -.402 -.357 -.523
+ 14HOH HW1 41 -.378 -.263 -.497
+ 14HOH HW2 42 -.418 -.411 -.441
+ 15HOH OW 43 .438 .392 -.363
+ 15HOH HW1 44 .520 .336 -.354
+ 15HOH HW2 45 .357 .334 -.359
+ 16HOH OW 46 -.259 .447 .737
+ 16HOH HW1 47 -.333 .493 .687
+ 16HOH HW2 48 -.208 .515 .790
+ 17HOH OW 49 .231 -.149 .483
+ 17HOH HW1 50 .265 -.072 .537
+ 17HOH HW2 51 .275 -.149 .393
+ 18HOH OW 52 -.735 -.521 -.172
+ 18HOH HW1 53 -.688 -.521 -.084
+ 18HOH HW2 54 -.783 -.608 -.183
+ 19HOH OW 55 .230 -.428 .538
+ 19HOH HW1 56 .204 -.332 .538
+ 19HOH HW2 57 .159 -.482 .583
+ 20HOH OW 58 .240 -.771 .886
+ 20HOH HW1 59 .254 -.855 .938
+ 20HOH HW2 60 .185 -.707 .941
+ 21HOH OW 61 .620 -.076 -.423
+ 21HOH HW1 62 .528 -.093 -.388
+ 21HOH HW2 63 .648 .016 -.397
+ 22HOH OW 64 .606 -.898 .123
+ 22HOH HW1 65 .613 -.814 .069
+ 22HOH HW2 66 .652 -.885 .211
+ 23HOH OW 67 -.268 .114 -.382
+ 23HOH HW1 68 -.286 .181 -.454
+ 23HOH HW2 69 -.271 .160 -.293
+ 24HOH OW 70 .122 .643 .563
+ 24HOH HW1 71 .077 .555 .580
+ 24HOH HW2 72 .121 .697 .647
+ 25HOH OW 73 -.020 -.095 .359
+ 25HOH HW1 74 .034 -.124 .439
+ 25HOH HW2 75 .010 -.005 .330
+ 26HOH OW 76 .027 -.266 .117
+ 26HOH HW1 77 .008 -.362 .138
+ 26HOH HW2 78 -.006 -.208 .192
+ 27HOH OW 79 -.173 .922 .612
+ 27HOH HW1 80 -.078 .893 .620
+ 27HOH HW2 81 -.181 .987 .537
+ 28HOH OW 82 -.221 -.754 .432
+ 28HOH HW1 83 -.135 -.752 .380
+ 28HOH HW2 84 -.207 -.707 .520
+ 29HOH OW 85 .113 .737 -.265
+ 29HOH HW1 86 .201 .724 -.220
+ 29HOH HW2 87 .100 .834 -.287
+ 30HOH OW 88 .613 -.497 .726
+ 30HOH HW1 89 .564 -.584 .735
+ 30HOH HW2 90 .590 -.454 .639
+ 31HOH OW 91 -.569 -.634 -.439
+ 31HOH HW1 92 -.532 -.707 -.497
+ 31HOH HW2 93 -.517 -.629 -.354
+ 32HOH OW 94 .809 .004 .502
+ 32HOH HW1 95 .849 .095 .493
+ 32HOH HW2 96 .709 .012 .508
+ 33HOH OW 97 .197 -.886 -.598
+ 33HOH HW1 98 .286 -.931 -.612
+ 33HOH HW2 99 .124 -.951 -.617
+ 34HOH OW 100 -.337 -.863 .190
+ 34HOH HW1 101 -.400 -.939 .203
+ 34HOH HW2 102 -.289 -.845 .276
+ 35HOH OW 103 -.675 -.070 -.246
+ 35HOH HW1 104 -.651 -.010 -.322
+ 35HOH HW2 105 -.668 -.165 -.276
+ 36HOH OW 106 .317 .251 -.061
+ 36HOH HW1 107 .388 .322 -.055
+ 36HOH HW2 108 .229 .290 -.033
+ 37HOH OW 109 -.396 -.445 -.909
+ 37HOH HW1 110 -.455 -.439 -.829
+ 37HOH HW2 111 -.411 -.533 -.955
+ 38HOH OW 112 -.195 -.148 .572
+ 38HOH HW1 113 -.236 -.171 .484
+ 38HOH HW2 114 -.213 -.222 .637
+ 39HOH OW 115 .598 .729 .270
+ 39HOH HW1 116 .622 .798 .202
+ 39HOH HW2 117 .520 .762 .324
+ 40HOH OW 118 -.581 .345 -.918
+ 40HOH HW1 119 -.667 .295 -.931
+ 40HOH HW2 120 -.519 .291 -.862
+ 41HOH OW 121 -.286 -.200 .307
+ 41HOH HW1 122 -.197 -.154 .310
+ 41HOH HW2 123 -.307 -.224 .212
+ 42HOH OW 124 .807 .605 -.397
+ 42HOH HW1 125 .760 .602 -.308
+ 42HOH HW2 126 .756 .550 -.463
+ 43HOH OW 127 -.468 .469 -.188
+ 43HOH HW1 128 -.488 .512 -.100
+ 43HOH HW2 129 -.390 .407 -.179
+ 44HOH OW 130 -.889 .890 -.290
+ 44HOH HW1 131 -.843 .806 -.319
+ 44HOH HW2 132 -.945 .924 -.365
+ 45HOH OW 133 -.871 .410 -.620
+ 45HOH HW1 134 -.948 .444 -.566
+ 45HOH HW2 135 -.905 .359 -.699
+ 46HOH OW 136 -.821 .701 .429
+ 46HOH HW1 137 -.795 .697 .525
+ 46HOH HW2 138 -.906 .650 .415
+ 47HOH OW 139 .076 .811 .789
+ 47HOH HW1 140 .175 .799 .798
+ 47HOH HW2 141 .052 .906 .810
+ 48HOH OW 142 .130 -.041 -.291
+ 48HOH HW1 143 .120 -.056 -.192
+ 48HOH HW2 144 .044 -.005 -.327
+ 49HOH OW 145 .865 .348 .195
+ 49HOH HW1 146 .924 .411 .146
+ 49HOH HW2 147 .884 .254 .166
+ 50HOH OW 148 -.143 .585 -.031
+ 50HOH HW1 149 -.169 .674 -.067
+ 50HOH HW2 150 -.145 .517 -.104
+ 51HOH OW 151 -.500 -.718 .545
+ 51HOH HW1 152 -.417 -.747 .497
+ 51HOH HW2 153 -.549 -.651 .489
+ 52HOH OW 154 .550 .196 .885
+ 52HOH HW1 155 .545 .191 .985
+ 52HOH HW2 156 .552 .292 .856
+ 53HOH OW 157 -.854 -.406 .477
+ 53HOH HW1 158 -.900 -.334 .425
+ 53HOH HW2 159 -.858 -.386 .575
+ 54HOH OW 160 .351 -.061 .853
+ 54HOH HW1 161 .401 -.147 .859
+ 54HOH HW2 162 .416 .016 .850
+ 55HOH OW 163 -.067 -.796 .873
+ 55HOH HW1 164 -.129 -.811 .797
+ 55HOH HW2 165 -.119 -.785 .958
+ 56HOH OW 166 -.635 -.312 -.356
+ 56HOH HW1 167 -.629 -.389 -.292
+ 56HOH HW2 168 -.687 -.338 -.436
+ 57HOH OW 169 .321 -.919 .242
+ 57HOH HW1 170 .403 -.880 .200
+ 57HOH HW2 171 .294 -1.001 .193
+ 58HOH OW 172 -.404 .735 .728
+ 58HOH HW1 173 -.409 .670 .803
+ 58HOH HW2 174 -.324 .794 .741
+ 59HOH OW 175 .461 -.596 -.135
+ 59HOH HW1 176 .411 -.595 -.221
+ 59HOH HW2 177 .398 -.614 -.059
+ 60HOH OW 178 -.751 -.086 .237
+ 60HOH HW1 179 -.811 -.148 .287
+ 60HOH HW2 180 -.720 -.130 .152
+ 61HOH OW 181 .202 .285 -.364
+ 61HOH HW1 182 .122 .345 -.377
+ 61HOH HW2 183 .192 .236 -.278
+ 62HOH OW 184 -.230 -.485 .081
+ 62HOH HW1 185 -.262 -.391 .071
+ 62HOH HW2 186 -.306 -.548 .069
+ 63HOH OW 187 .464 -.119 .323
+ 63HOH HW1 188 .497 -.080 .409
+ 63HOH HW2 189 .540 -.126 .258
+ 64HOH OW 190 -.462 .107 .426
+ 64HOH HW1 191 -.486 .070 .336
+ 64HOH HW2 192 -.363 .123 .430
+ 65HOH OW 193 .249 -.077 -.621
+ 65HOH HW1 194 .306 -.142 -.571
+ 65HOH HW2 195 .233 -.110 -.714
+ 66HOH OW 196 -.922 -.164 .904
+ 66HOH HW1 197 -.842 -.221 .925
+ 66HOH HW2 198 -.971 -.204 .827
+ 67HOH OW 199 .382 .700 .480
+ 67HOH HW1 200 .427 .610 .477
+ 67HOH HW2 201 .288 .689 .513
+ 68HOH OW 202 -.315 .222 -.133
+ 68HOH HW1 203 -.320 .259 -.041
+ 68HOH HW2 204 -.387 .153 -.145
+ 69HOH OW 205 .614 .122 .117
+ 69HOH HW1 206 .712 .100 .124
+ 69HOH HW2 207 .583 .105 .024
+ 70HOH OW 208 .781 .264 -.113
+ 70HOH HW1 209 .848 .203 -.070
+ 70HOH HW2 210 .708 .283 -.048
+ 71HOH OW 211 .888 -.348 -.667
+ 71HOH HW1 212 .865 -.373 -.761
+ 71HOH HW2 213 .949 -.417 -.628
+ 72HOH OW 214 -.511 .590 -.429
+ 72HOH HW1 215 -.483 .547 -.344
+ 72HOH HW2 216 -.486 .686 -.428
+ 73HOH OW 217 .803 -.460 .924
+ 73HOH HW1 218 .893 -.446 .882
+ 73HOH HW2 219 .732 -.458 .853
+ 74HOH OW 220 .922 .503 .899
+ 74HOH HW1 221 .897 .494 .803
+ 74HOH HW2 222 .970 .421 .930
+ 75HOH OW 223 .539 .064 .512
+ 75HOH HW1 224 .458 .065 .570
+ 75HOH HW2 225 .542 .147 .457
+ 76HOH OW 226 -.428 -.674 .041
+ 76HOH HW1 227 -.396 -.750 .098
+ 76HOH HW2 228 -.520 -.647 .071
+ 77HOH OW 229 .297 .035 .171
+ 77HOH HW1 230 .346 .119 .150
+ 77HOH HW2 231 .359 -.030 .216
+ 78HOH OW 232 -.927 .236 .480
+ 78HOH HW1 233 -.975 .277 .402
+ 78HOH HW2 234 -.828 .234 .461
+ 79HOH OW 235 -.786 .683 -.398
+ 79HOH HW1 236 -.866 .622 -.395
+ 79HOH HW2 237 -.705 .630 -.422
+ 80HOH OW 238 -.635 -.292 .793
+ 80HOH HW1 239 -.614 -.218 .728
+ 80HOH HW2 240 -.567 -.292 .866
+ 81HOH OW 241 .459 -.710 .741
+ 81HOH HW1 242 .388 -.737 .806
+ 81HOH HW2 243 .433 -.738 .648
+ 82HOH OW 244 -.591 -.065 .591
+ 82HOH HW1 245 -.547 -.001 .527
+ 82HOH HW2 246 -.641 -.013 .661
+ 83HOH OW 247 -.830 .549 .016
+ 83HOH HW1 248 -.871 .631 -.023
+ 83HOH HW2 249 -.766 .575 .089
+ 84HOH OW 250 .078 .556 -.476
+ 84HOH HW1 251 .170 .555 -.517
+ 84HOH HW2 252 .072 .630 -.409
+ 85HOH OW 253 .561 .222 -.715
+ 85HOH HW1 254 .599 .138 -.678
+ 85HOH HW2 255 .473 .241 -.671
+ 86HOH OW 256 .866 .454 .642
+ 86HOH HW1 257 .834 .526 .580
+ 86HOH HW2 258 .890 .373 .589
+ 87HOH OW 259 -.845 .039 .753
+ 87HOH HW1 260 -.917 .044 .684
+ 87HOH HW2 261 -.869 -.030 .822
+ 88HOH OW 262 -.433 -.689 .867
+ 88HOH HW1 263 -.488 -.773 .860
+ 88HOH HW2 264 -.407 -.660 .775
+ 89HOH OW 265 -.396 .590 -.870
+ 89HOH HW1 266 -.426 .495 -.863
+ 89HOH HW2 267 -.323 .606 -.804
+ 90HOH OW 268 -.005 .833 .377
+ 90HOH HW1 269 .037 .769 .441
+ 90HOH HW2 270 -.043 .782 .299
+ 91HOH OW 271 .488 -.477 .174
+ 91HOH HW1 272 .401 -.492 .221
+ 91HOH HW2 273 .471 -.451 .079
+ 92HOH OW 274 -.198 -.582 .657
+ 92HOH HW1 275 -.099 -.574 .671
+ 92HOH HW2 276 -.243 -.498 .688
+ 93HOH OW 277 -.472 .575 .078
+ 93HOH HW1 278 -.526 .554 .159
+ 93HOH HW2 279 -.381 .534 .087
+ 94HOH OW 280 .527 .256 .328
+ 94HOH HW1 281 .554 .197 .253
+ 94HOH HW2 282 .527 .351 .297
+ 95HOH OW 283 -.108 -.639 -.274
+ 95HOH HW1 284 -.017 -.678 -.287
+ 95HOH HW2 285 -.100 -.543 -.250
+ 96HOH OW 286 -.798 -.515 -.522
+ 96HOH HW1 287 -.878 -.538 -.467
+ 96HOH HW2 288 -.715 -.541 -.473
+ 97HOH OW 289 -.270 -.233 -.237
+ 97HOH HW1 290 -.243 -.199 -.327
+ 97HOH HW2 291 -.191 -.271 -.191
+ 98HOH OW 292 -.751 -.667 -.762
+ 98HOH HW1 293 -.791 -.623 -.681
+ 98HOH HW2 294 -.792 -.630 -.845
+ 99HOH OW 295 -.224 -.763 -.783
+ 99HOH HW1 296 -.219 -.682 -.724
+ 99HOH HW2 297 -.310 -.761 -.834
+ 100HOH OW 298 .915 .089 -.460
+ 100HOH HW1 299 .940 .069 -.555
+ 100HOH HW2 300 .987 .145 -.418
+ 101SOL OW 301 -.882 -.746 -.143
+ 101SOL HW1 302 -.981 -.740 -.133
+ 101SOL HW2 303 -.859 -.826 -.199
+ 102SOL OW 304 .705 -.812 .368
+ 102SOL HW1 305 .691 -.805 .467
+ 102SOL HW2 306 .789 -.863 .350
+ 103SOL OW 307 .410 .813 -.611
+ 103SOL HW1 308 .496 .825 -.561
+ 103SOL HW2 309 .368 .726 -.584
+ 104SOL OW 310 -.588 .386 -.600
+ 104SOL HW1 311 -.567 .460 -.536
+ 104SOL HW2 312 -.677 .403 -.643
+ 105SOL OW 313 .064 -.298 -.531
+ 105SOL HW1 314 .018 -.216 -.565
+ 105SOL HW2 315 .162 -.279 -.522
+ 106SOL OW 316 .367 -.762 .501
+ 106SOL HW1 317 .360 -.679 .445
+ 106SOL HW2 318 .371 -.842 .441
+ 107SOL OW 319 .566 .537 .865
+ 107SOL HW1 320 .578 .603 .791
+ 107SOL HW2 321 .612 .571 .948
+ 108SOL OW 322 -.610 -.514 .388
+ 108SOL HW1 323 -.560 -.437 .428
+ 108SOL HW2 324 -.705 -.512 .420
+ 109SOL OW 325 -.590 -.417 -.720
+ 109SOL HW1 326 -.543 -.404 -.633
+ 109SOL HW2 327 -.656 -.491 -.711
+ 110SOL OW 328 -.280 .639 .472
+ 110SOL HW1 329 -.311 .700 .545
+ 110SOL HW2 330 -.230 .691 .403
+ 111SOL OW 331 .354 -.352 -.533
+ 111SOL HW1 332 .333 -.396 -.620
+ 111SOL HW2 333 .451 -.326 -.530
+ 112SOL OW 334 .402 .751 -.264
+ 112SOL HW1 335 .470 .806 -.311
+ 112SOL HW2 336 .442 .663 -.237
+ 113SOL OW 337 -.275 .779 -.192
+ 113SOL HW1 338 -.367 .817 -.197
+ 113SOL HW2 339 -.215 .826 -.257
+ 114SOL OW 340 -.849 .105 -.092
+ 114SOL HW1 341 -.843 .190 -.144
+ 114SOL HW2 342 -.817 .029 -.149
+ 115SOL OW 343 .504 .050 -.122
+ 115SOL HW1 344 .462 -.007 -.192
+ 115SOL HW2 345 .438 .119 -.090
+ 116SOL OW 346 .573 .870 -.833
+ 116SOL HW1 347 .617 .959 -.842
+ 116SOL HW2 348 .510 .870 -.756
+ 117SOL OW 349 -.502 .862 -.817
+ 117SOL HW1 350 -.577 .862 -.883
+ 117SOL HW2 351 -.465 .770 -.808
+ 118SOL OW 352 -.653 .525 .275
+ 118SOL HW1 353 -.640 .441 .329
+ 118SOL HW2 354 -.682 .599 .335
+ 119SOL OW 355 .307 .213 -.631
+ 119SOL HW1 356 .284 .250 -.541
+ 119SOL HW2 357 .277 .118 -.637
+ 120SOL OW 358 .037 -.552 -.580
+ 120SOL HW1 359 .090 -.601 -.512
+ 120SOL HW2 360 .059 -.454 -.575
+ 121SOL OW 361 .732 .634 -.798
+ 121SOL HW1 362 .791 .608 -.874
+ 121SOL HW2 363 .704 .730 -.809
+ 122SOL OW 364 -.134 -.927 -.008
+ 122SOL HW1 365 -.180 -.934 -.097
+ 122SOL HW2 366 -.196 -.883 .058
+ 123SOL OW 367 .307 .063 .618
+ 123SOL HW1 368 .296 .157 .651
+ 123SOL HW2 369 .302 -.000 .695
+ 124SOL OW 370 -.240 .367 .374
+ 124SOL HW1 371 -.238 .291 .438
+ 124SOL HW2 372 -.288 .444 .414
+ 125SOL OW 373 -.839 .766 -.896
+ 125SOL HW1 374 -.824 .787 -.800
+ 125SOL HW2 375 -.869 .671 -.905
+ 126SOL OW 376 -.882 -.289 -.162
+ 126SOL HW1 377 -.902 -.245 -.250
+ 126SOL HW2 378 -.843 -.380 -.178
+ 127SOL OW 379 -.003 -.344 -.257
+ 127SOL HW1 380 .011 -.317 -.352
+ 127SOL HW2 381 .080 -.322 -.204
+ 128SOL OW 382 .350 .898 -.058
+ 128SOL HW1 383 .426 .942 -.010
+ 128SOL HW2 384 .385 .851 -.140
+ 129SOL OW 385 -.322 .274 .125
+ 129SOL HW1 386 -.383 .199 .148
+ 129SOL HW2 387 -.300 .326 .208
+ 130SOL OW 388 -.559 .838 .042
+ 130SOL HW1 389 -.525 .745 .057
+ 130SOL HW2 390 -.541 .865 -.053
+ 131SOL OW 391 -.794 -.529 .849
+ 131SOL HW1 392 -.787 -.613 .794
+ 131SOL HW2 393 -.732 -.460 .813
+ 132SOL OW 394 .319 .810 -.913
+ 132SOL HW1 395 .412 .846 -.908
+ 132SOL HW2 396 .313 .725 -.861
+ 133SOL OW 397 .339 .509 -.856
+ 133SOL HW1 398 .287 .426 -.873
+ 133SOL HW2 399 .416 .514 -.920
+ 134SOL OW 400 .511 .415 -.054
+ 134SOL HW1 401 .493 .460 .034
+ 134SOL HW2 402 .553 .480 -.117
+ 135SOL OW 403 -.724 .380 -.184
+ 135SOL HW1 404 -.769 .443 -.120
+ 135SOL HW2 405 -.631 .411 -.201
+ 136SOL OW 406 -.702 .207 -.385
+ 136SOL HW1 407 -.702 .271 -.308
+ 136SOL HW2 408 -.674 .255 -.468
+ 137SOL OW 409 .008 -.536 .200
+ 137SOL HW1 410 -.085 -.515 .169
+ 137SOL HW2 411 .018 -.635 .213
+ 138SOL OW 412 .088 -.061 .927
+ 138SOL HW1 413 .046 -.147 .900
+ 138SOL HW2 414 .182 -.058 .893
+ 139SOL OW 415 .504 -.294 .910
+ 139SOL HW1 416 .570 -.220 .919
+ 139SOL HW2 417 .548 -.373 .868
+ 140SOL OW 418 -.860 .796 -.624
+ 140SOL HW1 419 -.819 .764 -.538
+ 140SOL HW2 420 -.956 .769 -.627
+ 141SOL OW 421 .040 .544 -.748
+ 141SOL HW1 422 .125 .511 -.789
+ 141SOL HW2 423 .053 .559 -.650
+ 142SOL OW 424 .189 .520 -.140
+ 142SOL HW1 425 .248 .480 -.210
+ 142SOL HW2 426 .131 .591 -.181
+ 143SOL OW 427 -.493 -.912 -.202
+ 143SOL HW1 428 -.454 -.823 -.182
+ 143SOL HW2 429 -.483 -.932 -.299
+ 144SOL OW 430 .815 .572 .325
+ 144SOL HW1 431 .822 .483 .279
+ 144SOL HW2 432 .721 .606 .317
+ 145SOL OW 433 -.205 .604 -.656
+ 145SOL HW1 434 -.243 .535 -.594
+ 145SOL HW2 435 -.123 .568 -.700
+ 146SOL OW 436 .252 -.298 -.118
+ 146SOL HW1 437 .222 -.241 -.042
+ 146SOL HW2 438 .245 -.395 -.092
+ 147SOL OW 439 .671 .464 -.593
+ 147SOL HW1 440 .637 .375 -.623
+ 147SOL HW2 441 .697 .518 -.673
+ 148SOL OW 442 .930 -.184 -.397
+ 148SOL HW1 443 .906 -.202 -.492
+ 148SOL HW2 444 .960 -.090 -.387
+ 149SOL OW 445 .473 .500 .191
+ 149SOL HW1 446 .534 .580 .195
+ 149SOL HW2 447 .378 .531 .198
+ 150SOL OW 448 .159 -.725 -.396
+ 150SOL HW1 449 .181 -.786 -.320
+ 150SOL HW2 450 .169 -.774 -.482
+ 151SOL OW 451 -.515 -.803 -.628
+ 151SOL HW1 452 -.491 -.866 -.702
+ 151SOL HW2 453 -.605 -.763 -.646
+ 152SOL OW 454 -.560 .855 .309
+ 152SOL HW1 455 -.646 .824 .351
+ 152SOL HW2 456 -.564 .841 .210
+ 153SOL OW 457 -.103 -.115 -.708
+ 153SOL HW1 458 -.042 -.085 -.781
+ 153SOL HW2 459 -.141 -.204 -.730
+ 154SOL OW 460 -.610 -.131 -.734
+ 154SOL HW1 461 -.526 -.126 -.788
+ 154SOL HW2 462 -.633 -.227 -.716
+ 155SOL OW 463 .083 -.604 -.840
+ 155SOL HW1 464 .078 -.605 -.740
+ 155SOL HW2 465 .000 -.645 -.878
+ 156SOL OW 466 .688 -.200 -.146
+ 156SOL HW1 467 .632 -.119 -.137
+ 156SOL HW2 468 .740 -.196 -.232
+ 157SOL OW 469 .903 .086 .133
+ 157SOL HW1 470 .954 .087 .047
+ 157SOL HW2 471 .959 .044 .204
+ 158SOL OW 472 -.136 .135 .523
+ 158SOL HW1 473 -.063 .118 .456
+ 158SOL HW2 474 -.167 .048 .561
+ 159SOL OW 475 -.474 -.289 .477
+ 159SOL HW1 476 -.407 -.277 .403
+ 159SOL HW2 477 -.514 -.200 .500
+ 160SOL OW 478 .130 -.068 -.011
+ 160SOL HW1 479 .089 -.142 .042
+ 160SOL HW2 480 .194 -.017 .047
+ 161SOL OW 481 -.582 .927 .672
+ 161SOL HW1 482 -.522 .846 .674
+ 161SOL HW2 483 -.542 .996 .612
+ 162SOL OW 484 .830 -.589 -.440
+ 162SOL HW1 485 .825 -.556 -.345
+ 162SOL HW2 486 .744 -.570 -.486
+ 163SOL OW 487 .672 -.246 .154
+ 163SOL HW1 488 .681 -.236 .055
+ 163SOL HW2 489 .632 -.335 .175
+ 164SOL OW 490 -.212 -.142 -.468
+ 164SOL HW1 491 -.159 -.132 -.552
+ 164SOL HW2 492 -.239 -.052 -.434
+ 165SOL OW 493 -.021 .175 -.899
+ 165SOL HW1 494 .018 .090 -.935
+ 165SOL HW2 495 -.119 .177 -.918
+ 166SOL OW 496 .263 .326 .720
+ 166SOL HW1 497 .184 .377 .686
+ 166SOL HW2 498 .254 .311 .818
+ 167SOL OW 499 -.668 -.250 .031
+ 167SOL HW1 500 -.662 -.343 .068
+ 167SOL HW2 501 -.727 -.250 -.049
+ 168SOL OW 502 .822 -.860 -.490
+ 168SOL HW1 503 .862 -.861 -.582
+ 168SOL HW2 504 .832 -.768 -.450
+ 169SOL OW 505 .916 .910 .291
+ 169SOL HW1 506 .979 .948 .223
+ 169SOL HW2 507 .956 .827 .330
+ 170SOL OW 508 -.358 -.255 .044
+ 170SOL HW1 509 -.450 -.218 .051
+ 170SOL HW2 510 -.320 -.235 -.046
+ 171SOL OW 511 .372 -.574 -.372
+ 171SOL HW1 512 .359 -.481 -.406
+ 171SOL HW2 513 .288 -.626 -.385
+ 172SOL OW 514 -.248 -.570 -.573
+ 172SOL HW1 515 -.188 -.567 -.493
+ 172SOL HW2 516 -.323 -.506 -.560
+ 173SOL OW 517 -.823 -.764 .696
+ 173SOL HW1 518 -.893 -.811 .750
+ 173SOL HW2 519 -.764 -.832 .653
+ 174SOL OW 520 -.848 .236 -.891
+ 174SOL HW1 521 -.856 .200 -.984
+ 174SOL HW2 522 -.850 .160 -.826
+ 175SOL OW 523 .590 -.375 .491
+ 175SOL HW1 524 .632 -.433 .421
+ 175SOL HW2 525 .546 -.296 .447
+ 176SOL OW 526 -.153 .385 -.481
+ 176SOL HW1 527 -.080 .454 -.477
+ 176SOL HW2 528 -.125 .310 -.540
+ 177SOL OW 529 .255 -.514 .290
+ 177SOL HW1 530 .159 -.513 .263
+ 177SOL HW2 531 .267 -.461 .374
+ 178SOL OW 532 .105 -.849 -.136
+ 178SOL HW1 533 .028 -.882 -.082
+ 178SOL HW2 534 .190 -.879 -.094
+ 179SOL OW 535 .672 .203 -.373
+ 179SOL HW1 536 .762 .187 -.413
+ 179SOL HW2 537 .680 .208 -.274
+ 180SOL OW 538 .075 .345 .033
+ 180SOL HW1 539 -.017 .317 .004
+ 180SOL HW2 540 .106 .422 -.023
+ 181SOL OW 541 -.422 .856 -.464
+ 181SOL HW1 542 -.479 .908 -.527
+ 181SOL HW2 543 -.326 .868 -.488
+ 182SOL OW 544 .072 .166 .318
+ 182SOL HW1 545 .055 .249 .264
+ 182SOL HW2 546 .162 .129 .296
+ 183SOL OW 547 -.679 -.527 .119
+ 183SOL HW1 548 -.778 -.538 .121
+ 183SOL HW2 549 -.645 -.512 .212
+ 184SOL OW 550 .613 .842 -.431
+ 184SOL HW1 551 .669 .923 -.448
+ 184SOL HW2 552 .672 .762 -.428
+ 185SOL OW 553 -.369 -.095 -.903
+ 185SOL HW1 554 -.336 -.031 -.972
+ 185SOL HW2 555 -.303 -.101 -.828
+ 186SOL OW 556 .716 .565 -.154
+ 186SOL HW1 557 .735 .630 -.080
+ 186SOL HW2 558 .776 .485 -.145
+ 187SOL OW 559 -.412 -.642 -.229
+ 187SOL HW1 560 -.421 -.652 -.130
+ 187SOL HW2 561 -.316 -.649 -.255
+ 188SOL OW 562 .390 -.121 -.302
+ 188SOL HW1 563 .299 -.080 -.304
+ 188SOL HW2 564 .383 -.215 -.270
+ 189SOL OW 565 -.188 .883 -.608
+ 189SOL HW1 566 -.215 .794 -.645
+ 189SOL HW2 567 -.187 .951 -.681
+ 190SOL OW 568 -.637 .325 .449
+ 190SOL HW1 569 -.572 .251 .438
+ 190SOL HW2 570 -.617 .375 .533
+ 191SOL OW 571 .594 .745 .652
+ 191SOL HW1 572 .644 .830 .633
+ 191SOL HW2 573 .506 .747 .604
+ 192SOL OW 574 -.085 .342 -.220
+ 192SOL HW1 575 -.102 .373 -.314
+ 192SOL HW2 576 -.169 .305 -.182
+ 193SOL OW 577 -.132 -.928 -.345
+ 193SOL HW1 578 -.094 -.837 -.330
+ 193SOL HW2 579 -.140 -.945 -.444
+ 194SOL OW 580 .859 -.488 .016
+ 194SOL HW1 581 .813 -.473 .104
+ 194SOL HW2 582 .903 -.403 -.014
+ 195SOL OW 583 .661 -.072 -.909
+ 195SOL HW1 584 .615 .016 -.922
+ 195SOL HW2 585 .760 -.060 -.916
+ 196SOL OW 586 -.454 -.011 -.142
+ 196SOL HW1 587 -.550 -.022 -.169
+ 196SOL HW2 588 -.398 -.078 -.190
+ 197SOL OW 589 .859 -.906 .861
+ 197SOL HW1 590 .913 -.975 .909
+ 197SOL HW2 591 .827 -.837 .927
+ 198SOL OW 592 -.779 -.878 .087
+ 198SOL HW1 593 -.802 -.825 .005
+ 198SOL HW2 594 -.698 -.934 .068
+ 199SOL OW 595 -.001 -.293 .851
+ 199SOL HW1 596 -.072 -.305 .781
+ 199SOL HW2 597 .000 -.372 .911
+ 200SOL OW 598 .221 -.548 -.018
+ 200SOL HW1 599 .156 -.621 -.039
+ 200SOL HW2 600 .225 -.534 .080
+ 201SOL OW 601 .079 -.622 .653
+ 201SOL HW1 602 .078 -.669 .741
+ 201SOL HW2 603 .161 -.650 .602
+ 202SOL OW 604 .672 -.471 -.238
+ 202SOL HW1 605 .594 -.521 -.200
+ 202SOL HW2 606 .669 -.376 -.207
+ 203SOL OW 607 -.038 .192 -.635
+ 203SOL HW1 608 -.042 .102 -.591
+ 203SOL HW2 609 -.035 .181 -.734
+ 204SOL OW 610 .428 .424 .520
+ 204SOL HW1 611 .458 .352 .458
+ 204SOL HW2 612 .389 .384 .603
+ 205SOL OW 613 -.157 -.375 -.758
+ 205SOL HW1 614 -.250 -.400 -.785
+ 205SOL HW2 615 -.131 -.425 -.676
+ 206SOL OW 616 .317 .547 -.582
+ 206SOL HW1 617 .355 .488 -.510
+ 206SOL HW2 618 .357 .521 -.670
+ 207SOL OW 619 .812 -.276 .687
+ 207SOL HW1 620 .844 -.266 .593
+ 207SOL HW2 621 .733 -.338 .689
+ 208SOL OW 622 -.438 .214 -.750
+ 208SOL HW1 623 -.386 .149 -.695
+ 208SOL HW2 624 -.487 .277 -.689
+ 209SOL OW 625 -.861 .034 -.708
+ 209SOL HW1 626 -.924 -.038 -.739
+ 209SOL HW2 627 -.768 -.002 -.708
+ 210SOL OW 628 .770 -.532 .301
+ 210SOL HW1 629 .724 -.619 .318
+ 210SOL HW2 630 .861 -.535 .342
+ 211SOL OW 631 .618 -.295 -.578
+ 211SOL HW1 632 .613 -.213 -.521
+ 211SOL HW2 633 .707 -.298 -.623
+ 212SOL OW 634 -.510 .052 .168
+ 212SOL HW1 635 -.475 .011 .084
+ 212SOL HW2 636 -.600 .014 .188
+ 213SOL OW 637 -.562 .453 .691
+ 213SOL HW1 638 -.621 .533 .695
+ 213SOL HW2 639 -.547 .418 .784
+ 214SOL OW 640 -.269 .221 .882
+ 214SOL HW1 641 -.353 .220 .936
+ 214SOL HW2 642 -.267 .304 .826
+ 215SOL OW 643 .039 -.785 .300
+ 215SOL HW1 644 .138 -.796 .291
+ 215SOL HW2 645 -.001 -.871 .332
+ 216SOL OW 646 .875 -.216 .337
+ 216SOL HW1 647 .798 -.251 .283
+ 216SOL HW2 648 .843 -.145 .399
+ 1.86206 1.86206 1.86206
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <GroFile Name="Header">
+ <String Name="Title">Test system for solvate</String>
+ <Int Name="Number of atoms">12141</Int>
+ </GroFile>
+ </File>
+ <File Name="-p">
+ <String Name="Contents"><![CDATA[
+[ system ]
+; Name
+simple system in water
+
+[ molecules ]
+; Compound #mols
+A 1
+B 1
+HOH 1876
+SOL 2169
+]]></String>
+ </File>
+ </OutputFiles>
+</ReferenceData>
\ No newline at end of file
--- /dev/null
+Test system for solvate
+ 6
+ 1A A1 1 1.000 1.000 0.000
+ 1A A2 2 1.000 2.000 0.000
+ 1A A3 3 1.000 3.000 0.000
+ 2B B1 4 1.000 4.000 1.000
+ 2B B2 5 2.000 1.000 2.000
+ 2B B3 6 2.000 2.000 3.000
+ 5.00000 5.00000 5.00000
--- /dev/null
+[ system ]
+; Name
+simple system
+
+[ molecules ]
+; Compound #mols
+A 1
+B 1
using gmx::test::CommandLine;
using gmx::test::ConfMatch;
+using gmx::test::ExactTextMatch;
class SolvateTest : public gmx::test::CommandLineTestBase
{
runTest(CommandLine(cmdline));
}
+TEST_F(SolvateTest, update_Topology_Works)
+{
+ // use solvent box with 2 solvents, check that topology has been updated
+ const char *const cmdline[] = {
+ "solvate"
+ };
+ setInputFile("-cs", "mixed_solvent.gro");
+ setInputFile("-cp", "simple.gro");
+
+ // TODO: Consider adding a convenience method for this.
+ // Copies topology file to where it would be found as an output file, so the copied
+ // .top file is used as both input and output
+ std::string topFileName = gmx::test::TestFileManager::getInputFilePath("simple.top");
+ std::string modifiableTopFileName = fileManager().getTemporaryFilePath("simple.top");
+ gmx_file_copy(topFileName.c_str(), modifiableTopFileName.c_str(), true);
+ setOutputFile("-p", "simple.top", ExactTextMatch());
+
+ runTest(CommandLine(cmdline));
+}
+
} // namespace
}
}
}
- GMX_RELEASE_ASSERT(cudaSuccess == cudaPeekAtLastError(), "We promise to return with clean CUDA state!");
+
+ stat = cudaPeekAtLastError();
+ GMX_RELEASE_ASSERT(stat == cudaSuccess,
+ gmx::formatString("We promise to return with clean CUDA state, but non-success state encountered: %s: %s",
+ cudaGetErrorName(stat), cudaGetErrorString(stat)).c_str());
gpu_info->n_dev = ndev;
gpu_info->gpu_dev = devs;
if (getenv("GMX_OCL_DISABLE_FASTMATH") == nullptr)
{
compilerOptions += " -cl-fast-relaxed-math";
+
+ // Hint to the compiler that it can flush denorms to zero.
+ // In CUDA this is triggered by the -use_fast_math flag, equivalent with
+ // -cl-fast-relaxed-math, hence the inclusion on the conditional block.
+ compilerOptions += " -cl-denorms-are-zero";
}
if ((deviceVendorId == OCL_VENDOR_NVIDIA) && getenv("GMX_OCL_VERBOSE"))
}
}
+/* Returns the mass of atom atomIndex or 1 when setMassesToOne=true */
+static real getMass(const t_atoms &atoms,
+ int atomIndex,
+ bool setMassesToOne)
+{
+ if (!setMassesToOne)
+ {
+ return atoms.atom[atomIndex].m;
+ }
+ else
+ {
+ return 1;
+ }
+}
+
static void get_vsite_masses(const gmx_moltype_t *moltype,
const gmx_ffparams_t *ffparams,
+ bool setMassesToOne,
real *vsite_m,
int *n_nonlin_vsite)
{
{
const t_iparams *ip;
real inv_mass, coeff, m_aj;
- int a1, aj;
+ int a1;
ip = &ffparams->iparams[il->iatoms[i]];
assert(maxj <= 5);
for (j = 1; j < maxj; j++)
{
- cam[j] = moltype->atoms.atom[il->iatoms[i+1+j]].m;
+ int aj = il->iatoms[i + 1 + j];
+ cam[j] = getMass(moltype->atoms, aj, setMassesToOne);
if (cam[j] == 0)
{
- cam[j] = vsite_m[il->iatoms[i+1+j]];
+ cam[j] = vsite_m[aj];
}
if (cam[j] == 0)
{
gmx_fatal(FARGS, "In molecule type '%s' %s construction involves atom %d, which is a virtual site of equal or high complexity. This is not supported.",
*moltype->name,
interaction_function[ft].longname,
- il->iatoms[i+1+j]+1);
+ aj + 1);
}
}
inv_mass = 0;
for (j = 0; j < 3*ffparams->iparams[il->iatoms[i]].vsiten.n; j += 3)
{
- aj = il->iatoms[i+j+2];
- coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
+ int aj = il->iatoms[i + j + 2];
+ coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
if (moltype->atoms.atom[aj].ptype == eptVSite)
{
m_aj = vsite_m[aj];
}
static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
+ bool setMassesToOne,
verletbuf_atomtype_t **att_p,
int *natt_p,
int *n_nonlin_vsite)
ip = &mtop->ffparams.iparams[il->iatoms[i]];
a1 = il->iatoms[i+1];
a2 = il->iatoms[i+2];
- if (atoms->atom[a2].m > prop[a1].con_mass)
+ real mass1 = getMass(*atoms, a1, setMassesToOne);
+ real mass2 = getMass(*atoms, a2, setMassesToOne);
+ if (mass2 > prop[a1].con_mass)
{
- prop[a1].con_mass = atoms->atom[a2].m;
+ prop[a1].con_mass = mass2;
prop[a1].con_len = ip->constr.dA;
}
- if (atoms->atom[a1].m > prop[a2].con_mass)
+ if (mass1 > prop[a2].con_mass)
{
- prop[a2].con_mass = atoms->atom[a1].m;
+ prop[a2].con_mass = mass1;
prop[a2].con_len = ip->constr.dA;
}
}
* If this is not the case, we overestimate the displacement,
* which leads to a larger buffer (ok since this is an exotic case).
*/
- prop[a1].con_mass = atoms->atom[a2].m;
+ prop[a1].con_mass = getMass(*atoms, a2, setMassesToOne);
prop[a1].con_len = ip->settle.doh;
- prop[a2].con_mass = atoms->atom[a1].m;
+ prop[a2].con_mass = getMass(*atoms, a1, setMassesToOne);
prop[a2].con_len = ip->settle.doh;
- prop[a3].con_mass = atoms->atom[a1].m;
+ prop[a3].con_mass = getMass(*atoms, a1, setMassesToOne);
prop[a3].con_len = ip->settle.doh;
}
get_vsite_masses(&moltype,
&mtop->ffparams,
+ setMassesToOne,
vsite_m,
&n_nonlin_vsite_mol);
if (n_nonlin_vsite != nullptr)
}
else
{
- prop[a].mass = atoms->atom[a].m;
+ prop[a].mass = getMass(*atoms, a, setMassesToOne);
}
prop[a].type = atoms->atom[a].type;
prop[a].q = atoms->atom[a].q;
// Returns an (over)estimate of the energy drift for a single atom pair,
// given the kinetic properties, displacement variances and list buffer.
-static real energyDriftAtomPair(const atom_nonbonded_kinetic_prop_t *prop_i,
- const atom_nonbonded_kinetic_prop_t *prop_j,
+static real energyDriftAtomPair(bool isConstrained_i,
+ bool isConstrained_j,
real s2, real s2i_2d, real s2j_2d,
real r_buffer,
const pot_derivatives_t *der)
else
{
/* For constraints: adapt r and scaling for the Gaussian */
- if (prop_i->bConstr)
+ if (isConstrained_i)
{
real sh, sc;
rsh += sh;
sc_fac *= sc;
}
- if (prop_j->bConstr)
+ if (isConstrained_j)
{
real sh, sc;
lj.d2 = c6*ljDisp->d2 + c12*ljRep->d2;
lj.md3 = c6*ljDisp->md3 + c12*ljRep->md3;
- real pot_lj = energyDriftAtomPair(prop_i, prop_j,
+ real pot_lj = energyDriftAtomPair(prop_i->bConstr, prop_j->bConstr,
s2, s2i_2d, s2j_2d,
rlist - rlj,
&lj);
elec_qq.d2 = elec->d2 *prop_i->q*prop_j->q;
elec_qq.md3 = 0;
- real pot_q = energyDriftAtomPair(prop_i, prop_j,
+ real pot_q = energyDriftAtomPair(prop_i->bConstr, prop_j->bConstr,
s2, s2i_2d, s2j_2d,
rlist - rcoulomb,
&elec_qq);
return md3_pot + md3_sw;
}
+/* Returns the maximum reference temperature over all coupled groups */
+static real maxReferenceTemperature(const t_inputrec &ir)
+{
+ if (EI_MD(ir.eI) && ir.etc == etcNO)
+ {
+ /* This case should be handled outside calc_verlet_buffer_size */
+ gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
+ }
+
+ real maxTemperature = 0;
+ for (int i = 0; i < ir.opts.ngtc; i++)
+ {
+ if (ir.opts.tau_t[i] >= 0)
+ {
+ maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
+ }
+ }
+
+ return maxTemperature;
+}
+
+/* Returns the variance of the atomic displacement over timePeriod.
+ *
+ * Note: When not using BD with a non-mass dependendent friction coefficient,
+ * the return value still needs to be divided by the particle mass.
+ */
+static real displacementVariance(const t_inputrec &ir,
+ real temperature,
+ real timePeriod)
+{
+ real kT_fac;
+
+ if (ir.eI == eiBD)
+ {
+ /* Get the displacement distribution from the random component only.
+ * With accurate integration the systematic (force) displacement
+ * should be negligible (unless nstlist is extremely large, which
+ * you wouldn't do anyhow).
+ */
+ kT_fac = 2*BOLTZ*temperature*timePeriod;
+ if (ir.bd_fric > 0)
+ {
+ /* This is directly sigma^2 of the displacement */
+ kT_fac /= ir.bd_fric;
+ }
+ else
+ {
+ /* Per group tau_t is not implemented yet, use the maximum */
+ real tau_t = ir.opts.tau_t[0];
+ for (int i = 1; i < ir.opts.ngtc; i++)
+ {
+ tau_t = std::max(tau_t, ir.opts.tau_t[i]);
+ }
+
+ kT_fac *= tau_t;
+ /* This kT_fac needs to be divided by the mass to get sigma^2 */
+ }
+ }
+ else
+ {
+ kT_fac = BOLTZ*temperature*gmx::square(timePeriod);
+ }
+
+ return kT_fac;
+}
+
+/* Returns the largest sigma of the Gaussian displacement over all particle
+ * types. This ignores constraints, so is an overestimate.
+ */
+static real maxSigma(real kT_fac,
+ int natt,
+ const verletbuf_atomtype_t *att)
+{
+ assert(att);
+ real smallestMass = att[0].prop.mass;
+ for (int i = 1; i < natt; i++)
+ {
+ smallestMass = std::min(smallestMass, att[i].prop.mass);
+ }
+
+ return 2*std::sqrt(kT_fac/smallestMass);
+}
+
void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
const t_inputrec *ir,
int nstlist,
real nb_clust_frac_pairs_not_in_list_at_cutoff;
verletbuf_atomtype_t *att = nullptr;
- int natt = -1, i;
+ int natt = -1;
real elfac;
- real kT_fac, mass_min;
int ib0, ib1, ib;
real rb, rl;
real drift;
if (reference_temperature < 0)
{
- if (EI_MD(ir->eI) && ir->etc == etcNO)
- {
- /* This case should be handled outside calc_verlet_buffer_size */
- gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
- }
-
/* We use the maximum temperature with multiple T-coupl groups.
* We could use a per particle temperature, but since particles
* interact, this might underestimate the buffer size.
*/
- reference_temperature = 0;
- for (i = 0; i < ir->opts.ngtc; i++)
- {
- if (ir->opts.tau_t[i] >= 0)
- {
- reference_temperature = std::max(reference_temperature,
- ir->opts.ref_t[i]);
- }
- }
+ reference_temperature = maxReferenceTemperature(*ir);
}
/* Resolution of the buffer size */
/* Worst case assumption: HCP packing of particles gives largest distance */
particle_distance = std::cbrt(boxvol*std::sqrt(2)/mtop->natoms);
- get_verlet_buffer_atomtypes(mtop, &att, &natt, n_nonlin_vsite);
+ /* TODO: Obtain masses through (future) integrator functionality
+ * to avoid scattering the code with (or forgetting) checks.
+ */
+ const bool setMassesToOne = (ir->eI == eiBD && ir->bd_fric > 0);
+ get_verlet_buffer_atomtypes(mtop, setMassesToOne, &att, &natt, n_nonlin_vsite);
assert(att != nullptr && natt >= 0);
if (debug)
* For inertial dynamics (not Brownian dynamics) the mass factor
* is not included in kT_fac, it is added later.
*/
- if (ir->eI == eiBD)
- {
- /* Get the displacement distribution from the random component only.
- * With accurate integration the systematic (force) displacement
- * should be negligible (unless nstlist is extremely large, which
- * you wouldn't do anyhow).
- */
- kT_fac = 2*BOLTZ*reference_temperature*list_lifetime*ir->delta_t;
- if (ir->bd_fric > 0)
- {
- /* This is directly sigma^2 of the displacement */
- kT_fac /= ir->bd_fric;
-
- /* Set the masses to 1 as kT_fac is the full sigma^2,
- * but we divide by m in ener_drift().
- */
- for (i = 0; i < natt; i++)
- {
- att[i].prop.mass = 1;
- }
- }
- else
- {
- real tau_t;
-
- /* Per group tau_t is not implemented yet, use the maximum */
- tau_t = ir->opts.tau_t[0];
- for (i = 1; i < ir->opts.ngtc; i++)
- {
- tau_t = std::max(tau_t, ir->opts.tau_t[i]);
- }
-
- kT_fac *= tau_t;
- /* This kT_fac needs to be divided by the mass to get sigma^2 */
- }
- }
- else
- {
- kT_fac = BOLTZ*reference_temperature*gmx::square(list_lifetime*ir->delta_t);
- }
-
- mass_min = att[0].prop.mass;
- for (i = 1; i < natt; i++)
- {
- mass_min = std::min(mass_min, att[i].prop.mass);
- }
+ const real kT_fac = displacementVariance(*ir, reference_temperature,
+ list_lifetime*ir->delta_t);
if (debug)
{
fprintf(debug, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep.md1, ljRep.d2, ljRep.md3);
fprintf(debug, "Electro. -V' %9.2e V'' %9.2e\n", elec.md1, elec.d2);
fprintf(debug, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac));
- fprintf(debug, "mass_min %f\n", mass_min);
}
/* Search using bisection */
ib0 = -1;
/* The drift will be neglible at 5 times the max sigma */
- ib1 = static_cast<int>(5*2*std::sqrt(kT_fac/mass_min)/resolution) + 1;
+ ib1 = static_cast<int>(5*maxSigma(kT_fac, natt, att)/resolution) + 1;
while (ib1 - ib0 > 1)
{
ib = (ib0 + ib1)/2;
*rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution;
}
+
+/* Returns the pairlist buffer size for use as a minimum buffer size
+ *
+ * Note that this is a rather crude estimate. It is ok for a buffer
+ * set for good energy conservation or RF electrostatics. But it is
+ * too small with PME and the buffer set with the default tolerance.
+ */
+static real minCellSizeFromPairlistBuffer(const t_inputrec &ir)
+{
+ return ir.rlist - std::max(ir.rvdw, ir.rcoulomb);
+}
+
+real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ real chanceRequested)
+{
+ if (!EI_DYNAMICS(ir.eI) || (EI_MD(ir.eI) && ir.etc == etcNO))
+ {
+ return minCellSizeFromPairlistBuffer(ir);
+ }
+
+ /* We use the maximum temperature with multiple T-coupl groups.
+ * We could use a per particle temperature, but since particles
+ * interact, this might underestimate the displacements.
+ */
+ const real temperature = maxReferenceTemperature(ir);
+
+ const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0);
+
+ verletbuf_atomtype_t *att = nullptr;
+ int natt = -1;
+ get_verlet_buffer_atomtypes(&mtop, setMassesToOne, &att, &natt, nullptr);
+
+ const real kT_fac = displacementVariance(ir, temperature,
+ ir.nstlist*ir.delta_t);
+
+ /* Resolution of the cell size */
+ real resolution = 0.001;
+
+ /* Search using bisection, avoid 0 and start at 1 */
+ int ib0 = 0;
+ /* The chance will be neglible at 10 times the max sigma */
+ int ib1 = int(10*maxSigma(kT_fac, natt, att)/resolution) + 1;
+ real cellSize = 0;
+ while (ib1 - ib0 > 1)
+ {
+ int ib = (ib0 + ib1)/2;
+ cellSize = ib*resolution;
+
+ /* We assumes atom are distributed uniformly over the cell width.
+ * Once an atom has moved by more than the cellSize (as passed
+ * as the buffer argument to energyDriftAtomPair() below),
+ * the chance of crossing the boundary of the neighbor cell
+ * thus increases as 1/cellSize with the additional displacement
+ * on to of cellSize. We thus create a linear interaction with
+ * derivative = -1/cellSize. Using this in the energyDriftAtomPair
+ * function will return the chance of crossing the next boundary.
+ */
+ const pot_derivatives_t boundaryInteraction = { 1/cellSize, 0, 0 };
+
+ real chance = 0;
+ for (int i = 0; i < natt; i++)
+ {
+ const atom_nonbonded_kinetic_prop_t &propAtom = att[i].prop;
+ real s2_2d;
+ real s2_3d;
+ get_atom_sigma2(kT_fac, &propAtom, &s2_2d, &s2_3d);
+
+ real chancePerAtom = energyDriftAtomPair(propAtom.bConstr, false,
+ s2_2d + s2_3d, s2_2d, 0,
+ cellSize,
+ &boundaryInteraction);
+
+ if (propAtom.bConstr)
+ {
+ /* energyDriftAtomPair() uses an unlimited Gaussian displacement
+ * distribution for constrained atoms, whereas they can
+ * actually not move more than the COM of the two constrained
+ * atoms plus twice the distance from the COM.
+ * Use this maximum, limited displacement when this results in
+ * a smaller chance (note that this is still an overestimate).
+ */
+ real massFraction = propAtom.con_mass/(propAtom.mass + propAtom.con_mass);
+ real comDistance = propAtom.con_len*massFraction;
+
+ real chanceWithMaxDistance =
+ energyDriftAtomPair(false, false,
+ s2_3d, 0, 0,
+ cellSize - 2*comDistance,
+ &boundaryInteraction);
+ chancePerAtom = std::min(chancePerAtom, chanceWithMaxDistance);
+ }
+
+ /* Take into account the line density of the boundary */
+ chancePerAtom /= cellSize;
+
+ chance += att[i].n*chancePerAtom;
+ }
+
+ /* Note: chance is for every nstlist steps */
+ if (chance > chanceRequested*ir.nstlist)
+ {
+ ib0 = ib;
+ }
+ else
+ {
+ ib1 = ib;
+ }
+ }
+
+ sfree(att);
+
+ return cellSize;
+}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2013,2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
int *n_nonlin_vsite,
real *rlist);
+/* Determines the mininum cell size based on atom displacement
+ *
+ * The value returned is the minimum size for which the chance that
+ * an atom crosses to non nearest-neighbor cells is <= chanceRequested
+ * within ir.nstlist steps.
+ * Without T-coupling, SD or BD, we can not estimate atom displacements
+ * and fall back to the, crude, estimate of using the pairlist buffer size.
+ *
+ * Note: Like the Verlet buffer estimate, this estimate is based on
+ * non-interacting atoms and constrained atom-pairs. Therefore for
+ * any system that is not an ideal gas, this will be an overestimate.
+ *
+ * Note: This size increases (very slowly) with system size.
+ */
+real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ real chanceRequested);
+
/* Struct for unique atom type for calculating the energy drift.
* The atom displacement depends on mass and constraints.
* The energy jump for given distance depend on LJ type and q.
}
else if (i == F_COM_PULL)
{
- md->bEner[i] = (ir->bPull && pull_have_potential(ir->pull_work));
+ md->bEner[i] = ((ir->bPull && pull_have_potential(ir->pull_work)) ||
+ ir->bRot);
}
else if (i == F_ECONSERVED)
{
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/nb_verlet.h"
#include "gromacs/mdlib/nbnxn_consts.h"
-#include "gromacs/mdlib/nbnxn_gpu_common_utils.h"
#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
const nbnxn_pairlist_t *h_plist,
int iloc)
{
- if (canSkipWork(nb, iloc))
- {
- return;
- }
-
char sbuf[STRLEN];
- bool bDoTime = nb->bDoTime;
+ bool bDoTime = (nb->bDoTime && h_plist->nsci > 0);
cudaStream_t stream = nb->stream[iloc];
cu_plist_t *d_plist = nb->plist[iloc];
#include "gromacs/mdlib/nb_verlet.h"
#include "gromacs/mdlib/nbnxn_consts.h"
#include "gromacs/mdlib/nbnxn_gpu.h"
-#include "gromacs/mdlib/nbnxn_gpu_common_utils.h"
#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
#include "gromacs/mdlib/nbnxn_gpu_jit_support.h"
#include "gromacs/mdtypes/interaction_const.h"
const nbnxn_pairlist_t *h_plist,
int iloc)
{
- if (canSkipWork(nb, iloc))
- {
- return;
- }
-
char sbuf[STRLEN];
- bool bDoTime = nb->bDoTime;
+ // Timing accumulation should happen only if there was work to do
+ // because getLastRangeTime() gets skipped with empty lists later
+ // which leads to the counter not being reset.
+ bool bDoTime = (nb->bDoTime && h_plist->nsci > 0);
cl_command_queue stream = nb->stream[iloc];
cl_plist_t *d_plist = nb->plist[iloc];
static void reset_nblist(t_nblist *nl)
{
+ GMX_RELEASE_ASSERT(nl, "Should only reset valid nblists");
+
nl->nri = -1;
nl->nrj = 0;
if (nl->jindex)
{
int n, i;
- if (fr->bQMMM)
+ if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
{
/* only reset the short-range nblist */
reset_nblist(fr->QMMMlist);
#include <vector>
#include "gromacs/commandline/filenm.h"
+#include "gromacs/domdec/collect.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/ewald/pme.h"
&state->s, state_global, observablesHistory,
state->f);
- if (confout != nullptr && MASTER(cr))
+ if (confout != nullptr)
{
- GMX_RELEASE_ASSERT(bX, "The code below assumes that (with domain decomposition), x is collected to state_global in the call above.");
- /* With domain decomposition the call above collected the state->s.x
- * into state_global->x. Without DD we copy the local state pointer.
- */
- if (!DOMAINDECOMP(cr))
+ if (DOMAINDECOMP(cr))
{
+ /* If bX=true, x was collected to state_global in the call above */
+ if (!bX)
+ {
+ gmx::ArrayRef<gmx::RVec> globalXRef = MASTER(cr) ? gmx::makeArrayRef(state_global->x) : gmx::EmptyArrayRef();
+ dd_collect_vec(cr->dd, &state->s, state->s.x, globalXRef);
+ }
+ }
+ else
+ {
+ /* Copy the local state pointer */
state_global = &state->s;
}
- if (ir->ePBC != epbcNONE && !ir->bPeriodicMols && DOMAINDECOMP(cr))
+ if (MASTER(cr))
{
- /* Make molecules whole only for confout writing */
- do_pbc_mtop(fplog, ir->ePBC, state->s.box, top_global,
- as_rvec_array(state_global->x.data()));
- }
+ if (ir->ePBC != epbcNONE && !ir->bPeriodicMols && DOMAINDECOMP(cr))
+ {
+ /* Make molecules whole only for confout writing */
+ do_pbc_mtop(fplog, ir->ePBC, state->s.box, top_global,
+ as_rvec_array(state_global->x.data()));
+ }
- write_sto_conf_mtop(confout,
- *top_global->name, top_global,
- as_rvec_array(state_global->x.data()), nullptr, ir->ePBC, state->s.box);
+ write_sto_conf_mtop(confout,
+ *top_global->name, top_global,
+ as_rvec_array(state_global->x.data()), nullptr, ir->ePBC, state->s.box);
+ }
}
}
step = 0;
- // Ensure the extra per-atom state array gets allocated
- state_global->flags |= (1<<estCGP);
+ if (MASTER(cr))
+ {
+ // In CG, the state is extended with a search direction
+ state_global->flags |= (1<<estCGP);
+
+ // Ensure the extra per-atom state array gets allocated
+ state_change_natoms(state_global, state_global->natoms);
+
+ // Initialize the search direction to zero
+ for (RVec &cg_p : state_global->cg_p)
+ {
+ cg_p = { 0, 0, 0 };
+ }
+ }
/* Create 4 states on the stack and extract pointers that we will swap */
em_state_t s0 {}, s1 {}, s2 {}, s3 {};
gmx_sumd(1, &minstep, cr);
}
- minstep = GMX_REAL_EPS/sqrt(minstep/(3*state_global->natoms));
+ minstep = GMX_REAL_EPS/sqrt(minstep/(3*top_global->natoms));
if (stepsize < minstep)
{
}
/* Send energies and positions to the IMD client if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
+ if (MASTER(cr) && do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle))
{
IMD_send_positions(inputrec->imd);
}
* However, we should only do it if we did NOT already write this step
* above (which we did if do_x or do_f was true).
*/
+ /* Note that with 0 < nstfout != nstxout we can end up with two frames
+ * in the trajectory with the same step number.
+ */
do_x = !do_per_step(step, inputrec->nstxout);
do_f = (inputrec->nstfout > 0 && !do_per_step(step, inputrec->nstfout));
const rvec x[],
rvec *xp);
+/*! \brief Margin for checking PBC distances compared to half the box size in pullCheckPbcWithinGroups() */
+static constexpr real c_pullGroupPbcMargin = 0.9;
+
+/*! \brief Checks whether all groups that use a reference atom are within PBC restrictions
+ *
+ * Groups that use a reference atom for determining PBC should have all their
+ * atoms within half the box size from the PBC atom. The box size is used
+ * per dimension for rectangular boxes, but can be a combination of
+ * dimensions for triclinic boxes, depending on which dimensions are
+ * involved in the pull coordinates a group is involved in.
+ *
+ * Should be called without MPI parallelization and after pull_calc_coms()
+ * has been called at least once.
+ *
+ * \param[in] pull The pull data structure
+ * \param[in] x The coordinates
+ * \param[in] pbc Information struct about periodicity
+ * \returns -1 when all groups obey PBC or the first group index that fails PBC
+ */
+int pullCheckPbcWithinGroups(const pull_t &pull,
+ const rvec *x,
+ const t_pbc &pbc);
/*! \brief Returns if we have pull coordinates with potential pulling.
*
make_cyl_refgrps(cr, pull, md, pbc, t, x);
}
}
+
+using BoolVec = gmx::BasicVector<bool>;
+
+/* Returns whether the pull group obeys the PBC restrictions */
+static bool pullGroupObeysPbcRestrictions(const pull_group_work_t &group,
+ const BoolVec &dimUsed,
+ const rvec *x,
+ const t_pbc &pbc,
+ const rvec &x_pbc)
+{
+ /* Determine which dimensions are relevant for PBC */
+ BoolVec dimUsesPbc = { false, false, false };
+ bool pbcIsRectangular = true;
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsed[d])
+ {
+ dimUsesPbc[d] = true;
+ /* All non-zero dimensions of vector v are involved in PBC */
+ for (int d2 = d + 1; d2 < pbc.ndim_ePBC; d2++)
+ {
+ assert(d2 < DIM);
+ if (pbc.box[d2][d] != 0)
+ {
+ dimUsesPbc[d2] = true;
+ pbcIsRectangular = false;
+ }
+ }
+ }
+ }
+
+ rvec marginPerDim = {};
+ real marginDistance2 = 0;
+ if (pbcIsRectangular)
+ {
+ /* Use margins for dimensions independently */
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ marginPerDim[d] = c_pullGroupPbcMargin*pbc.hbox_diag[d];
+ }
+ }
+ else
+ {
+ /* Check the total distance along the relevant dimensions */
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsesPbc[d])
+ {
+ marginDistance2 += c_pullGroupPbcMargin*gmx::square(0.5)*norm2(pbc.box[d]);
+ }
+ }
+ }
+
+ auto localAtomIndices = group.atomSet.localIndex();
+ for (gmx::index indexInSet = 0; indexInSet < localAtomIndices.size(); indexInSet++)
+ {
+ rvec dx;
+ pbc_dx(&pbc, x[indexInSet], x_pbc, dx);
+
+ bool atomIsTooFar = false;
+ if (pbcIsRectangular)
+ {
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsesPbc[d] && (dx[d] < -marginPerDim[d] ||
+ dx[d] > marginPerDim[d]))
+ {
+ atomIsTooFar = true;
+ }
+ }
+ }
+ else
+ {
+ real pbcDistance2 = 0;
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsesPbc[d])
+ {
+ pbcDistance2 += gmx::square(dx[d]);
+ }
+ }
+ atomIsTooFar = (pbcDistance2 > marginDistance2);
+ }
+ if (atomIsTooFar)
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+int pullCheckPbcWithinGroups(const pull_t &pull,
+ const rvec *x,
+ const t_pbc &pbc)
+{
+ if (pbc.ePBC == epbcNONE)
+ {
+ return -1;
+ }
+
+ /* Determine what dimensions are used for each group by pull coordinates */
+ std::vector<BoolVec> dimUsed(pull.group.size(), { false, false, false });
+ for (size_t c = 0; c < pull.coord.size(); c++)
+ {
+ const t_pull_coord &coordParams = pull.coord[c].params;
+ for (int groupIndex = 0; groupIndex < coordParams.ngroup; groupIndex++)
+ {
+ for (int d = 0; d < DIM; d++)
+ {
+ if (coordParams.dim[d] &&
+ !(coordParams.eGeom == epullgCYL && groupIndex == 0))
+ {
+ dimUsed[coordParams.group[groupIndex]][d] = true;
+ }
+ }
+ }
+ }
+
+ /* Check PBC for every group that uses a PBC reference atom treatment */
+ for (size_t g = 0; g < pull.group.size(); g++)
+ {
+ const pull_group_work_t &group = pull.group[g];
+ if (group.epgrppbc == epgrppbcREFAT &&
+ !pullGroupObeysPbcRestrictions(group, dimUsed[g], x, pbc, pull.comm.rbuf[g]))
+ {
+ return g;
+ }
+ }
+
+ return -1;
+}
}
static inline bool gmx_simdcall
-anyTrue(SimdDIBool a) { return _mm_movemask_epi8(_mm_shuffle_epi32(a.simdInternal_, _MM_SHUFFLE(1, 0, 1, 0))) != 0; }
+anyTrue(SimdDIBool a) { return _mm_movemask_epi8(a.simdInternal_) != 0; }
static inline SimdDInt32 gmx_simdcall
selectByMask(SimdDInt32 a, SimdDIBool mask)
#include <cmath>
+#include <array>
+
#include "gromacs/math/utilities.h"
#include "gromacs/simd/simd.h"
#include "gromacs/utility/basedefinitions.h"
TEST_F(SimdFloatingpointTest, anyTrueB)
{
- SimdBool eq;
+ alignas(GMX_SIMD_ALIGNMENT) std::array<real, GMX_SIMD_REAL_WIDTH> mem {};
- /* this test is a bit tricky since we don't know the simd width.
- * We cannot check for truth values for "any" element beyond the first,
- * since that part of the data will not be used if simd width is 1.
- */
- eq = rSimd_c4c6c8 == setSimdRealFrom3R(c4, 0, 0);
- EXPECT_TRUE(anyTrue(eq));
+ // Test the false case
+ EXPECT_FALSE(anyTrue(setZero() < load<SimdReal>(mem.data())));
- eq = rSimd_c0c1c2 == rSimd_c3c4c5;
- EXPECT_FALSE(anyTrue(eq));
+ // Test each bit (these should all be true)
+ for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
+ {
+ mem.fill(0.0);
+ mem[i] = 1.0;
+ EXPECT_TRUE(anyTrue(setZero() < load<SimdReal>(mem.data()))) << "Not detecting true in element " << i;
+ }
}
TEST_F(SimdFloatingpointTest, blend)
*/
#include "gmxpre.h"
+#include <array>
+
#include "gromacs/simd/simd.h"
#include "gromacs/utility/basedefinitions.h"
TEST_F(SimdIntegerTest, anyTrue)
{
- SimdIBool eq;
+ alignas(GMX_SIMD_ALIGNMENT) std::array<std::int32_t, GMX_SIMD_REAL_WIDTH> mem {};
- /* See comment in floatingpoint.cpp. We should only check the first element here,
- * since the SIMD width could be 1 as a special case.
- */
- eq = (iSimd_5_7_9 == setSimdIntFrom3I(5, 0, 0));
- EXPECT_TRUE(anyTrue(eq));
+ // Test the false case
+ EXPECT_FALSE(anyTrue(setZero() < load<SimdInt32>(mem.data())));
- eq = (iSimd_1_2_3 == iSimd_4_5_6);
- EXPECT_FALSE(anyTrue(eq));
+ // Test each bit (these should all be true)
+ for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
+ {
+ mem.fill(0);
+ mem[i] = 1;
+ EXPECT_TRUE(anyTrue(setZero() < load<SimdInt32>(mem.data()))) << "Not detecting true in element " << i;
+ }
}
TEST_F(SimdIntegerTest, blend)
writer->writeLine("TNG support: disabled");
#endif
#if GMX_USE_HWLOC
- writer->writeLine(formatString("Hwloc support: hwloc-%d.%d.%d",
- HWLOC_API_VERSION>>16,
- (HWLOC_API_VERSION>>8) & 0xFF,
- HWLOC_API_VERSION & 0xFF));
+ writer->writeLine(formatString("Hwloc support: hwloc-%s", HWLOC_VERSION));
#else
writer->writeLine("Hwloc support: disabled");
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff