Skip to content
Permalink
Browse files

Merge #3073

3073: Use steepest descent without thermostat in tutorial 02 r=fweik a=jngrad

Fixes #3030

Description of changes:
- cleanup and document the steepest descent algorithm
- reflect recent changes to the Python interface in the tutorials: `mindist()` renamed to `min_dist()`, periodicity takes a list of booleans
- in tutorial 02 charged system:
    - set Langevin thermostat after energy minimization
    - replace adaptive force-cap minimization strategy with steepest descent

Co-authored-by: Jean-Noël Grad <jgrad@icp.uni-stuttgart.de>
  • Loading branch information...
bors and jngrad committed Aug 12, 2019
2 parents da6b58d + 980210f commit 60f4f7a685c4461bf2403624f31dbf1c7552a8ed
@@ -96,7 +96,7 @@ For example, ::
7.0
>>> system.analysis.dist_to(pos=[0, 0, 0])
1.4142135623730951
>>> system.analysis.mindist()
>>> system.analysis.min_dist()
1.0


@@ -13,7 +13,7 @@
"source": [
"## 1 Introduction\n",
"\n",
"This tutorial introduces some of the basic features of **ESPResSo** for charged systems by constructing a simulation script for a simple salt crystal. In the subsequent task, we use a more realistic force-field for a NaCl crystal. Finally, we introduce constraints and 2D-Electrostatics to simulate a molten salt in a parallel plate capacitor. We assume that the reader is familiar with the basic concepts of Python and MD simulations. Compile espresso with the following features in your <tt>myconfig.hpp</tt> to be set throughout the whole tutorial:\n",
"This tutorial introduces some of the basic features of **ESPResSo** for charged systems by constructing a simulation script for a simple salt crystal. In the subsequent task, we use a more realistic force-field for a NaCl crystal. Finally, we introduce constraints and 2D-Electrostatics to simulate a molten salt in a parallel plate capacitor. We assume that the reader is familiar with the basic concepts of Python and MD simulations. Compile **ESPResSo** with the following features in your <tt>myconfig.hpp</tt> to be set throughout the whole tutorial:\n",
"\n",
"```c++\n",
"#define EXTERNAL_FORCES\n",
@@ -90,11 +90,10 @@
"equally sized, monovalent salt.\n",
"\n",
"The simulation engine itself is modified by changing the\n",
"espressomd.System() properties. We create an instance <tt>system</tt> and\n",
"<tt>espressomd.System()</tt> properties. We create an instance <tt>system</tt> and\n",
"set the box length, periodicity and time step. The skin depth <tt>skin</tt> \n",
"is a parameter for the link--cell system which tunes its\n",
"performance, but shall not be discussed here. Further, we activate the Langevin thermostat\n",
"for our NVT ensemble with temperature <tt>temp</tt> and friction coefficient <tt>gamma</tt>. "
"performance, but shall not be discussed here."
]
},
{
@@ -107,10 +106,9 @@
"box_l = (n_part / density)**(1. / 3.)\n",
"system = System(box_l=[box_l, box_l, box_l])\n",
"system.seed = 42\n",
"system.periodicity = [1, 1, 1]\n",
"system.periodicity = [True, True, True]\n",
"system.time_step = time_step\n",
"system.cell_system.skin = 0.3\n",
"system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)"
"system.cell_system.skin = 0.3"
]
},
{
@@ -186,11 +184,11 @@
"With randomly positioned particles, we most likely have huge overlap and the strong repulsion will\n",
"cause the simulation to crash. The next step in our script therefore is a suitable LJ equilibration.\n",
"This is known to be a tricky part of a simulation and several approaches exist to reduce the particle overlap.\n",
"Here, we use a highly damped system (large gamma in the thermostat) and cap the forces of the LJ interaction.\n",
"We use <tt>system.analysis.mindist</tt> to get the minimal distance between all particles pairs. This value\n",
"is used to progressively increase the force capping. This results in a slow increase of the force capping at\n",
"strong overlap. At the end, we reset our thermostat to the target values and deactivate the force cap by setting \n",
"it to zero."
"Here, we use the steepest descent algorithm and cap the maximal particle displacement per integration step\n",
"to 1% of sigma.\n",
"We use <tt>system.analysis.min_dist()</tt> to get the minimal distance between all particles pairs. This value\n",
"is used to stop the minimization when particles are far away enough from each other. At the end, we activate the\n",
"Langevin thermostat for our NVT ensemble with temperature <tt>temp</tt> and friction coefficient <tt>gamma</tt>."
]
},
{
@@ -202,21 +200,16 @@
"# Lennard-Jones Equilibration\n",
"max_sigma = max(lj_sigmas.values())\n",
"min_dist = 0.0\n",
"cap = 10.0\n",
"# Warmup Helper: Cold, highly damped system\n",
"system.thermostat.set_langevin(kT=temp * 0.1, gamma=gamma * 50.0, seed=42)\n",
"system.minimize_energy.init(f_max=0, gamma=10.0, max_steps=10,\n",
" max_displacement=max_sigma * 0.01)\n",
"\n",
"while min_dist < max_sigma:\n",
" # Warmup Helper: Cap max. force, increase slowly for overlapping particles\n",
" min_dist = system.analysis.min_dist([types[\"Anion\"], types[\"Cation\"]],\n",
" [types[\"Anion\"], types[\"Cation\"]])\n",
" cap += min_dist\n",
" system.force_cap = cap\n",
" system.integrator.run(10)\n",
" system.minimize_energy.minimize()\n",
"\n",
"# Don't forget to reset thermostat, timestep and force cap\n",
"system.thermostat.set_langevin(kT=temp, gamma=gamma)\n",
"system.force_cap = 0"
"# Set thermostat\n",
"system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)"
]
},
{
@@ -159,10 +159,9 @@
"system.seed = 42\n",
"box_volume = numpy.prod([box_l, box_l, box_z])\n",
"\n",
"system.periodicity = [1, 1, 1]\n",
"system.periodicity = [True, True, True]\n",
"system.time_step = time_step\n",
"system.cell_system.skin = 0.3\n",
"system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)"
"system.cell_system.skin = 0.3"
]
},
{
@@ -258,9 +257,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Next is the Lennard-Jones Equilibration: Here we use an alternative way to get rid of the overlap: **ESPResSo** features a routine for energy\n",
"minimization with similar features as in the manual implementation used before. Basically it \n",
"caps the forces and limits the displacement during integration."
"Next is the Lennard-Jones equilibration, followed by the thermostat:"
]
},
{
@@ -270,12 +267,15 @@
"outputs": [],
"source": [
"energy = system.analysis.energy()\n",
"print(\"Before Minimization: E_total=\", energy['total'])\n",
"print(\"Before Minimization: E_total = {:.3e}\".format(energy['total']))\n",
"system.minimize_energy.init(f_max=10, gamma=10, max_steps=2000,\n",
" max_displacement=0.1)\n",
" max_displacement=0.01)\n",
"system.minimize_energy.minimize()\n",
"energy = system.analysis.energy()\n",
"print(\"After Minimization: E_total=\", energy['total'])"
"print(\"After Minimization: E_total = {:.3e}\".format(energy['total']))\n",
"\n",
"# Set thermostat\n",
"system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)"
]
},
{
@@ -424,7 +424,7 @@
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"plt.ion()\n"
"plt.ion()"
]
},
{
@@ -54,10 +54,9 @@
box_l = (n_part / density)**(1. / 3.)
system = espressomd.System(box_l=[box_l] * 3)
system.seed = system.cell_system.get_state()['n_nodes'] * [1234]
system.periodicity = [1, 1, 1]
system.periodicity = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)

# Place particles
for i in range(int(n_ionpairs)):
@@ -97,21 +96,15 @@ def combination_rule_sigma(rule, sig1, sig2):
print("\n--->Lennard-Jones Equilibration")
max_sigma = max(lj_sigmas.values())
min_dist = 0.0
cap = 10.0
# Warmup Helper: Cold, highly damped system
system.thermostat.set_langevin(kT=temp * 0.1, gamma=gamma * 50.0)
system.minimize_energy.init(f_max=0, gamma=10, max_steps=10,
max_displacement=max_sigma * 0.01)

while min_dist < max_sigma:
# Warmup Helper: Cap max. force, increase slowly for overlapping particles
min_dist = system.analysis.min_dist([types["Anion"], types["Cation"]], [
types["Anion"], types["Cation"]])
cap += min_dist
system.force_cap = cap
system.integrator.run(10)

# Don't forget to reset thermostat, timestep and force cap
system.thermostat.set_langevin(kT=temp, gamma=gamma)
system.force_cap = 0
system.minimize_energy.minimize()
min_dist = system.analysis.min_dist()

# Set thermostat
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)

print("\n--->Tuning Electrostatics")
p3m = electrostatics.P3M(prefactor=l_bjerrum, accuracy=1e-3)
@@ -60,7 +60,7 @@
# Setup System
box_l = (n_ionpairs * sum(masses.values()) / density)**(1. / 3.)
system.box_l = [box_l, box_l, box_l]
system.periodicity = [1, 1, 1]
system.periodicity = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)
@@ -64,7 +64,7 @@
box_volume = box_l * box_l * box_z
elc_gap = box_z * 0.15
system.box_l = [box_l, box_l, box_z + elc_gap]
system.periodicity = [1, 1, 1]
system.periodicity = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)
@@ -65,7 +65,7 @@
box_volume = box_l * box_l * box_z
elc_gap = box_z * 0.15
system.box_l = [box_l, box_l, box_z + elc_gap]
system.periodicity = [1, 1, 1]
system.periodicity = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)
@@ -19,7 +19,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import espressomd
from espressomd import assert_features, electrostatics, electrostatic_extensions
from espressomd import assert_features, electrostatics
from espressomd import visualization_opengl
import numpy
from threading import Thread
@@ -97,7 +97,7 @@ def decreaseTemp():
# Setup System
box_l = (n_ionpairs * sum(masses.values()) / density)**(1. / 3.)
system.box_l = [box_l, box_l, box_l]
system.periodicity = [1, 1, 1]
system.periodicity = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)
@@ -143,7 +143,7 @@
"metadata": {},
"outputs": [],
"source": [
"system.periodicity = [1, 1, 1]"
"system.periodicity = [True, True, True]"
]
},
{
@@ -28,6 +28,7 @@

#include <utils/math/sqr.hpp>

#include <boost/algorithm/clamp.hpp>
#include <boost/mpi/collectives/all_reduce.hpp>
#include <boost/mpi/operations.hpp>

@@ -40,32 +41,17 @@
#define MINIMIZE_ENERGY_TRACE(A)
#endif

struct MinimizeEnergyParameters {
double f_max;
double gamma;
int max_steps;
double max_displacement;
};

/** Currently active steepest descent instance */
static MinimizeEnergyParameters *params = nullptr;

/* Signum of val */
template <typename T> int sgn(T val) { return (T(0) < val) - (val < T(0)); }

bool steepest_descent_step(const ParticleRange &particles) {
// Maximal force encountered on node
double f_max = -std::numeric_limits<double>::max();
// and globally

// Positional increments
double dp, dp2, dp2_max = -std::numeric_limits<double>::max();
auto f_max = -std::numeric_limits<double>::max();

// Iteration over all local particles

for (auto &p : particles) {
auto f = 0.0;

dp2 = 0.0;
// For all Cartesian coordinates
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
@@ -80,12 +66,10 @@ bool steepest_descent_step(const ParticleRange &particles) {
// Square of force on particle
f += Utils::sqr(p.f.f[j]);

// Positional increment
dp = params->gamma * p.f.f[j];
if (fabs(dp) > params->max_displacement)
// Crop to maximum allowed by user
dp = sgn<double>(dp) * params->max_displacement;
dp2 += Utils::sqr(dp);
// Positional increment, crop to maximum allowed by user
auto const dp = boost::algorithm::clamp(params->gamma * p.f.f[j],
-params->max_displacement,
params->max_displacement);

// Move particle
p.r.p[j] += dp;
@@ -103,9 +87,8 @@ bool steepest_descent_step(const ParticleRange &particles) {
auto const l = dq.norm();
if (l > 0.0) {
auto const axis = dq / l;
auto const angle = (std::abs(l) > params->max_displacement)
? sgn(l) * params->max_displacement
: l;
auto const angle = boost::algorithm::clamp(l, -params->max_displacement,
params->max_displacement);

// Rotate the particle around axis dq by amount l
local_rotate_particle(p, axis, angle);
@@ -116,7 +99,6 @@ bool steepest_descent_step(const ParticleRange &particles) {
#endif
// Note maximum force/torque encountered
f_max = std::max(f_max, f);
dp2_max = std::max(dp2_max, dp2);
}

set_resort_particles(Cells::RESORT_LOCAL);
@@ -126,8 +108,6 @@ bool steepest_descent_step(const ParticleRange &particles) {
auto const f_max_global =
mpi::all_reduce(comm_cart, f_max, mpi::maximum<double>());

// Return true, if the maximum force/torque encountered is below the user
// limit.
return (sqrt(f_max_global) < params->f_max);
}

@@ -24,9 +24,45 @@

#include "ParticleRange.hpp"

void minimize_energy();
/** Parameters for the steepest descent algorithm */
struct MinimizeEnergyParameters {
/** Maximal particle force
*
* If the maximal force experienced by particles in the system (in any
* direction) is inferior to this threshold, minimization stops.
*/
double f_max;
/** Dampening constant */
double gamma;
/** Maximal number of integration steps */
int max_steps;
/** Maximal particle displacement
*
* Maximal distance that a particle can travel during one integration step,
* in one direction.
*/
double max_displacement;
};

/** Steepest descent initializer
*
* Sets the parameters in @ref MinimizeEnergyParameters
*/
void minimize_energy_init(double f_max, double gamma, int max_steps,
double max_displacement);

/** Steepest descent main integration loop
*
* Integration stops when the maximal force is lower than the user limit
* @ref MinimizeEnergyParameters::f_max "f_max" or when the maximal number
* of steps @ref MinimizeEnergyParameters::max_steps "max_steps" is reached.
*/
void minimize_energy();

/** Steepest descent integrator
* @return whether the maximum force/torque encountered is below the user
* limit @ref MinimizeEnergyParameters::f_max "f_max".
*/
bool steepest_descent_step(const ParticleRange &particles);

#endif /* __MINIMIZE_ENERGY */

0 comments on commit 60f4f7a

Please sign in to comment.
You can’t perform that action at this time.