Skip to content

Commit 57698d4

Browse files
committed
updated file
1 parent 6d56974 commit 57698d4

4 files changed

Lines changed: 39 additions & 49 deletions

File tree

figures/US-system-unbiased.png

-28.7 KB
Loading

files/tutorial7/solution/free-sampling.lmp

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -11,26 +11,26 @@ pair_style lj/cut 3.822
1111
pair_modify shift yes
1212
boundary p p p
1313

14-
region myreg block 25 25 5 5 25 25
14+
region myreg block -25 25 -5 5 -25 25
1515
create_box 1 myreg
1616
create_atoms 1 random 200 34134 myreg overlap 3 maxtry 50
1717

1818
mass * 39.95
1919
pair_coeff * * ${epsilon} ${sigma}
2020

21-
variable U atom ${U0}*atan((x+${x0})/${dlt})-${U0}*atan((x${x0})/${dlt})
22-
variable F atom ${U0}/((x${x0})^2/${dlt}^2+1)/${dlt}-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}
21+
variable U atom ${U0}*atan((x+${x0})/${dlt})-${U0}*atan((x-${x0})/${dlt})
22+
variable F atom ${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt}-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}
2323
fix myadf all addforce v_F 0.0 0.0 energy v_U
2424

2525
fix mynve all nve
26-
fix mylgv all langevin 119.8 119.8 50 1530917
26+
fix mylgv all langevin 119.8 119.8 50 30917
2727

28-
region mymes block ${x0} ${x0} INF INF INF INF
28+
region mymes block -${x0} ${x0} INF INF INF INF
2929
variable n_center equal count(all,mymes)
3030
thermo_style custom step temp etotal v_n_center
3131
thermo 10000
3232

33-
dump viz all image 50000 myimage*.ppm type type &
33+
dump viz all image 50000 myimage-*.ppm type type &
3434
shiny 0.1 box yes 0.01 view 180 90 zoom 7 size 1600 400
3535
dump_modify viz backcolor white acolor 1 cyan adiam 1 1 boxcolor black
3636

lammps-tutorials.tex

Lines changed: 32 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -3618,7 +3618,7 @@ \subsubsection{Method 1: Free sampling}
36183618
Taking the derivative of the potential with respect to $x$, we obtain the expression
36193619
for the force that will be imposed on the atoms:
36203620
\begin{equation}
3621-
F= \dfrac{U_0}{\delta} \left[ \dfrac{1}{(x-x_0)^2/\delta^2+1}
3621+
F = \dfrac{U_0}{\delta} \left[ \dfrac{1}{(x-x_0)^2/\delta^2+1}
36223622
- \dfrac{1}{(x+x_0)^2/\delta^2+1} \right].
36233623
\label{eq:F}
36243624
\end{equation}
@@ -3638,80 +3638,70 @@ \subsubsection{Method 1: Free sampling}
36383638
\label{fig:potential}
36393639
\end{figure}
36403640

3641-
Let us apply energy minimization to the system, and then impose the force $F(x)$
3642-
to all of the atoms in the simulation using the \textit{addforce} command. Add
3643-
the following lines to \flecmd{input.lmp}:
3641+
Let us impose the force $F(x)$ to the atoms in the simulation
3642+
using the \lmpcmd{fix addforce} command. Add the following
3643+
lines to \flecmd{free-energy.lmp}:
36443644
\begin{lstlisting}
3645-
minimize 1e-4 1e-6 100 1000
3646-
reset_timestep 0
3647-
36483645
variable U atom ${U0}*atan((x+${x0})/${dlt})&
36493646
-${U0}*atan((x-${x0})/${dlt})
3650-
variable F atom &
3651-
${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt}&
3647+
variable F atom ${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt}&
36523648
-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}
36533649
fix myadf all addforce v_F 0.0 0.0 energy v_U
36543650
\end{lstlisting}
3655-
Finally, let us combine the \textit{fix nve} with a \textit{Langevin} thermostat
3656-
and run a molecular dynamics simulation:
3651+
Let us combine the \lmpcmd{fix nve} with a \lmpcmd{fix langevin} thermostat:
36573652
\begin{lstlisting}
36583653
fix mynve all nve
3659-
fix mylgv all langevin 119.8 119.8 50 1530917
3654+
fix mylgv all langevin 119.8 119.8 50 30917
36603655
\end{lstlisting}
3661-
With these two commands, the MD simulation
3662-
is effectively in the NVT ensemble: constant number of atoms $N$, constant volume
3656+
When combining these two commands, the MD simulation is effectively
3657+
in the NVT ensemble: constant number of atoms $N$, constant volume
36633658
$V$, and constant temperature $T$.
36643659

3665-
To make sure that $1\,\text{ns}$ is long enough, we will record the evolution of
3666-
the number of atoms in the central (energetically unfavorable) region called \textit{mymes}
3667-
using the \textit{n\_center} variable:
3660+
To make sure that the equilibration time is sufficient, we will record the evolution of
3661+
the number of atoms in the central (energetically unfavorable) region called \lmpcmd{mymes}
3662+
using the \lmpcmd{n\_center} variable:
36683663
\begin{lstlisting}
36693664
region mymes block -${x0} ${x0} INF INF INF INF
36703665
variable n_center equal count(all,mymes)
36713666
thermo_style custom step temp etotal v_n_center
36723667
thermo 10000
36733668

3674-
dump mydmp all image 5000 dump.*.jpg type type &
3675-
shiny 0.1 box yes 0.02 view 0 90 zoom 1.9
3676-
dump_modify mydmp backcolor white &
3677-
acolor A1 cyan adiam A1 1 boxcolor black
3669+
dump viz all image 50000 myimage-*.ppm type type &
3670+
shiny 0.1 box yes 0.01 view 180 90 zoom 7 size 1600 400
3671+
dump_modify viz backcolor white acolor 1 cyan adiam 1 1 &
3672+
boxcolor black
36783673
\end{lstlisting}
3679-
A \textit{dump image} was also added for system visualization.
3674+
A \lmpcmd{dump image} was also added for system visualization.
36803675

36813676
Finally, let us perform an equilibration of 500000 steps
36823677
in total, using a timestep of $2\,\text{fs}$ (i.e.~a total duration of $1\,\text{ns}$):
36833678
\begin{lstlisting}
36843679
timestep 2.0
36853680
run 500000
36863681
\end{lstlisting}
3687-
Run the simulation with LAMMPS. To ensure that the equilibration of $1\,\text{ns}$ is long
3688-
enough, let us have a look at the evolution of the number of atoms in the central region,
3689-
$n_\mathrm{center}$. We can see that $n_\mathrm{center}$ reaches
3690-
its equilibrium value (which is close to 0) after about $0.1\,\text{ns}$
3691-
(Fig.\,\ref{fig:US-density-evolution}). See also a snapshot of the equilibrated
3692-
system in Fig.\,\ref{fig:US-system-unbiased}.
3682+
Run the simulation with LAMMPS. The number of atoms in the
3683+
central region, $n_\mathrm{center}$, reaches
3684+
its equilibrium value after approximately $0.1\,\text{ns}$
3685+
(Fig.\,\ref{fig:US-density-evolution}). A snapshot of the
3686+
equilibrated system is shown in Fig.\,\ref{fig:US-system-unbiased}.
36933687

36943688
\paragraph{Run and data acquisition}
3695-
Once the system is equilibrated, let us record the density profile of
3696-
the atoms along the $x$ axis using
3697-
the \textit{ave/chunk} command. A total of 10 density profiles will be printed.
3698-
Add the following line to \flecmd{input.lmp}:
3689+
3690+
Once the system is equilibrated, we will record the density profile of
3691+
the atoms along the $x$ axis using the \lmpcmd{ave/chunk} command. A total
3692+
of 10 density profiles will be printed. Add the following line to \flecmd{free-energy.lmp}:
36993693
\begin{lstlisting}
3700-
unfix myat
3701-
undump mydmp
37023694
reset_timestep 0
37033695

37043696
compute cc1 all chunk/atom bin/1d x 0.0 1.0
37053697
fix myac all ave/chunk 10 400000 4000000 &
3706-
cc1 density/number file density_profile.dat
3707-
dump mydmp all atom 200000 dump.lammpstrj
3698+
cc1 density/number file free-sampling.dat
37083699

3709-
thermo 100000
37103700
run 4000000
37113701
\end{lstlisting}
3712-
The step count is reset to 0 using \textit{reset\_timestep} to synchronize
3713-
with the output times of \textit{fix density/number}. Feel free to increase the
3714-
duration of the last run for a better resolved density profile.
3702+
The step count is reset to 0 using \lmpcmd{reset\_timestep} to synchronize it
3703+
with the output times of \lmpcmd{fix density/number}. Increase the
3704+
duration of the last run to obtain a better-resolved density profile.
37153705

37163706
\begin{figure}
37173707
\centering
@@ -3740,8 +3730,8 @@ \subsubsection{Method 1: Free sampling}
37403730
\begin{figure}
37413731
\centering
37423732
\includegraphics[width=\linewidth]{US-system-unbiased}
3743-
\caption{Snapshot of the system. The density of atoms is lowest in the central
3744-
part of the box, \textit{mymes}, due to the additional force $F$. Here,
3733+
\caption{Snapshot of the system. The density of atoms is lowest in the central
3734+
part of the box, \lmpcmd{mymes}, due to the additional force $F$. Here,
37453735
$U_0 = 0.36~\text{kcal/mol}$, $\delta = 0.5~\mathrm{\AA{}}$, and $x_0 = 5~\mathrm{\AA{}}$.}
37463736
\label{fig:US-system-unbiased}
37473737
\end{figure}

0 commit comments

Comments
 (0)