Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Run NpT sanity checks at integration start #4165

Merged
merged 4 commits into from
Mar 29, 2021
Merged

Run NpT sanity checks at integration start #4165

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d3afabc
core: Remove unreachable code
jngrad Mar 19, 2021
1616de2
core: Run NpT sanity checks at integration start
jngrad Mar 19, 2021
c7a88ff
tests: Check NpT sanity checks
jngrad Mar 19, 2021
900e065
Merge branch 'python' into fix-4139
kodiakhq[bot] Mar 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 28 additions & 29 deletions src/core/npt.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,28 @@ void mpi_bcast_nptiso_geom_barostat() {
mpi_call_all(mpi_bcast_nptiso_geom_barostat_worker);
}

void integrator_npt_coulomb_dipole_sanity_checks(nptiso_struct const &params) {
#ifdef ELECTROSTATICS
if (params.dimension < 3 && !params.cubic_box && coulomb.prefactor > 0) {
throw std::runtime_error("If electrostatics is being used you must "
"use the cubic box NpT.");
}
#endif

#ifdef DIPOLES
if (params.dimension < 3 && !params.cubic_box && dipole.prefactor > 0) {
throw std::runtime_error("If magnetostatics is being used you must "
"use the cubic box NpT.");
}
#endif

#if defined(ELECTROSTATICS) && defined(CUDA)
if (coulomb.method == COULOMB_P3M_GPU) {
throw std::runtime_error("NpT virial cannot be calculated on P3M GPU");
}
#endif
}

void nptiso_init(double ext_pressure, double piston, bool xdir_rescale,
bool ydir_rescale, bool zdir_rescale, bool cubic_box) {

Expand Down Expand Up @@ -114,27 +136,14 @@ void nptiso_init(double ext_pressure, double piston, bool xdir_rescale,
new_nptiso.non_const_dim = 2;
}

/* Sanity Checks */
#ifdef ELECTROSTATICS
if (new_nptiso.dimension < 3 && !cubic_box && coulomb.prefactor > 0) {
throw std::runtime_error("If electrostatics is being used you must "
"use the cubic box npt.");
}
#endif

#ifdef DIPOLES
if (new_nptiso.dimension < 3 && !cubic_box && dipole.prefactor > 0) {
throw std::runtime_error("If magnetostatics is being used you must "
"use the cubic box npt.");
}
#endif

if (new_nptiso.dimension == 0 || new_nptiso.non_const_dim == -1) {
throw std::runtime_error(
"You must enable at least one of the x y z components "
"as fluctuating dimension(s) for box length motion!");
}

integrator_npt_coulomb_dipole_sanity_checks(new_nptiso);

nptiso = new_nptiso;

mpi_bcast_nptiso_geom_barostat();
Expand All @@ -145,14 +154,7 @@ void npt_ensemble_init(const BoxGeometry &box) {
if (integ_switch == INTEG_METHOD_NPT_ISO) {
/* prepare NpT-integration */
nptiso.inv_piston = 1 / nptiso.piston;
if (nptiso.dimension == 0) {
throw std::runtime_error(
"INTERNAL ERROR: npt integrator was called but "
"dimension not yet set. This should not happen.");
}

nptiso.volume = pow(box.length()[nptiso.non_const_dim], nptiso.dimension);

if (recalc_forces) {
nptiso.p_inst = 0.0;
nptiso.p_vir = Utils::Vector3d{};
Expand All @@ -163,14 +165,11 @@ void npt_ensemble_init(const BoxGeometry &box) {

void integrator_npt_sanity_checks() {
if (integ_switch == INTEG_METHOD_NPT_ISO) {
if (nptiso.piston <= 0.0) {
runtimeErrorMsg() << "npt on, but piston mass not set";
}
#if defined(ELECTROSTATICS) && defined(CUDA)
if (coulomb.method == COULOMB_P3M_GPU) {
runtimeErrorMsg() << "npt on, but virial cannot be calculated on P3M GPU";
try {
integrator_npt_coulomb_dipole_sanity_checks(nptiso);
} catch (std::runtime_error const &err) {
runtimeErrorMsg() << err.what();
}
#endif
}
}

Expand Down
6 changes: 1 addition & 5 deletions src/core/thermostat.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,11 +110,7 @@ void thermo_init() {
#endif
#ifdef NPT
if (thermo_switch & THERMO_NPT_ISO) {
if (nptiso.piston == 0.0) {
thermo_switch = (thermo_switch ^ THERMO_NPT_ISO);
} else {
npt_iso.recalc_prefactors(nptiso.piston, time_step);
}
npt_iso.recalc_prefactors(nptiso.piston, time_step);
}
#endif
if (thermo_switch & THERMO_BROWNIAN)
Expand Down
47 changes: 38 additions & 9 deletions testsuite/python/integrator_npt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#
import unittest as ut
import unittest_decorators as utx
import tests_common
import numpy as np
import espressomd
import espressomd.integrate
Expand Down Expand Up @@ -115,31 +116,59 @@ def test_integrator_recovery(self):
np.copy(system.part[:].pos),
positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]]))

def run_with_p3m(self, p3m):
def run_with_p3m(self, p3m, **npt_kwargs):
system = self.system
np.random.seed(42)
# set up particles
system.box_l = [6] * 3
system.part.add(pos=np.random.uniform(0, system.box_l[0], (11, 3)),
q=np.sign(np.arange(-5, 6)))
system.part.add(pos=np.random.uniform(0, system.box_l[0], (11, 3)))
if espressomd.has_features("P3M"):
system.part[:].q = np.sign(np.arange(-5, 6))
if espressomd.has_features("DP3M"):
system.part[:].dip = tests_common.random_dipoles(11)
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=1, sigma=1, cutoff=2**(1 / 6), shift=0.25)
system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42)
system.integrator.set_isotropic_npt(ext_pressure=0.001, piston=0.001)
system.integrator.set_steepest_descent(
f_max=10, gamma=0.1, max_displacement=0.01)
system.integrator.run(100)
system.integrator.set_vv()
# combine NpT with a P3M algorithm
system.actors.add(p3m)
system.integrator.run(100)
system.integrator.run(20)
system.integrator.set_isotropic_npt(ext_pressure=0.001, piston=0.001,
**npt_kwargs)
system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42)
system.integrator.run(20)

@utx.skipIfMissingFeatures(["DP3M"])
def test_dp3m_exception(self):
# NpT is compatible with DP3M CPU (only cubic box)
import espressomd.magnetostatics
dp3m = espressomd.magnetostatics.DipolarP3M(
prefactor=1.0, accuracy=1e-2, mesh=3 * [36], cao=7, r_cut=1.0,
alpha=2.995, tune=False)
with self.assertRaisesRegex(RuntimeError, 'If magnetostatics is being used you must use the cubic box NpT'):
self.run_with_p3m(dp3m, cubic_box=False, direction=(0, 1, 1))
self.tearDown()
try:
self.run_with_p3m(dp3m)
except Exception as err:
self.fail(f'integrator raised ValueError("{err}")')

@utx.skipIfMissingFeatures(["P3M"])
def test_p3m_exception(self):
# NpT is compatible with P3M CPU
# NpT is compatible with P3M CPU (only cubic box)
import espressomd.electrostatics
p3m = espressomd.electrostatics.P3M(
prefactor=1.0, accuracy=1e-2, mesh=3 * [8], cao=3, r_cut=0.36,
alpha=5.35, tune=False)
with self.assertRaisesRegex(RuntimeError, 'If electrostatics is being used you must use the cubic box NpT'):
self.run_with_p3m(p3m, cubic_box=False, direction=(0, 1, 1))
self.tearDown()
try:
self.run_with_p3m(p3m)
except Exception as err:
self.fail(f'run_with_p3m() raised ValueError("{err}")')
self.fail(f'integrator raised ValueError("{err}")')

@utx.skipIfMissingGPU()
@utx.skipIfMissingFeatures(["P3M"])
Expand All @@ -149,7 +178,7 @@ def test_p3mgpu_exception(self):
p3m = espressomd.electrostatics.P3MGPU(
prefactor=1.0, accuracy=1e-2, mesh=3 * [24], cao=2, r_cut=0.24,
alpha=8.26, tune=False)
with self.assertRaises(Exception):
with self.assertRaisesRegex(RuntimeError, 'NpT virial cannot be calculated on P3M GPU'):
self.run_with_p3m(p3m)


Expand Down