Skip to content

Commit a9eafba

Browse files
committed
fixed umbrella sampling part
1 parent 3a6f5f6 commit a9eafba

2 files changed

Lines changed: 44 additions & 47 deletions

File tree

files/tutorial7/solution/umbrella-sampling.lmp

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@ variable epsilon equal 0.238
44
variable U0 equal 10*${epsilon}
55
variable dlt equal 1.0
66
variable x0 equal 10.0
7-
87
variable k equal 1.5
98

109
units real

lammps-tutorials.tex

Lines changed: 44 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -3764,6 +3764,7 @@ \subsubsection{Method 1: Free sampling}
37643764
In such case, employing an enhanced sampling method is recommended, as done in the next section.
37653765

37663766
\subsubsection{Method 2: Umbrella sampling}
3767+
37673768
Umbrella sampling is a biased molecular dynamics method in which additional forces
37683769
are added to a chosen atom to force it to explore the more unfavorable areas of
37693770
the system \cite{kastner2011umbrella, allen2017computer, frenkel2023understanding}.
@@ -3793,9 +3794,11 @@ \subsubsection{Method 2: Umbrella sampling}
37933794
pair_modify shift yes
37943795
boundary p p p
37953796
\end{lstlisting}
3796-
The only difference with the previous case is the larger value
3797+
The first difference with the previous case is the larger value
37973798
for the repulsive potential $U_0$, which will make the central area
3798-
of the system very unliky to be visited by free particles.
3799+
of the system very unliky to be visited by free particles. The second
3800+
difference is the addition of the $k$ variable that will be used for
3801+
the biasing potential.
37993802

38003803
Let us create a box with 2 atom types, with a single particle of type 2,
38013804
by adding the following lines to \lmpcmd{umbrella-sampling.lmp}:
@@ -3839,6 +3842,12 @@ \subsubsection{Method 2: Umbrella sampling}
38393842
run 50000
38403843
\end{lstlisting}
38413844

3845+
So far, our code resembles that of Method 1, except for the additional particle
3846+
of type 2. Particles of type 1 and 2 are identical, having the same mass and
3847+
Lennard-Jones parameters. However, the particle of type 2 will additionally
3848+
be exposed to the biasing potential $V$, forcing it to explore
3849+
the central part of box (Fig.\,\ref{fig:US-system-biased}).
3850+
38423851
\begin{figure}
38433852
\centering
38443853
\includegraphics[width=\linewidth]{US-system-biased}
@@ -3851,59 +3860,48 @@ \subsubsection{Method 2: Umbrella sampling}
38513860
\label{fig:US-system-biased}
38523861
\end{figure}
38533862

3854-
\clearpage
3855-
3856-
3857-
\begin{lstlisting}
3858-
3859-
fix mynve all nve
3860-
fix mylgv all langevin 119.8 119.8 50 1530917
3861-
timestep 2.0
3862-
thermo 100000
3863-
run 500000
3864-
reset_timestep 0
3865-
3866-
dump mydmp all image 1000000 dump.*.jpg type &
3867-
type shiny 0.1 box yes 0.02 view 0 90 zoom 1.9
3868-
dump_modify mydmp backcolor white &
3869-
acolor A1 cyan adiam A1 1 &
3870-
acolor A2 pink adiam A2 1 &
3871-
boxcolor black
3872-
\end{lstlisting}
3873-
So far, this code resembles that of Method 1, except for the additional particle
3874-
of type A2. Particles of type A1 and A2
3875-
are identical, having the same mass and Lennard-Jones parameters. However, the
3876-
particle of type A2 will additionally be exposed to the biasing potential
3877-
$V$ (Fig.\,){fig:US-system-biased}).
3878-
3879-
Let us create a loop with 50 steps, and move progressively the center of the
3880-
bias potential by an increment of 0.1 nm. Add the following lines to \flecmd{input.lmp}:
3863+
Now, let us create a loop with 50 steps, and move progressively the center of the
3864+
bias potential by an increment of 0.2 nm. Add the following lines to \flecmd{umbrella-sampling.lmp}:
38813865
\begin{lstlisting}
38823866
variable a loop 50
38833867
label loop
3884-
variable xdes equal ${a}-25
3868+
3869+
variable xdes equal 2*${a}-50
38853870
variable xave equal xcm(topull,x)
3886-
fix mytth topull spring tether ${k} &
3887-
${xdes} 0 0 0
3888-
run 200000
3889-
fix myat1 all ave/time 10 10 100 v_xave &
3890-
v_xdes file data/position.${a}.dat
3891-
run 1000000
3871+
fix mytth topull spring tether ${k} ${xdes} 0 0 0
3872+
3873+
run 20000
3874+
3875+
fix myat1 all ave/time 10 10 100 &
3876+
v_xave v_xdes file umbrella-sampling.${a}.dat
3877+
3878+
run 100000
38923879
unfix myat1
38933880
next a
38943881
jump SELF loop
38953882
\end{lstlisting}
3896-
A folder called \flrcmd{data/} needs to be created within \flrcmd{BiasedSampling/}.
38973883
The spring command serves to impose the additional harmonic potential $V$ with
3898-
the spring constant $k$. Note that the value of $k$ should be chosen with care,
3899-
if $k$ is too small, the particle won't follow the biasing potential center,
3900-
and if $k$ is too large, there will be no overlapping between the different windows.
3901-
The center of the harmonic potential $x_\text{des}$ successively takes values
3902-
from -25 to 25. For each value of $x_\text{des}$, an equilibration step of 0.4 ns
3903-
is performed, followed by a step of 2 ns during which the position along $x$ of
3904-
the particle is saved in data files (one data file per value of $x_\text{des}$).
3905-
You can always increase the duration of the runs for better samplings. Run the
3906-
\flecmd{input.lmp} file using LAMMPS.
3884+
the previously defined spring constant $k$. The center of the harmonic
3885+
potential, $x_\text{des}$, successively takes values
3886+
from -48 to 50. For each value of $x_\text{des}$, an equilibration step of 40\,ps
3887+
is performed, followed by a step of 200\,ps during which the position along $x$ of
3888+
the particle of type 2, $x_\text{ave}$, is saved in data files named \flecmd{umbrella-sampling.i.dat},
3889+
where i goes from 1 to 50. Run the \flecmd{umbrella-sampling.lmp} file using LAMMPS.
3890+
3891+
\begin{note}
3892+
The value of $k$ should be chosen with care:
3893+
if $k$ is too small the particle won't follow the biasing potential,
3894+
and if $k$ is too large there will be no overlapping between
3895+
the different windows, leading to poor reconstruction of the free energy profile.
3896+
\end{note}
3897+
3898+
3899+
3900+
\clearpage
3901+
3902+
3903+
3904+
39073905

39083906

39093907

0 commit comments

Comments
 (0)