-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathparticle_simulation.jl
2702 lines (2221 loc) · 92.7 KB
/
particle_simulation.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.16.0
using Markdown
using InteractiveUtils
# ╔═╡ a82f22e4-f35b-461a-b481-1dff43722e44
using StaticArrays
# ╔═╡ d42f842d-6c2a-40db-b0c4-e936244a9e7c
using BenchmarkTools
# ╔═╡ 99b5c818-a825-4939-849e-1cade802f63d
using Measurements
# ╔═╡ d6564250-f646-40de-9463-a956af1a5b1d
using ForwardDiff
# ╔═╡ d2f4d622-8e35-4be3-b421-39b28a748cab
using CellListMap
# ╔═╡ f049ab19-7ecf-4c65-bf6d-21352c1fe767
using FastPow
# ╔═╡ 7c792b6b-b6ee-4e30-88d5-d0b8064f2734
begin
using Plots
plot_font = "Computer Modern"
default(
fontfamily=plot_font,
linewidth=2, framestyle=:box, label=:none, grid=false,
size=(400,400)
)
end
# ╔═╡ febe8c06-b3aa-4db1-a3ea-fdc2a81bdebd
using Printf
# ╔═╡ a756dd18-fac6-4527-944e-c16d8cc4bf95
begin
using PlutoUI
TableOfContents()
end
# ╔═╡ a87fad48-73c1-4a08-a6f1-aae759b3c6fc
md"""
# Particle Simulations with Julia
Leandro Martínez
Institute of Chemistry - University of Campinas
[http://m3g.iqm.unicamp.br](http://m3g.iqm.unicamp.br) -
[https://github.com/m3g](https://github.com/m3g)
"""
# ╔═╡ 172227c2-b27a-40db-91f4-9566c2f6cf52
md"""
# Outline
- Elements of a particle simulation
- Benchmarking vs. a conventional compiled language (Fortran)
- Exploring the generic character of functions
- Differentiable simulations and parameter fitting
- Using cell lists
- An efficient and generic cell list implementation
- The Packmol strategy
- Benchmarking vs. NAMD
- Remarks
"""
# ╔═╡ 4b484cf6-4888-4f04-b3fd-94862822b0c0
md"""
# Defining the type of particle
We define a simple vector in 2D space, with coordinates `x` and `y`. The vector will be defined with the aid of the `StaticArrays` package, which provides convenient constructors for this type of variable, and all the arithmetics. The memory layout of a vector of these vectors is identical to that of a `N×M` matrix, where `N` is the dimensio nf the space (2D here) and `M` is the number of vector. Julia is column-major, thus this is the most efficient memory layout for this type of computation.
"""
# ╔═╡ 8c444ee4-8c77-413a-bbeb-9e5ae2428876
struct Vec2D{T} <: FieldVector{2,T}
x::T
y::T
end
# ╔═╡ a0de01b5-779a-48c0-8d61-12b02a5f527e
md"""
For convenience, here we will also define a function that returns a random vector, given a range of coordinates:
"""
# ╔═╡ 532eb3bc-5522-4348-afa5-5336ec6752c7
md"""
In defining the function above we took care of making it generic for the type and dimension of the vectors desired, such that we do not need to redefine it later when peforming simulations with different data structures.
"""
# ╔═╡ dc5f7484-3dc3-47a7-ad4a-30f97fc14d11
md"""
## Force between a pair of particles
Initially, the energy function will be a soft potential, which is zero for distances greater than a cutoff, and increasing quadratically for distances smaller than the cutoff:
If $d = ||\vec{y}-\vec{x}||$ is the norm of the relative position of two positions, we have:
$$u(\vec{x},\vec{y},c)=
\begin{cases}
(d-c)^2 &\textrm{if} & d\leq c \\
0 & \textrm{if} & d > c \\
\end{cases}$$
for which the forces are
$$\vec{f_x}(\vec{x},\vec{y},c)=
\begin{cases}
2(d-c)\frac{(\vec{y}-\vec{x})}{d} &\textrm{if} & d\leq c \\
\vec{0} & \textrm{if} & d > c \\
\end{cases}$$
and
$$\vec{f_y} = -\vec{f_x}$$.
"""
# ╔═╡ ab3ff21b-bf82-4d8c-abd1-c27418956ed8
md"""
The standard Julia `LinearAlgebra` library provides a `norm` function, and there is no reason not to use it (although a manual definition of the same function can also be easily implemented):
"""
# ╔═╡ 7a1db355-bba9-4322-9fb4-a6d7b7bdd60d
import LinearAlgebra: norm
# ╔═╡ cc49ef04-08d8-42bb-9170-9db64e275a51
md"""
The energy and force functions are clear to read:
"""
# ╔═╡ d00e56f2-9d3a-4dd3-80eb-3201eff39b96
md"""
And for a unidimensional case, with a defined cutoff, look like:
"""
# ╔═╡ b5c09cd3-6063-4a36-96cd-2d128aa11b82
const cutoff = 5.
# ╔═╡ 7719b317-e85b-4583-b401-a8614d4b2373
md"""
The function that will compute the force over all pairs will just *naively* run over all (non-repeated) the pairs. The function `forces!` will receive as a parameter the function that computes the force between pairs, such that this pairwise function can be changed later.
Inside `forces!`, the `force_pair` function will receive four parameters: the indexes of the particles and their positions. We will use the indexes later.
"""
# ╔═╡ 144119ad-ab88-4165-883a-f2fc2464a838
md"""
Let us create some particle positions to explain how the function will be called.
"""
# ╔═╡ dd9e4332-2908-40ef-b461-6b571df56cf4
md"""
The function `force_pair`, will be passed to the function that computes the forces to all pairs as *closure*, which will capture the value of the cutoff. The closure also allows us to ignore the indexes of the particles, which are expected by the inner implementation of the function inside `forces`. For example:
"""
# ╔═╡ 017cd6a8-712d-4ec5-b31f-7f1f507764f2
md"""
The third argument of `forces!` function is a *closure*, which can be read as: it is the function that *given* `(i,j,x,y)`, returns `fₓ(x,y,cutoff)`. Thus, it is consistent with the internal call of `fₓ` of `forces!`, and *closes over* the additional parameter `cutoff` required for the computation.
"""
# ╔═╡ b5206dd5-1f46-4437-929b-efd68393b12b
md"""
# Performing a particle simulation
Now, given the function that computes the forces, we can perform a particle simulation. We will use a simple Euler integration scheme, and the algorithm will be:
1. Compute forces at time $t$ from positions $x$:
$f(t) = f(x)$
2. Update the positions (using $a = f/m$):
$x(t + dt) = x(t) + v(t)dt + a(t)dt^2/2$
3. Update the velocities:
$v(t+dt) = v(t) + a(t)dt$
4. Goto 1.
## The actual simulation code is as short:
"""
# ╔═╡ eb5dc224-1491-11ec-1cae-d51c93cd292c
function md(
x0::Vector{T},
v0::Vector{T},
mass,dt,nsteps,isave,forces!
) where T
x = copy(x0)
v = copy(v0)
a = similar(x0)
f = similar(x0)
trajectory = [ copy(x0) ] # will store the trajectory
for step in 1:nsteps
# Compute forces
forces!(f,x)
# Accelerations
@. a = f / mass
# Update positions
@. x = x + v*dt + a*dt^2/2
# Update velocities
@. v = v + a*dt
# Save if required
if mod(step,isave) == 0
println("Saved trajectory at step: ",step)
push!(trajectory,copy(x))
end
end
return trajectory
end
# ╔═╡ 594ba1d6-2dae-4f20-9546-f52fac17c2f0
md"""
By using a parametric type of input (i. e. `Vector{T}`) we can guarantee that an error will be thrown if the positions and velocities are not provided as the same type of variable.
The `@.` notation is very common in Julia and means that the computation will be performed element-wise.
"""
# ╔═╡ 66c7d960-7e05-4613-84e8-2a40fe40dc3d
md"""
## Let us run the simulation!
"""
# ╔═╡ e717a8d9-ccfb-4f89-b2a2-f244f108b48d
md"""
Here we generate random positions and velocities, and use masses equal to `1.0` for all particles.
"""
# ╔═╡ eab7b195-64d5-4587-8687-48a673ab091b
md"""
## Using periodic boundary conditions
Our particles just explode, since they have initial random velocities and there are only repulsive interactions.
We have a more interesting dynamics if we use periodic boundary conditions. To do so, we will update how the forces are computed.
"""
# ╔═╡ 34dc72dc-4864-47c0-b730-183f67e7aea3
md"""
## Wrapping of coordinates
The following function defines how to wrap the coordinates on the boundaries, for a square or cubic box of side `side`:
"""
# ╔═╡ 02d9bf3b-708c-4293-b198-9043b334ff7e
md"""
This allows writting the force computation now as:
"""
# ╔═╡ 0a5282ee-c88a-4bcc-aca2-477f28e9e04d
md"""
Our box has a side of 100:
"""
# ╔═╡ b2a4a505-47ff-40bb-9a6d-a08d91c53217
const side = 100.
# ╔═╡ fcff6973-012a-40fc-a618-f6262266287a
md"""
To run the simulation with the new periodic forces, we use the same `md` function, just passing the new `fₓ` function in the *closure* definition:
"""
# ╔═╡ 14867ffd-cde5-43f8-8399-01169ee29b73
md"""
A relevant detail here is that we could use the same `fₓ` name for the function, because it receives the `side` of the box as a parameter, and multiple-dispatch then chooses the correct method automaticaly.
"""
# ╔═╡ c4798182-de75-4b59-8be7-f7cf1051364d
md"""
While plotting the trajectory, we will wrap the coordinates:
"""
# ╔═╡ 22fb0386-d4fa-47b9-ac31-decf2731cbc1
md"""
## Benchmarking
"""
# ╔═╡ 8e23a3ea-3039-4a5f-b37f-c4710153938e
md"""
Benchmarkming in Julia can be done with the `@time` macro or the macros from the `BenchmarkTools` package. Compilation occurs on the first call to each method, and the macros from `BenchmarkTools` discount the compilation time automatically.
"""
# ╔═╡ 2a3e7257-63ad-4761-beda-cec18b91f99c
md"""
Something of the order of `200ms` and `200KiB` of allocations does not seem bad, but it doesn't mean anything either. What is interesting to point here is just that this code, compared to ahead-of-time compiled language like Fortran, is completely comparable in terms of performance, as [this benchmark](https://github.com/m3g/2021_FortranCon/tree/main/benchmark_vs_fortran) shows.
"""
# ╔═╡ 49b1f040-929a-4238-acd9-6554757b592c
md"""
# Exploring generics
## Running the simulations in 3D
Not much is needed to just run the simulation in three dimensions. We only need to define our 3D vector:
"""
# ╔═╡ 26d5c6e9-a903-4792-a0e0-dec1a2e86a01
struct Vec3D{T} <: FieldVector{3,T}
x::T
y::T
z::T
end
# ╔═╡ 2aef27bd-dea6-4a93-9d0f-b9249c9dd2cd
md"""
That is enough such that we can run the simulations in 3D:
"""
# ╔═╡ b6dcb9a3-59e3-4eae-9399-fb072c704f1a
md"""
## Automatic error propagation
Performing simulations in different dimensions is not the most interesting, or most useful property of generic programming. We can, for example, propagate the error in the positions of the particles, simply by defining a type of particle that carries both the position and the cumulative error.
A small example of how that can be done is shown. First, we create a type of variable that carries both the coordinates and the uncertainty on the coordinates:
"""
# ╔═╡ e4657169-1bb2-4d4a-ac9d-adc80499d07d
struct MyMeasurement{T}
x::T
Δx::T
end
# ╔═╡ 5d395353-5681-4780-983e-902fdb89eaf2
md"""
and we will overload the printing of this variables to make things prettier:
"""
# ╔═╡ e9376a4b-3d60-42eb-8681-cd2bcec13fe8
Base.show(io::IO,m::MyMeasurement) = println(io," $(m.x) ± $(m.Δx)")
# ╔═╡ c5cdd71f-5ded-482f-9205-c13f01a14d0b
m = MyMeasurement(1.0,0.1)
# ╔═╡ a0993a1f-60a6-45b5-815e-676c49a9f049
md"""
Now we define the arithmetics for this type of variable. For example, the sum of two `MyMeasurement`s sums the uncertainties, but so do the subtraction. The other uncertainties are also propagaged linearly, according to the first derivative of their operations relative to the values:
"""
# ╔═╡ 4e7f8db4-b5cc-4a3e-9fa7-e62d8f2a36ac
begin
import Base: -, +, *, /, ^, sqrt
+(m1::MyMeasurement,m2::MyMeasurement) = MyMeasurement(m1.x+m2.x,m1.Δx+m2.Δx)
-(m1::MyMeasurement,m2::MyMeasurement) = MyMeasurement(m1.x-m2.x,m1.Δx+m2.Δx)
*(α,m::MyMeasurement) = MyMeasurement(α*m.x,sign(α)*α*m.Δx)
*(m::MyMeasurement,α) = α*m
/(m::MyMeasurement,α) = inv(α)*m
sqrt(m::MyMeasurement{T}) where T = MyMeasurement(sqrt(m.x),inv(2*sqrt(m.x))*m.Δx)
^(m::MyMeasurement{T},n) where T = MyMeasurement{T}(m.x^n,n*m.x^(n-1)*m.Δx)
end
# ╔═╡ f87e4036-8f82-41c7-90c1-daa5f677488d
function random_vec(::Type{VecType},range) where VecType
dim = length(VecType)
T = eltype(VecType)
p = VecType(
range[begin] + rand(T)*(range[end]-range[begin]) for _ in 1:dim
)
return p
end
# ╔═╡ df33b999-4a42-4133-bf59-5a65240790cf
function energy(x::T,y::T,cutoff) where T
Δv = y - x
d = norm(Δv)
if d > cutoff
energy = zero(T)
else
energy = (d - cutoff)^2
end
return energy
end
# ╔═╡ 0f52365d-34f4-46ed-923e-3ea31c6db0ca
function fₓ(x::T,y::T,cutoff) where T
Δv = y - x
d = norm(Δv)
if d > cutoff
fₓ = zero(T)
else
fₓ = 2*(d - cutoff)*(Δv/d)
end
return fₓ
end
# ╔═╡ f58769a6-a656-42a3-8bc6-c204d4cfd897
function forces!(f::Vector{T},x::Vector{T},fₓ::F) where {T,F}
fill!(f,zero(T))
n = length(x)
for i in 1:n-1
for j in i+1:n
fᵢ = fₓ(i,j,x[i],x[j])
f[i] += fᵢ
f[j] -= fᵢ
end
end
return f
end
# ╔═╡ beeb3335-5c49-47de-a1d3-3eef5f9479f1
function wrap(x,side)
x = rem(x,side)
if x >= side/2
x -= side
elseif x < -side/2
x += side
end
return x
end
# ╔═╡ 0967b90d-ac88-476d-a57a-7c38dfa82204
function fₓ(x::T,y::T,cutoff,side) where T
Δv = wrap.(y - x, side)
d = norm(Δv)
if d > cutoff
fₓ = zero(T)
else
fₓ = 2*(d - cutoff)*(Δv/d)
end
return fₓ
end
# ╔═╡ 36da3e92-000c-4d4b-9abf-4cd588b3a354
md"""
With such definitions, we can operate over variables of type `MyMeasurement`, propagating automatically the uncertainty along the operations:
"""
# ╔═╡ 70eb2e0a-a5c8-4975-8f6c-589035bea29c
sqrt((2*(m + 4*m)^2/3))
# ╔═╡ b4646a29-3efd-4bd1-bffc-3575559de937
md"""
We can also define a 2D (or 3D) vectors of values with uncertainties, without changing the previous definitions of these:
"""
# ╔═╡ d32743c0-fc80-406f-83c5-4528e439589a
x = Vec2D(MyMeasurement(1.0,0.1),MyMeasurement(2.0,0.2))
# ╔═╡ 98478246-5940-4828-a8f1-9c9fa990676d
md"""
Now operations on this vector propagate the uncertainties of the componentes as well:
"""
# ╔═╡ 310f247e-3fe8-4621-ae0b-b5ee38d2ee89
2*x .+ sqrt.(x)
# ╔═╡ 8267220a-f06e-4761-b310-00f8ba44e4b1
md"""
Propagating uncertainties in more general scenarios requires the definition of other propagation rules. Also, one might want to consider the correlation between variables, which makes the propagation rules more complicated and expensive.
Fortunately, there are some package that provide the error propagation in more general scenarios, by defining the proper progagation rules.
Here, we use the `Measurements` package.
"""
# ╔═╡ ce916139-221a-462e-877f-88212663c05e
md"""
### Using `Measurements`
"""
# ╔═╡ 8e2903be-4975-4e14-84ed-6e712f47fe47
md"""
Using `Measurments` we do not need to change anything in our previous code, but only redefine the content of our vectors, which will now carry in each coordinate the position and the error in the position, accumulated from an initial uncertainty:
"""
# ╔═╡ 418f31bb-81d5-459b-b402-4fd4e3f4ab27
md"""
We need to redefine your initial random vector generator only:
"""
# ╔═╡ 05402cbd-78c6-4234-8680-c351c8c37778
function random_vec(::Type{Vec2D{Measurement{T}}},range,Δ) where T
p = Vec2D(
range[begin] + rand(T)*(range[end]-range[begin]) ± rand()*Δ,
range[begin] + rand(T)*(range[end]-range[begin]) ± rand()*Δ
)
return p
end
# ╔═╡ 356ac5a4-c94e-42cb-a085-0198b29c7e52
x0 = [ random_vec(Vec2D{Float64},(0,100)) for _ in 1:100]
# ╔═╡ d23b4a92-055e-4ed7-bd46-8a3c59312993
f = similar(x0)
# ╔═╡ e6e29d1e-9a93-49db-a358-6b66f0bc3433
forces!(
f,
x0,
(i,j,x,y) -> fₓ(x,y,cutoff) # closure
)
# ╔═╡ 3755a4f3-1842-4de2-965e-d294c06c54c7
trajectory = md((
x0 = [random_vec(Vec2D{Float64},(-50,50)) for _ in 1:100 ],
v0 = [random_vec(Vec2D{Float64},(-1,1)) for _ in 1:100 ],
mass = [ 1.0 for _ in 1:100 ],
dt = 0.1,
nsteps = 1000,
isave = 10,
forces! = (f,x) -> forces!(f,x, (i,j,p1,p2) -> fₓ(p1,p2,cutoff))
)...)
# ╔═╡ 985b4ffb-7964-4b50-8c2f-e5f45f352500
trajectory_periodic = md((
x0 = [random_vec(Vec2D{Float64},(-50,50)) for _ in 1:100 ],
v0 = [random_vec(Vec2D{Float64},(-1,1)) for _ in 1:100 ],
mass = [ 10.0 for _ in 1:100 ],
dt = 0.1,
nsteps = 1000,
isave = 10,
forces! = (f,x) -> forces!(f,x,(i,j,p1,p2) -> fₓ(p1,p2,cutoff,side))
)...)
# ╔═╡ 1ad401b5-20b2-489b-b2aa-92f729b1d725
@benchmark md($(
x0 = [random_vec(Vec2D{Float64},-50:50) for _ in 1:100 ],
v0 = [random_vec(Vec2D{Float64},-1:1) for _ in 1:100 ],
mass = [ 1.0 for _ in 1:100 ],
dt = 0.1,
nsteps = 1000,
isave = 10,
forces! = (f,x) -> forces!(f,x, (i,j,p1,p2) -> fₓ(p1,p2,cutoff,side))
)...)
# ╔═╡ 0546ee2d-b62d-4c7a-8172-ba87b3c1aea4
trajectory_periodic_3D = md((
x0 = [random_vec(Vec3D{Float64},-50:50) for _ in 1:100 ],
v0 = [random_vec(Vec3D{Float64},-1:1) for _ in 1:100 ],
mass = [ 1.0 for _ in 1:100 ],
dt = 0.1,
nsteps = 1000,
isave = 10,
forces! = (f,x) -> forces!(f,x,(i,j,p1,p2) -> fₓ(p1,p2,cutoff,side))
)...)
# ╔═╡ 4e97f24c-c237-4117-bc57-e4e88c8fb8d2
md"""
Which generates random positions carrying an initial uncertainty we defined:
"""
# ╔═╡ b31da90d-7165-42de-b18d-90584affea03
random_vec(Vec2D{Measurement{Float64}},(-50,50),1e-5)
# ╔═╡ 5f37640b-ffd9-4877-a78c-a699b2671919
md"""
That given, the same simulation codes can be used to run the particles simulations while propagating the uncertinties of the coordinates of each position:
"""
# ╔═╡ 1d6eedfd-d013-4557-9cf2-103f8fb7b72a
md"""
The trajectory, of course, looks the same (except that we ran less steps, because propagating the error is expensive):
"""
# ╔═╡ c003a61d-a434-4d7b-9214-5b52aa044248
md"""
But now we have an estimate of the error of the positions, propagated from the initial uncertainty:
"""
# ╔═╡ 63eb391f-0238-434a-bc3a-2fa8ed41448e
md"""
### Planetary motion
Perhaps this is more interesting to see in a planetary trajectory:
"""
# ╔═╡ 7b9bb0fd-34a5-42e1-bc35-7259447b73d0
function gravitational_force(i,j,x,y,mass)
G = 0.00049823382528 # MKm³ / (10²⁴kg days²)
dr = y - x
r = norm(dr)
return G*mass[i]*mass[j]*dr/r^3
end
# ╔═╡ 6a4e0e2e-75c5-4cab-987d-3d6b62f9bb06
md"""
Note that now we need the indexes of the particles to be able to pass the information of their masses.
A set of planetary positions and velocities is something that we have to obtain [experimentaly](https://nssdc.gsfc.nasa.gov/planetary/factsheet/). Here, the distance units $10^6$ km, and time is in days. Thus, velocities are in MKm per day.
The uncertainty of the positions will be taken as the diameter of each planet. In this illustrative example we will not add uncertainties to the velcities.
"""
# ╔═╡ c91862dd-498a-4712-8e3d-b77e088cd470
planets_x0 = [
Vec2D( 0.0 ± 1.39 , 0. ± 1.39 ), # "Sun"
Vec2D( 57.9 ± 4.879e-3, 0. ± 4.879e-3), # "Mercury"
Vec2D(108.2 ± 12.104e-3, 0. ± 12.104e-3), # "Venus"
Vec2D(149.6 ± 12.756e-3, 0. ± 12.756e-3), # "Earth"
Vec2D(227.9 ± 6.792e-3, 0. ± 6.792e-3), # "Mars"
]
# ╔═╡ a08d6e6d-ddc4-40aa-b7c4-93ea03191415
planets_v0 = [
Vec2D(0. ± 0., 0.0 ± 0.), # "Sun"
Vec2D(0. ± 0., 4.10 ± 0.), # "Mercury"
Vec2D(0. ± 0., 3.02 ± 0.), # "Venus"
Vec2D(0. ± 0., 2.57 ± 0.), # "Earth"
Vec2D(0. ± 0., 2.08 ± 0.) # "Mars"
]
# ╔═╡ a356e2cc-1cb1-457a-986c-998cf1efe008
md"""
And the masses are given in units of $10^{24}$ kg:
"""
# ╔═╡ 57141f7c-9261-4dc5-98e4-b136a15f86fc
const masses = [ 1.99e6, 0.330, 4.87, 5.97, 0.642 ]
# ╔═╡ 055e32d7-073c-40db-a267-750636b9f786
md"""
Let us see the planets orbiting the sun:
"""
# ╔═╡ aaa97ce4-a5ff-4332-89a2-843cee2e5b6d
trajectory_planets = md((
x0 = planets_x0,
v0 = planets_v0,
mass = masses,
dt = 1, # days
nsteps = 2*365, # two Earth years
isave = 1, # save every day
forces! = (f,x) -> forces!(
f,x, (i,j,p1,p2) -> gravitational_force(i,j,p1,p2,masses)
)
)...)
# ╔═╡ 93697e4d-369b-48e9-8b28-a0ff58604d02
md"""
If you are wandering why the errors oscilate, it is because the trajectories are periodic. Whenever all possible trajectories starting from within the uncertainty interval cross each other, the error of the predicted position is independent on the initial coordinates. Thus, the derivative of the position is zero relative to initial position, and so it the propagated uncertainty when using a linear propagation rule.
"""
# ╔═╡ c4344e64-aa22-4328-a97a-71e44bcd289f
md"""
One thing I don't like, though, is that in two years the Earth did not complete two revolutions around the Sun. Something is wrong with our data. Can we improve that?
"""
# ╔═╡ 827bda6f-87d4-4d36-8d89-f144f4595240
md"""
## We can differentiate everything!
Perhaps astoningshly (at least for me), our simulation is completely differentiable. That means that we can tune the parameters of the simulation, and the data, using optimization algorithms that require derivatives.
Here we speculate that what was wrong with our data was that the initial position of the Earth was somewhat out of place. That caused the Earth orbit to be slower than it should.
We will define, then, an objective function which returns the displacement of the Earth relative to its initial position (at day one) after two years. Our goal is that after one year the Earth returns to its initial position.
"""
# ╔═╡ 1ff4077a-4742-4c5e-a8d6-c4699469a683
md"""
First, we define a function that executes a simulation of *two years* of an Earth orbit, starting from a given position for the Earth `x` coordinate as a parameter. We will be careful in making all other coordinates of the same type of `x`, so that the generality of the type of variable being used is kept consistent:
"""
# ╔═╡ 4a75498d-8f4e-406f-8b01-f6a5f153919f
function earth_orbit(x::T=149.6,nsteps=2*3650,isave=20) where T
x0 = [
Vec2D( zero(T), zero(T)), # "Sun"
Vec2D( x, zero(T)) # "Earth"
]
v0 = [
Vec2D( zero(T), zero(T)), # "Sun"
Vec2D( zero(T), 2.57*one(T)), # "Earth"
]
masses = [ 1.99e6, 5.97 ]
trajectory = md((
x0 = x0,
v0 = v0,
mass = masses,
dt = 0.1, # days
nsteps = nsteps, # one Earth year
isave = isave,
forces! = (f,x) -> forces!(f,x,
(i,j,p1,p2) -> gravitational_force(i,j,p1,p2,masses)
)
)...)
return trajectory
end
# ╔═╡ 3ae783ce-d06e-4cc2-b8a3-94512e8f1490
md"""
Now we define our objective function, consisting of the norm of the difference between the initial and final coordinates of the Earth after two year (what we want is that the Earth returns to its initial position):
"""
# ╔═╡ 13e7da81-8581-4f32-9fdb-2599dd36a12c
function error_in_orbit(x::T=149.6) where T
nsteps = 2*3650 # two years
traj = earth_orbit(x,nsteps,nsteps) # Save last point only
return norm(traj[end][2]-[x,0.])
end
# ╔═╡ 4870b1f3-3134-4ddc-a59d-fa806b456a23
md"""
We can see that our current data results in a significant error:
"""
# ╔═╡ fda6171c-9675-4f2e-b226-7ccf100529cd
error_in_orbit()
# ╔═╡ a862f8a3-0131-4644-bc90-246bf3120790
md"""
We want to minimize this error, and it turns out that your simulation is fully differentiable. We will use here the `ForwardDiff` automatic differentiation package:
"""
# ╔═╡ eee3ac4b-4ddb-4699-b6e6-f0ffcc562c07
md"""
Which can be used just as it it to compute the derivative of the error in the orbit relative to the initial `x` position of the Earth:
"""
# ╔═╡ 107aec28-ecb5-4007-95e5-25d0a7f0c465
ForwardDiff.derivative(error_in_orbit,149.6)
# ╔═╡ 1394e4c6-c371-47c0-8ca8-f0830d63d8ec
md"""
To minimize the error in the orbit we will write a simple stepest descent algorithm. Many packages are available for optimization, but here we will keep things simpler also to illustrate that writting the optimizer in Julia is a valid alternative:
"""
# ╔═╡ 535716e6-9c1c-4324-a4cd-b1214df3c01d
function gradient_descent(x,f,g,tol,maxtrial)
itrial = 0
step = 1.0
fx = f(x)
gx = g(x)
while (abs(gx) > tol) && (itrial < maxtrial) && (step > 1e-10)
xtrial = x - gx*step
ftrial = f(xtrial)
if ftrial > fx
step = step / 2
else
x = xtrial
fx = ftrial
gx = g(x)
step = step * 2
end
itrial += 1
end
return x, gx, itrial
end
# ╔═╡ b8edfb4e-6780-4ce7-94c1-4073ff7fa832
md"""
The derivative of our error can be computed by *closing over* the `error_in_orbit` function:
"""
# ╔═╡ b8320f78-323c-49a9-a9f9-2748d19ecb35
error_derivative(x) = ForwardDiff.derivative(error_in_orbit,x)
# ╔═╡ 92737d73-676c-4c96-a321-831ecaf37690
md"""
And now we can call the `gradient_descent` function directly:
"""
# ╔═╡ 931a9c5f-8f91-4e88-956b-50c0efc9c58b
best_x0 = gradient_descent(149.6,error_in_orbit,error_derivative,1e-4,1000)
# ╔═╡ b5b96082-efde-464f-bcd4-f2e0a84befcd
md"""
The result is reasonable: the error in the orbit has significantly being disminished:
"""
# ╔═╡ 7658a32c-d3da-4ec9-9d96-0d30bb18f08c
error_in_orbit(best_x0[1])
# ╔═╡ e61981d5-5448-45e9-81dc-320ac87ba813
md"""
Let us see our trajectory now with the new initial condition:
"""
# ╔═╡ 31e1bb51-c531-4c4a-8634-5caafb7e9e51
earth_traj_0 = earth_orbit(149.6)
# ╔═╡ b0b81da4-6788-45c4-b618-188a02b5e09c
earth_traj_best = earth_orbit(best_x0[1])
# ╔═╡ 47c205c3-ceae-4e12-9ade-753df1608deb
md"""
The dark blue dot is the corrected trajectory, and the light blue dot is the original one. Therefore, we were able to optimize the *initial point* of the trajectory with a gradient-based method. This concept can be used for adjusting parameters in simulations of many kinds (particle simulations or differential equations in general).
"""
# ╔═╡ 826693ff-9a9b-46b1-aeb3-767a5e6f9441
md"""
# Accelerating with CellListMap.jl
"""
# ╔═╡ d231842d-9b7a-4711-b71b-5d54041ebc1f
md"""
[`CellListMap.jl`](https://m3g.github.io/CellListMap.jl/stable/) is package aiming an efficient implementation of [cell lists](https://en.wikipedia.org/wiki/Cell_lists). Cell lists are practical algorithm to reduce the cost of computing short-ranged distances between particles. The package provides a general interface to compute any distance-dependent property, as potential energies and forces, nearest-neighbour lists, distribution functions, etc. It accepts systems with general (triclinic) periodic boundary conditions, in two and three dimensions.
The most simple cell list algorithm is relatively simple. Many optimizations can be done, however, on the construction of the lists, on the handling of periodic conditions, minimization of the number of unnecessary distance computations, and the parallelization of the construction of the lists and the mapping of the property to be evaluated.
"""
# ╔═╡ 53cedd26-3742-4c23-a8b8-8a1f2bdfa135
md"""
## The naive algorithm is too slow O(n²)
"""
# ╔═╡ 889f837d-2e26-4261-b276-5fd91efdda6a
md"""
With ~1k, particles, the number of pairs of particles is already of the order of hundreds of thousands. The naive O(n²) algorithm is already too slow. Typical simulations involve tenths of thousands to millions of particles.
"""
# ╔═╡ 670a01e3-82f8-4c7f-8577-852081d91ed7
md"""
Here, we will simulate 1000 particles to start:
"""
# ╔═╡ fce1e3e0-cdf7-453f-b913-964c10fa85a6
const n_large = 1000
# ╔═╡ 69a92ac6-833c-4605-b3d0-9400e4572886
md"""
Our previous system had 100 particles in a square of side 100. We will keep the density constant:
"""
# ╔═╡ 542a9ef5-d9ee-49bd-9d31-61e28b80b5cb
const box_side = sqrt(n_large / (100/100^2))
# ╔═╡ 8bada25c-b586-42b4-851d-232ccca8a456
md"""
We only need to generate the coordinates and run:
"""
# ╔═╡ 7600c6dc-769e-4c77-8526-281a1bcec079
x0_large = [ Vec2D(box_side*rand(),box_side*rand()) for _ in 1:n_large ]
# ╔═╡ 29dbc47b-3697-4fdf-8f34-890ab4d0cdae
t_naive = @elapsed trajectory_periodic_large = md((
x0 = x0_large,
v0 = [random_vec(Vec2D{Float64},(-1,1)) for _ in 1:n_large ],
mass = [ 10.0 for _ in 1:n_large ],
dt = 0.1,
nsteps = 1000,
isave = 10,
forces! = (f,x) -> forces!(f,x,(i,j,p1,p2) -> fₓ(p1,p2,cutoff,box_side))
)...)
# ╔═╡ 0ee7fc18-f41f-4179-a75e-1e1d56b2db29
md"""
Running time of naive algorithm: $t_naive seconds
"""
# ╔═╡ 0d0374ed-5150-40e6-b5a4-9a344b6ca47a
md"""
## Using cell lists
"""
# ╔═╡ f7cf613e-be9d-4f62-a778-cc4375eb99df
md"""
In cell lists, the particles are classified in cells before any distance computation. The distances are computed only for particles of vicinal cells. If the side of the cells is much smaller than the side of the complete system, the number of computations is drastically reduced.
"""
# ╔═╡ 5be87c6f-5c31-4d14-a8cb-4e63ef39d538
begin
function cell_list_picture()
function square(c,side)
x = [ c[1]-side/2, c[1]+side/2, c[1]+side/2, c[1]-side/2, c[1]-side/2]
y = [ c[2]-side/2, c[2]-side/2, c[2]+side/2, c[2]+side/2, c[2]-side/2]
return x, y
end
plt = plot()
x,y=square([5,5],2)
plot!(
plt,x,y,seriestype=[:shape],
linewidth=2,fillalpha=0.05,color="green",label=""
)
x,y=square([5,5],6)
plot!(
plt,x,y,seriestype=[:shape],
linewidth=2,fillalpha=0.05,color="orange",label=""
)
lines = collect(2:2:8)
vline!(plt,lines,color="gray",label="",style=:dash)
hline!(plt,lines,color="gray",label="",style=:dash)
px = [ 0.1 + 9.8*rand() for i in 1:100 ]
py = [ 0.1 + 9.8*rand() for i in 1:100 ]
scatter!(plt,px,py,label="",alpha=0.20,color="blue")
fontsize=8
annotate!(plt,3,3,text("(i-1,j-1)",fontsize,:Courier))
annotate!(plt,5,3,text("(i-1,j)",fontsize,:Courier))
annotate!(plt,7,3,text("(i-1,j+1)",fontsize,:Courier))
annotate!(plt,3,5,text("(i,j-1)",fontsize,:Courier))
annotate!(plt,5,5,text("(i,j)",fontsize,:Courier))
annotate!(plt,7,5,text("(i,j+1)",fontsize,:Courier))
annotate!(plt,3,7,text("(i+1,j-1)",fontsize,:Courier))
annotate!(plt,5,7,text("(i+1,j)",fontsize,:Courier))
annotate!(plt,7,7,text("(i+1,j+1)",fontsize,:Courier))
plot!(
plt,size=(400,400),
xlim=(1.3,8.7),xticks=:none,
ylim=(1.3,8.7),yticks=:none,
framestyle=:box,
xlabel="x",ylabel="y",grid=false
)
return plt
end
cell_list_picture()
end
# ╔═╡ 0c07edd3-c0a1-4f72-a16a-74badb7a6123
md"""
To use `CellListMap.jl` we need to setup our system, by providing the data on the box properties and the cutoff of the interactions:
"""
# ╔═╡ 4fc5ef4d-e072-41f7-aef9-b42730c8313c
box = Box([box_side,box_side],cutoff)
# ╔═╡ 19c5cc9d-8304-4e36-a3ea-a1151f28f71d
md"""
The particles are then classified in the cells. Virtual (ghost) particles are created at the boundaries to handle peridic boundary conditions and avoid having to wrap coordinates during the pairwise computation stage:
"""
# ╔═╡ 7dcadd85-2986-4e42-aa84-67128a8f666d
cl = CellList(x0_large,box)
# ╔═╡ 0b5c6ede-bceb-499a-a9a8-3c6a75ed340a
md"""
We only need to implement the function that has to be evaluated *if$ the particles are closer than the cutoff. This function will only be called in that case. Here, the function will update the force vector:
"""
# ╔═╡ 91b5eac1-4799-4a72-ac6a-e2b117b787d5
function fpair_cl(x,y,i,j,d2,f,box::Box)
Δv = y - x
d = sqrt(d2)
fₓ = 2*(d - box.cutoff)*(Δv/d)
f[i] += fₓ
f[j] -= fₓ
return f
end
# ╔═╡ 0f86ab3c-29aa-472b-8194-228c736ee940
md"""
The function that computes the forces in our simulation will, then, consist of an update of the cell lists followed by a call to the `map_pairwise!` function of `CellListMap.jl`, which takes as arguments the function to be mapped (`fpair_cl` here), the initial value of the forces vector `f`, and the system properties. We run only the serial version in this example:
"""
# ╔═╡ 0b8a2292-c0d6-44e4-b560-32d9d579a008
function forces_cl!(f::Vector{T},x,box::Box,cl::CellList,fpair::F) where {T,F}
fill!(f,zero(T))
cl = UpdateCellList!(x,box,cl,parallel=false)
map_pairwise!(
(x,y,i,j,d2,f) -> fpair(x,y,i,j,d2,f,box),
f, box, cl, parallel=false
)
return f
end
# ╔═╡ d6585cca-78bf-41d1-aea3-01d9831d76cb
md"""
With a proper definition of the function to compute forces, we can now run again the simulation:
"""
# ╔═╡ 1b7b7d48-79d2-4317-9045-5b7e7bd073e5
t_cell_lists = @elapsed trajectory_cell_lists = md((
x0 = x0_large,
v0 = [random_vec(Vec2D{Float64},(-1,1)) for _ in 1:n_large ],
mass = [ 10.0 for _ in 1:n_large ],
dt = 0.1,
nsteps = 1000,
isave = 10,
forces! = (f,x) -> forces_cl!(f,x,box,cl,fpair_cl)
)...)
# ╔═╡ 3f9dad58-294c-405c-bfc4-67855bb1e825
md"""
Running time of CellListMap: $t_cell_lists seconds (on the second run - compilation takes about 2 seconds).
"""
# ╔═╡ 6d61b58f-b88f-48f4-8bdd-0bb1a8bc1c82
md"""
Even for a small system like this one, the speedup is significant (of about $(round(Int,t_naive/t_cell_lists)) times here).
"""
# ╔═╡ 76b8695e-64dc-44bc-8938-ce22c4a9e4d0
md"""
## Energy minimization: the Packmol strategy
"""
# ╔═╡ 372637ff-9305-4d45-bf6e-e6531dadbd14
md"""
### A Lennard-Jones potential energy
Molecular dynamics simulations usually involve computing, for each pair of atoms, a Lennard-Jones function of the form: