Skip to content

Commit

Permalink
replaced LJ with WCA in charged_system-1 tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
reinaual committed Oct 11, 2019
1 parent ec1e785 commit ec2e064
Showing 1 changed file with 15 additions and 21 deletions.
36 changes: 15 additions & 21 deletions doc/tutorials/02-charged_system/02-charged_system-1.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
"#define EXTERNAL_FORCES\n",
"#define MASS\n",
"#define ELECTROSTATICS\n",
"#define LENNARD_JONES\n",
"#define WCA\n",
"```"
]
},
Expand All @@ -46,7 +46,7 @@
"plt.ion()\n",
"\n",
"# Print enabled features\n",
"required_features = [\"EXTERNAL_FORCES\", \"MASS\", \"ELECTROSTATICS\", \"LENNARD_JONES\"]\n",
"required_features = [\"EXTERNAL_FORCES\", \"MASS\", \"ELECTROSTATICS\", \"WCA\"]\n",
"espressomd.assert_features(required_features)\n",
"print(espressomd.features())\n",
"\n",
Expand All @@ -67,12 +67,8 @@
"types = {\"Anion\": 0, \"Cation\": 1}\n",
"numbers = {\"Anion\": n_ionpairs, \"Cation\": n_ionpairs}\n",
"charges = {\"Anion\": -1.0, \"Cation\": 1.0}\n",
"lj_sigmas = {\"Anion\": 1.0, \"Cation\": 1.0}\n",
"lj_epsilons = {\"Anion\": 1.0, \"Cation\": 1.0}\n",
"\n",
"WCA_cut = 2.**(1. / 6.)\n",
"lj_cuts = {\"Anion\": WCA_cut * lj_sigmas[\"Anion\"],\n",
" \"Cation\": WCA_cut * lj_sigmas[\"Cation\"]}"
"wca_sigmas = {\"Anion\": 1.0, \"Cation\": 1.0}\n",
"wca_epsilons = {\"Anion\": 1.0, \"Cation\": 1.0}"
]
},
{
Expand Down Expand Up @@ -139,12 +135,11 @@
"metadata": {},
"source": [
"Before we can really start the simulation, we have to specify the\n",
"interactions between our particles. We already defined the Lennard-Jones parameters at the beginning,\n",
"interactions between our particles. We already defined the WCA parameters at the beginning,\n",
"what is left is to specify the combination rule and to iterate over particle type pairs. For simplicity, \n",
"we implement only the *Lorentz-Berthelot* rules. \n",
"We pass our interaction pair to <tt>system.non_bonded_inter[\\*,\\*]</tt> and set the \n",
"pre-calculated LJ parameters <tt>epsilon</tt>, <tt>sigma</tt> and <tt>cutoff</tt>. With <tt>shift=\"auto\"</tt>,\n",
"we shift the interaction potential to the cutoff so that $U_\\mathrm{LJ}(r_\\mathrm{cutoff})=0$."
"pre-calculated WCA parameters <tt>epsilon</tt> and <tt>sigma</tt>."
]
},
{
Expand All @@ -165,14 +160,13 @@
" else:\n",
" return ValueError(\"No combination rule defined\")\n",
"\n",
"# Lennard-Jones interactions parameters\n",
"# WCA interactions parameters\n",
"for s in [[\"Anion\", \"Cation\"], [\"Anion\", \"Anion\"], [\"Cation\", \"Cation\"]]:\n",
" lj_sig = combination_rule_sigma(\"Berthelot\", lj_sigmas[s[0]], lj_sigmas[s[1]])\n",
" lj_cut = combination_rule_sigma(\"Berthelot\", lj_cuts[s[0]], lj_cuts[s[1]])\n",
" lj_eps = combination_rule_epsilon(\"Lorentz\", lj_epsilons[s[0]], lj_epsilons[s[1]])\n",
" wca_sig = combination_rule_sigma(\"Berthelot\", wca_sigmas[s[0]], wca_sigmas[s[1]])\n",
" wca_eps = combination_rule_epsilon(\"Lorentz\", wca_epsilons[s[0]], wca_epsilons[s[1]])\n",
"\n",
" system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(\n",
" epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift=\"auto\")"
" system.non_bonded_inter[types[s[0]], types[s[1]]].wca.set_params(\n",
" epsilon=wca_eps, sigma=wca_sig)"
]
},
{
Expand All @@ -182,7 +176,7 @@
"## 3 Equilibration\n",
"\n",
"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",
"cause the simulation to crash. The next step in our script therefore is a suitable WCA 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 the steepest descent algorithm and cap the maximal particle displacement per integration step\n",
"to 1% of sigma.\n",
Expand All @@ -197,8 +191,8 @@
"metadata": {},
"outputs": [],
"source": [
"# Lennard-Jones Equilibration\n",
"max_sigma = max(lj_sigmas.values())\n",
"# WCA Equilibration\n",
"max_sigma = max(wca_sigmas.values())\n",
"min_dist = 0.0\n",
"system.minimize_energy.init(f_max=0, gamma=10.0, max_steps=10,\n",
" max_displacement=max_sigma * 0.01)\n",
Expand Down Expand Up @@ -423,7 +417,7 @@
"\n",
"To make your life easier, don't try to equilibrate randomly positioned particles,\n",
"but set them up in a crystal structure close to equilibrium. If you do it right,\n",
"you don't even need the Lennard-Jones equilibration. \n",
"you don't even need the WCA equilibration. \n",
"To speed things up, don't go further than 1000 particles and use a P$^3$M accuracy of $10^{-2}$.\n",
"Your RDF should look like the plot in figure 2. If you get stuck,\n",
"you can look at the solution script <tt>/doc/tutorials/02-charged_system/scripts/nacl_units.py</tt> (or <tt>nacl_units_vis.py</tt> with visualization)."
Expand Down

0 comments on commit ec2e064

Please sign in to comment.