|
| 1 | +<!doctype html> |
| 2 | +<html lang="en" data-theme="light"> |
| 3 | +<head> |
| 4 | + <meta charset="UTF-8"> |
| 5 | + <meta name="viewport" content="width=device-width, initial-scale=1"> |
| 6 | + |
| 7 | + <link rel="stylesheet" href="/libs/katex/katex.min.css"> |
| 8 | + <link rel="stylesheet" href="/libs/highlight/styles/github.min.css"> |
| 9 | + |
| 10 | + |
| 11 | +<link rel="stylesheet" href="/css/franklin.css" /> |
| 12 | +<link rel="stylesheet" href="/css/icons.css" /> |
| 13 | +<link rel="icon" href="/assets/rndimg.jpg" /> |
| 14 | + |
| 15 | +<!-- Font Awesome Icons --> |
| 16 | +<link |
| 17 | + rel="stylesheet" |
| 18 | + href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.5.2/css/all.min.css" |
| 19 | +/> |
| 20 | + |
| 21 | + |
| 22 | +<link rel="stylesheet" href="/css/nav.css" /> |
| 23 | +<link rel="stylesheet" href="/css/sidebar.css" /> |
| 24 | +<link rel="stylesheet" href="/css/cards.css" /> |
| 25 | +<link rel="stylesheet" href="/css/showcase.css" /> |
| 26 | +<link rel="stylesheet" href="/css/software.css" /> |
| 27 | +<link rel="stylesheet" href="/css/posts.css" /> |
| 28 | + |
| 29 | + |
| 30 | + <title>Your Title</title> |
| 31 | +</head> |
| 32 | +<body> |
| 33 | + |
| 34 | + |
| 35 | + <div class="masthead"> |
| 36 | + <nav> |
| 37 | + <div class="author-name-nav"> |
| 38 | + <li><a href="/">Brian C. Ferrari</a></li> |
| 39 | + </div> |
| 40 | + <div class="theme-switch"> |
| 41 | + <button id="theme-toggle" onclick="toggleTheme()"> |
| 42 | + <i class="fas fa-moon"></i> |
| 43 | + </button> |
| 44 | + </div> |
| 45 | + <a class="burger-nav" onclick="myFunction()"> |
| 46 | + <i class="fa-2x fa fa-bars"></i> |
| 47 | + </a> |
| 48 | + <div class="desktop-nav"><li> |
| 49 | + <a href="/pages/aboutMe">About Me</a> |
| 50 | +</li> |
| 51 | + |
| 52 | +<li> |
| 53 | + <a href="/pages/software">Software</a> |
| 54 | +</li> |
| 55 | + |
| 56 | +<li> |
| 57 | + <a href="/pages/blog">Blog</a> |
| 58 | +</li> |
| 59 | + |
| 60 | +</div> |
| 61 | + <div id="hidden" class="mobile-nav"><li> |
| 62 | + <a href="/pages/aboutMe">About Me</a> |
| 63 | +</li> |
| 64 | + |
| 65 | +<li> |
| 66 | + <a href="/pages/software">Software</a> |
| 67 | +</li> |
| 68 | + |
| 69 | +<li> |
| 70 | + <a href="/pages/blog">Blog</a> |
| 71 | +</li> |
| 72 | + |
| 73 | +</div> |
| 74 | + </nav> |
| 75 | +</div> |
| 76 | + |
| 77 | + <div class="academic-sidebar"> |
| 78 | + <div class="author_avatar"><img |
| 79 | + src="/assets/figs/me.jpg" |
| 80 | + class="author_avatar" |
| 81 | + alt="Brian C. Ferrari" |
| 82 | + /></div> |
| 83 | + |
| 84 | + <h3>Brian C. Ferrari</h3> |
| 85 | + |
| 86 | + <!-- Pronouns are a class --> |
| 87 | + <p class="pronouns"> |
| 88 | + (he/him) |
| 89 | + </p> |
| 90 | + |
| 91 | + <!-- Bio is just text --> |
| 92 | + Chemistry PhD Candidate at Leiden University |
| 93 | + <br /><br /> |
| 94 | + |
| 95 | + <br><br> |
| 96 | +<li> |
| 97 | + <a href="https://www.linkedin.com/in/brian-ferrari-001b8a163/" target=_blank> |
| 98 | + <i class="fa fa-brands fa-linkedin" aria-hidden="true"></i> |
| 99 | + </a> |
| 100 | +</li> |
| 101 | +<li> |
| 102 | + <a href="https://orcid.org/0000-0002-7416-8629" target=_blank> |
| 103 | + <i class="fa fa-brands fa-orcid" aria-hidden="true"></i> |
| 104 | + </a> |
| 105 | +</li> |
| 106 | +<li> |
| 107 | + <a href="https://www.researchgate.net/profile/Brian-Ferrari-2?ev=hdr_xprf" target=_blank> |
| 108 | + <i class="fa fa-brands fa-researchgate" aria-hidden="true"></i> |
| 109 | + </a> |
| 110 | +</li> |
| 111 | +<li> |
| 112 | + <a href="https://github.com/Cavenfish" target=_blank> |
| 113 | + <i class="fa fa-brands fa-github" aria-hidden="true"></i> |
| 114 | + </a> |
| 115 | +</li> |
| 116 | +<li> |
| 117 | + <a href="https://scholar.google.com/citations?user=col48NcAAAAJ&hl=en" target=_blank> |
| 118 | + <i class="fa fa-graduation-cap" aria-hidden="true"></i> |
| 119 | + </a> |
| 120 | +</li> |
| 121 | + |
| 122 | +</div> |
| 123 | + |
| 124 | + |
| 125 | +<!-- Content appended here --> |
| 126 | +<div class="franklin-content"> |
| 127 | +<h1 id="canonical_velocity_rescaling_thermostat"><a href="#canonical_velocity_rescaling_thermostat" class="header-anchor">Canonical Velocity Rescaling Thermostat</a></h1> |
| 128 | +<p>This thermostat is the most widely used thermostat for molecular dynamics (MD) simulations. The <a href="https://pubs.aip.org/aip/jcp/article-abstract/126/1/014101/186581/Canonical-sampling-through-velocity-rescaling?redirectedFrom=fulltext">manuscript</a> it was published in is one of the highest cited papers related to MD simulations. Due to the massive success of this thermostat there is a wide array of resources explaining the thermostat and the math behind it in simple terms. However, none of these resources show how to implement the thermostat in a programming language. Making the source code of MD packages the only option, but they are often very difficult to follow. Since I've implemented this thermostat in Julia for my MD package (<a href="https://github.com/Cavenfish/YetAnotherSimulationSuite.jl">YASS</a>), I wanted to put out a short explanation of implementing it. </p> |
| 129 | +<h2 id="brief_math_explanation"><a href="#brief_math_explanation" class="header-anchor">Brief Math Explanation</a></h2> |
| 130 | +<p>The temperature of a system is related to the total kineteic energy of a system by,</p> |
| 131 | +\( |
| 132 | +K = \frac{1}{2} k_B N_f T |
| 133 | +\) |
| 134 | +<p>where \(K\) is the kinetic energy, \(k_B\) is Boltzmann constant, \(N_f\) is the degrees of freedom, and \(T\) is the temperature. A velocity rescaling thermostat takes advantage of this experssion by scaling the velocities to increase or decrease the temperature of the system. Typically through a scaling factor (\(\alpha\)),</p> |
| 135 | +\( |
| 136 | +\alpha = \sqrt{\frac{K_\text{target}}{K\text{current}}} |
| 137 | +\) |
| 138 | +<p>where \(K_\text{target}\) is the target kinetic energy, and \(K\text{current}\) is the kinetic energy of the current timeframe of the simulation. The CVR thermostat enforces a canonical sampling for \(K_\text{target}\) to ensure proper thermodynamic proporties in the simulation. The explicit derivation of this can be found in the original manuscript, or even other online resources. For the sake of just learning to implement the thermostat only these two concepts are necessary.</p> |
| 139 | +<h2 id="coding_it_up"><a href="#coding_it_up" class="header-anchor">Coding it Up</a></h2> |
| 140 | +<p>We first need a function that calculates the temperature of our system from the kinetic energy. We can calculate this by taking \(\sum_i m_i v_i^2 = k_B N_f T\), where the \(\frac{1}{2}\) is dropped from both sides.</p> |
| 141 | +<pre><code class="language-julia">function getTemp(v::AbstractVector, m::AbstractVector, kB::Float64) |
| 142 | + # Leave out the 1/2 to get 2Ekin for T calc |
| 143 | + N = length(m) |
| 144 | + Nf = 3N - 3 |
| 145 | + v2 = [i'i for i in v] |
| 146 | + Ekin = sum(m .* v2) |
| 147 | + |
| 148 | + Ekin / (kB * Nf) |
| 149 | +end</code></pre> |
| 150 | +<p>The next thing we need is a function for sampling from a normal and gamma distributions. You could implement this from scratch but its better to use the implementations in <a href="https://juliastats.org/Distributions.jl/stable/">Distributions.jl</a> since they are well tested and highly performant. These two functions are <code>Normal</code> and <code>Gamma</code> in the code blocks below.</p> |
| 151 | +<p>The CVR thermostat needs to have non-zero velocities to scale (or else you still have zero after scaling), so if your system has all zero velocities (or \(T=0\)) then you skip apply the thermostat and let the system gain some amount of kinetic energy due to the potential. This means adding a conditional return in your CVR function to avoid wasteful computations. In this example I use a short circuit expression (which I love but some hate), so I'll quickly explain them. These expressions take advantage of the and/or logic by using the <code>&&</code> and <code>||</code> operators to write compact if-then statements. </p> |
| 152 | +<pre><code class="language-julia"># This shows a typical if statment within a for loop |
| 153 | +for i = 1:100 |
| 154 | + if i % 10 == 0 |
| 155 | + println(i) |
| 156 | + end |
| 157 | +end |
| 158 | + |
| 159 | +# This does the same thing but using a short circuit expression |
| 160 | +for i = 1:100 |
| 161 | + i % 10 == 0 && println(i) |
| 162 | +end |
| 163 | + |
| 164 | +# Again, the same thing but using || |
| 165 | +for i = 1:100 |
| 166 | + i % 10 != 0 || println(i) |
| 167 | +end</code></pre> |
| 168 | +<p>Getting back to the CVR thermostat, once we have a non-zero kinetic energy we can rescale the velocities to try to reach our desired temperature.</p> |
| 169 | +<pre><code class="language-julia">function CVR!( |
| 170 | + v::AbstractVector, m::AbstractVector, |
| 171 | + kB::Float64, tau::Float64, Ttar::Float64 |
| 172 | +) |
| 173 | + N = length(m) |
| 174 | + Tsim = getTemp(v, m, kB) |
| 175 | + |
| 176 | + Tsim == 0.0 && return |
| 177 | + |
| 178 | + Nf = 3N - 3 |
| 179 | + n = Normal(0.0,1.0) |
| 180 | + R1 = rand(n) |
| 181 | + |
| 182 | + R2 = if (Nf-1)%2 == 0 |
| 183 | + alpha = (Nf - 1)/2 |
| 184 | + G = Gamma(alpha, 1) |
| 185 | + |
| 186 | + 2 .* rand(G) |
| 187 | + else |
| 188 | + alpha = (Nf - 2)/2 |
| 189 | + G = Gamma(alpha, 1) |
| 190 | + |
| 191 | + 2 .* rand(G) .+ (rand(n) .^2) |
| 192 | + end |
| 193 | + |
| 194 | + tau > 0.1 ? c1 = exp(-1/tau) : c1 = 0.0 |
| 195 | + |
| 196 | + K = 0.5 * Nf * kB * Tsim |
| 197 | + sigma = 0.5 * Nf * kB * Ttar |
| 198 | + |
| 199 | + Knew = @. ( |
| 200 | + K + (1-c1) * (sigma * (R2 + R1^2) / Nf - K) + |
| 201 | + 2 * R1 * sqrt(K * sigma / Nf * (1-c1) * c1) |
| 202 | + ) |
| 203 | + alpha2 = @. sqrt(Knew / K) |
| 204 | + |
| 205 | + @. v *= alpha2 |
| 206 | +end</code></pre> |
| 207 | +<div class="page-foot"> |
| 208 | + Website built with <a href="https://franklinjl.org/">Franklin.jl</a> |
| 209 | +</div> |
| 210 | +</div><!-- CONTENT ENDS HERE --> |
| 211 | + |
| 212 | + <script src="/libs/katex/katex.min.js"></script> |
| 213 | +<script src="/libs/katex/contrib/auto-render.min.js"></script> |
| 214 | +<script>renderMathInElement(document.body)</script> |
| 215 | + |
| 216 | + |
| 217 | + |
| 218 | + <script src="/libs/highlight/highlight.min.js"></script> |
| 219 | +<script>hljs.highlightAll();hljs.configure({tabReplace: ' '});</script> |
| 220 | + |
| 221 | + |
| 222 | + <script src="/libs/buttons/nav.js"></script> |
| 223 | + <script src="/libs/buttons/theme.js"></script> |
| 224 | + </body> |
| 225 | +</html> |
0 commit comments