-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmonteswitch_mod.f95
More file actions
6761 lines (6463 loc) · 284 KB
/
Copy pathmonteswitch_mod.f95
File metadata and controls
6761 lines (6463 loc) · 284 KB
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
!-----------------------------------------------------------------------
!
! Copyright (c) 2016 Tom L. Underwood
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject to
! the following conditions:
!
! The above copyright notice and this permission notice shall be
! included in all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
! BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
! ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
! CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
!-----------------------------------------------------------------------
!
! By extracting lines beginning with the regular expression '\!\! ?'
! (ignoring leading whitespace), and then removing matches to the
! regular expression, html documentation corresponding to this
! source code will be created.
!
! By extracting lines beginning with the regular expression '\!init_params ?'
! (ignoring leading whitespace), and then removing matches to the
! regular expression, a template of an 'init_params' file will be created,
! i.e. a file corresponding to the 'filename_params' argument in the
! subroutine 'initialise_from_files(filename_params,filename_lattice)'.
! However, this template does not include the contributions from the
! 'energy_mod' module, which must be added to the end of the template.
!
!! <html>
!! <head>
!! <title> monteswitch_mod Documentation </title>
!! </head>
!! <body>
!!
!! <h1> <code>monteswitch_mod</code> Documentation </h1>
!!
!! <h2> Author </h2>
!! <p> Tom Underwood </p>
!!
!! <h2> General Notes </h2>
!!
!! <h3> Description </h3>
!! <p>
!! This module contains the nuts and bolts for lattice switch Monte Carlo simulations. The variables in this module define the
!! 'state' of the simulation, which includes - among other things - variables describing the microstate of the
!! system at the current time and variables determining in what manner information will be output during a simulation.
!! The variables associated with particle interactions are inherited from the module <code>interactions_mod</code>.
!! </p>
!! <p>
!! Particles are moved according to a random walk, with a maximum move of <code>part_step</code> in any Cartesian
!! direction. The volume is moved, if volume moves are used at all, in a manner according the flag
!! <code>vol_dynamics</code>. Details of the various options are provided below.
!! </p>
!! <p>
!! This module terminates the current program with a non-zero exit status of 1 for most errors. If the system has
!! melted, and <code>melt_option="stop"</code> then the exit status is 2. The exit status is also 2 if the system has 'exploded',
!! i.e., it has negative volume. If the energy has diverged due to precision-related
!! errors then the exit status is 3. Note that exit statuses are not part of the Fortran standard, and may not work for
!! all operating systems or compilers.
!! </p>
!!
!! <h3> Definitions and conventions </h3>
!! <p>
!! <ul>
!! <li>
!! The free energy difference is the free energy of phase 1 minus that of phase 2.
!! </li>
!! <li>
!! A "move" is one of either a: particle move, in which one particle is moved; a lattice move, in which
!! the lattice is switched; and a volume move, in which the volume of the unit cell is altered. Note that
!! if the move is rejected then it is still deemed to have taken place. For an NVT ensemble with lattice
!! switch dynamics, the following move "cycle" is performed: particle move, lattice move. For an NPT
!! ensemble with lattice switch dynamics the move cycle is: particle move, lattice move, (volume move,
!! lattice move), where the set of moves in the brackets occur on average <code>vol_freq</code> times per
!! sweep. A different particle is used for each cycle; the particles are either cycled through one by one,
!! or selected randomly for each cycle, according to the flag <code>part_select</code>.
!! If the flag <code>enable_lattice_moves</code> is set to <code>.false</code>
!! then lattice switch dynamics are not used, and lattice moves are never attempted.
!! </li>
!! <li>
!! A "sweep" comprises a cycle containing <code>n_part</code> particle moves, where <code>n_part</code> is the number
!! of particles in the system. When the number of sweeps is printed to stdout
!! or a file, specifically it is the number of <i>completed</i> sweeps. Tasks performed periodically on the scale of
!! sweeps, e.g. updating the weight function, outputing data to files, checkpointing, are performed at the <i>end</i> of the
!! sweep. The following table gives a breakdown of the number of moves of each type per sweep for different
!! combinations of flags.
!! <table border="1">
!! <tr>
!! <td> <b> Dynamics </b> </td>
!! <td> <b> Particle moves </b> </td>
!! <td> <b> Volume moves </b> </td>
!! <td> <b> Lattice moves </b> </td>
!! <td> <b> Total moves </b> </td>
!! </tr>
!! <tr>
!! <td>
!! <code>enable_part_moves=.true.</code>
!! <code>enable_vol_moves=.true.</code>
!! <code>enable_lattice_moves=.true.</code>
!! </td>
!! <td> <code>n_part</code> </td>
!! <td> <code>v_freq</code> (on average) </td>
!! <td> <code>(n_part+vol_freq)</code> (on average)</td>
!! <td> <code>2*(n_part+vol_freq)</code> (on average)</td>
!! </tr>
!! <tr>
!! <td>
!! <code>enable_part_moves=.true.</code>
!! <code>enable_vol_moves=.false.</code>
!! <code>enable_lattice_moves=.true.</code>
!! </td>
!! <td> <code>n_part</code> </td>
!! <td> 0 </td>
!! <td> <code>n_part</code> </td>
!! <td> <code>2*n_part</code> </td>
!! </tr>
!! </tr>
!! <tr>
!! <td>
!! <code>enable_part_moves=.true.</code>
!! <code>enable_vol_moves=.true.</code>
!! <code>enable_lattice_moves=.false.</code>
!! </td>
!! <td> <code>n_part</code> </td>
!! <td> <code>v_freq</code> (on average) </td>
!! <td> 0 (on average)</td>
!! <td> <code>(n_part+vol_freq)</code> (on average)</td>
!! </tr>
!! </tr>
!! <tr>
!! <td>
!! <code>enable_part_moves=.true.</code>
!! <code>enable_vol_moves=.false.</code>
!! <code>enable_lattice_moves=.false.</code>
!! </td>
!! <td> <code>n_part</code> </td>
!! <td> 0 </td>
!! <td> 0 </td>
!! <td> <code>n_part</code></td>
!! </tr>
!! <tr>
!! <td>
!! <code>enable_part_moves=.false.</code>
!! <code>enable_vol_moves=.false.</code>
!! <code>enable_lattice_moves=.true.</code>
!! </td>
!! <td> 0 </td>
!! <td> 0 </td>
!! <td> <code>n_part</code></td>
!! <td> <code>n_part</code></td>
!! </tr>
!! <tr>
!! <td>
!! <code>enable_part_moves=.false.</code>
!! <code>enable_vol_moves=.true.</code>
!! <code>enable_lattice_moves=.true.</code>
!! </td>
!! <td> 0 </td>
!! <td> <code>v_freq</code> (on average) </td>
!! <td> <code>(n_part+vol_freq)</code> (on average)</td>
!! <td> <code>(n_part+vol_freq)</code> (on average)</td>
!! </tr>
!! </table>
!! </li>
!! <li>
!! The order parameter for a microstate is defined as <code>E_1-E_2+P*(V_1-V_2)-log(V_1/V_2)/beta</code>,
!! where <code>E_1</code> is the energy of the lattice-type-1 microstate corresponding to the current displacements,
!! <code>E_2</code> is the energy of the lattice-type-2 microstate, <code>V_1</code> is the volume of the lattice-type-1
!! microstate, <code>V_2</code> is the volume of the lattice-type-2 microstate, and <code>P</code> and <code>beta</code>
!! have their usual significance.
!! </li>
!! <li>
!! The weight function is such that the probability of microstate 'sigma' being sampled using multicanonical
!! sampling is exp(-beta*E(sigma))exp(eta(sigma)), where eta(sigma) is the value of the weight fuunction for
!! sigma.
!! </li>
!! </ul>
!! </p>
!!
!! <h3> Warnings </h3>
!! <p>
!! <ul>
!! <li>
!! The program may crash without any error messages if the Fortran intrinsic function <code>mod(a,p)</code> is
!! called with <code>p=0</code>. This could occur if any of the variables which describe a 'period' in
!! sweeps is 0. Therefore do not set any of these to 0 if they are in use. Depending on how I've coded, it
!! may even be the case that the program crashes if any of these are 0 and they aren't 'in use'. I will have
!! to test this at some point; to be safe just don't use 0 for such variables.
!! </li>
!! <li>
!! There are no default values for any of the variables. It is therefore recommended that their values are all
!! explicitly stated at initialisation: one would definitely not want a strange value of, say,
!! <code>stop_sweeps</code> creeping into the simulation because its value wasn't initialised.
!! </li>
!! <li>
!! In a similar vein to the above point, if one, say, tells the simulation to not evaluate equilibrium
!! quantities, then the values of the variables associated with equilibrium quantities will be nonsense, and
!! do not imply a bug in the program.
!! </li>
!! </ul>
!! </p>
!!
module monteswitch_mod
!! <h2> Dependencies </h2>
!! <p>
!! <ul>
!! <li> <code> kinds_mod </code> </li>
!! <li> <code> rng_mod </code> </li>
!! <li> <code> metropolis_mod </code> </li>
!! <li> <code> interactions_mod </code> </li>
!! </ul>
!! </p>
use kinds_mod
use rng_mod
use metropolis_mod
use interactions_mod
implicit none
!! <h2> Variables </h2>
!! <p>
!! What follows is a list of the module variables. Note also that variables are inherited from <code>interactions_mod</code>
!! which contribute to the specification of the state of the system, such as neighbour lists.
!! Variables can be loosely categorised as either user-defined or 'internal'. In theory a user need never tamper with
!! internal variables - and I would not recommend it unless they <i>really</i> know what they are doing. Internal
!! variables' descriptions are in red.
!! </p>
!!
!! <h3> Variables determining the ensemble and dynamics</h3>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! <td> <b> Description </b> </td>
!! </tr>
!! <tr>
!! <td> <code>enable_multicanonical<code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag which determines whether multicanonical sampling is used. If set to <code>.true.</code> then
!! multicanonical sampling is used via the weight function <code>eta</code>; if set to <code>.false</code>
!! then Boltzmann sampling is used.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>beta<code> </td>
!! <td> <code>real(rk)</code> </td>
!! <td>
!! The thermodynamic beta. This should be greater than zero - setting this to zero will yield division by
!! zero errors in some quantities, e.g. the free energy difference between phases
!! </td>
!! </tr>
!! <tr>
!! <td> <code>P<code> </td>
!! <td> <code>real(rk)</code> </td>
!! <td> Pressure of the system. </td>
!! </tr>
!! <tr>
!! <td> <code>enable_lattice_moves<code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag which which determines whether lattice switch dynamics is used. If set to <code>.true.</code>,
!! then after every particle or volume move a lattice switch is attempted. If set to
!! <code>.false.</code> then conventional dynamics are used: the underlying lattice is never changed
!! during the course of a simulation.
!! </td>
!! <tr>
!! <td> <code>enable_part_moves<code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag which which determines whether particles will be moved: if set to <code>.true.</code>,
!! then particles will be moved. This exists mainly for debugging purposes.
!! </td>
!! <tr>
!! <td> <code>enable_vol_moves<code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag which which determines whether volume moves are used, and determines the ensemble to be used:
!! if set to <code>.true.</code> then we have an NPT ensemble; if set to <code>.false.</code> then we
!! have an NVT ensemble. Note that the Metropolis algorithm differs slightly between NVT and NPT
!! ensembles.
!! </td>
!! <tr>
!! <td> <code>part_select<code> </td>
!! <td> <code>character(len=30)</code> </td>
!! <td>
!! Flag which which determines how a particle is to be selected at the
!! beginning of each cycle of moves. Possible values are as follows:
!! <code>"cycle"</code> results in selecting particles in order, i.e. particle
!! 1 is used first, then particle 2,..., then particle <code>n_part</code>, then
!! particle 1, etc.; <code>"rand"</code> selects a particle at the beginning of each
!! cycle at random.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>part_step<code> </td>
!! <td> <code>real(rk)</code> </td>
!! <td>
!! Maximum displacement change in any Cartesian direction which any particle can be
!! moved by during a particle move. Recall that random walk particle dynamics is used.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>enable_COM_frame<code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag which determines whether we use the centre-of-mass frame. If set to <code>.true.</code>
!! the particle displacements are ammended after every accepted particle move such that the
!! average displacement is 0.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>vol_dynamics<code> </td>
!! <td> <code>character(len=30)</code> </td>
!! <td>
!! Flag which determines how the supercell is altered for a volume move. Note
!! that this is only important for NPT ensembles. Possible values are as follows:
!! <code>"FVM"</code> uses a random walk in ln(V) with a maximum step size of <code>vol_step<code>,
!! and a fixed aspect ratio, i.e. the proportions of the supercell dimensions are unchanged by a
!! volume move; <code>"UVM"</code> uses a random walk to expand/contract each dimension of the supercell.
!! </td>
!! <tr>
!! <td> <code>vol_freq<code> </td>
!! <td> <code>integer(ik)</code> </td>
!! <td>
!! The number of volume moves which will occur per sweep on average.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>vol_step<code> </td>
!! <td> <code>real(rk)</code> </td>
!! <td>
!! Maximum step size for a volume move random walk. The volume can change by a maximum of <code>vol_step<code>
!! each volume move.
!! </td>
!! </tr>
!! </table>
logical :: enable_multicanonical
real(rk) :: beta
real(rk) :: P
logical :: enable_lattice_moves
logical :: enable_part_moves
logical :: enable_vol_moves
character(len=30) :: part_select
character(len=*), parameter :: part_select_cycle="cycle"
character(len=*), parameter :: part_select_rand="rand"
real(rk) :: part_step
logical :: enable_COM_frame
character(len=30) :: vol_dynamics
character(len=*), parameter :: vol_dynamics_FVM="FVM"
character(len=*), parameter :: vol_dynamics_UVM="UVM"
integer(ik) :: vol_freq
real(rk) :: vol_step
!! <h3> Variables determining the length of the simulation and the equilibration time </h3>
!! <p>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! <td> <b> Description </b> </td>
!! </tr>
!! <tr>
!! <td> <code>stop_sweeps</code> </td>
!! <td> <code>integer(ik)</code> </td>
!! <td>
!! The number of sweeps to perform during this simulation. If this is set to 0 then the Monte Carlo loop
!! is by-passed, but tasks performed periodically during the Monte Carlo loop, e.g. updating the weight
!! function and checking whether or not the system has melted/exploded, are done once. Furthermore, if it is set to
!! 0, the equilibrium properties - the free energy difference between phases, the volume and enthalpy of
!! each phase, and the associated uncertainties of these quantities - are recalculated if they are to be
!! calculated.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>equil_sweeps</code> </td>
!! <td> <code>integer(ik)</code> </td>
!! <td> The number of sweeps after which equilibration is assumed to have taken place. </td>
!! </tr>
!! </table>
integer(ik) :: stop_sweeps
integer(ik) :: equil_sweeps
!! <h3> Variables pertaining to safety checks on the simulation </h3>
!! <p>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! <td> <b> Description </b> </td>
!! </tr>
!! <tr>
!! <td> <code>enable_melt_checks</code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag determining whether the simulation will periodically check if the particle displacements have
!! exceeded 'safe' values, i.e. if the system has melted, and if the system has a negative volume, i.e.,
!! the system has 'exploded'.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>melt_sweeps</code> </td>
!! <td> <code>integer(ik)</code> </td>
!! <td>
!! If <code>enable_melt_checks=.true.</code>, the particle displacements are checked for melting/exploding every
!! <code>melt_sweeps</code> sweeps.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>melt_threshold</code> </td>
!! <td> <code>real(rk)</code> </td>
!! <td>
!! If any particle displacement is of magnitude <code>>melt_threshold</code> in any Cartesian direction then the system
!! is deemed to have melted.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>melt_option</code> </td>
!! <td> <code>character(len=30)</code> </td>
!! <td>
!! Flag which determines what happens if the system melts. Options are as follows. <code>"zero_1"</code> and
!! <code>"zero_2"</code> move the system to the zero displacement microstate at the current volume in lattice types 1
!! and 2 respectively, and then proceeds. <code>"zero_random"</code> does the same but picks a random lattice type, and
!! <code>"zero_current"</code> picks the current lattice.
!! In all of these cases, the system waits for <code>equil_sweeps</code> after the displacements are 'reset' for the
!! system to equilibrate before variables which should not be updated before equilibration are updated. Also, the current
!! block with regards to the calculation of equilibrium quantities is disregarded.
!! <code>"stop"</code> prints an error message to stderr if the system melts and returns an exit status of 2.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>enable_divergence_checks</code> </td>
!! <td> <code>logical</code> </td>
!! <td>
!! Flag determining whether the simulation will periodically check if the energy is correct. For particle moves, the energy
!! is not calculated exactly for the trial particle displacements because it is very computationally expensive and unnecessary.
!! Instead the energy <i>change</i> for the trial displacements is calculated, which involves consideration of the energy
!! stored in the interactions between the particle which has been moved, and its neighbours. This is far less demanding to
!! calculate. If the move is accepted this change is ammended to the total energy. However, over time it is possible that this
!! 'running total' approach will yield incorrect energies due to the finite precision of the computer. Hence one should periodically
!! recalculate the energy exactly. This is done during volume moves. If <code>enable_divergence_checks=.true.</code> this is also
!! done every <code>divergence_sweeps</code> sweeps, after which the recalculated energy is compared to the 'current' energy, and
!! an error is flagged if they are different - outwith a tolerence of <code>divergence_tol</code>. Note that the order parameter
!! is also ammended after the energy (for lattice 1 and 2) is recalculated.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>divergence_sweeps</code> </td>
!! <td> <code>integer(ik)</code> </td>
!! <td>
!! If <code>enable_divergence_checks=.true.</code>, the energy is checked every <code>divergence_sweeps</code> sweeps.
!! </td>
!! </tr>
!! <tr>
!! <td> <code>divergence_tol</code> </td>
!! <td> <code>real(rk)</code> </td>
!! <td>
!! The difference between the current energy and the 'true' energy of the system which would count as an error. This should
!! be positive.
!! </td>
!! </tr>
!! </table>
logical :: enable_melt_checks
integer(ik) :: melt_sweeps
real(rk) :: melt_threshold
character(len=30) :: melt_option
character(len=30), parameter :: melt_option_zero_1="zero_1"
character(len=30), parameter :: melt_option_zero_2="zero_2"
character(len=30), parameter :: melt_option_zero_random="zero_random"
character(len=30), parameter :: melt_option_zero_current="zero_current"
character(len=30), parameter :: melt_option_stop="stop"
logical :: enable_divergence_checks
integer(ik) :: divergence_sweeps
real(rk) :: divergence_tol
!! <h3> Variables determining the nature of the output </h3>
!! <p>
!! The flags <code>output_file_X</code>, where <code>X</code> is the name of a simulation variable, when set to <code>.true</code>,
!! will result in a line consisting of "<code>X: </code>" followed by the number of completed sweeps, followed by the current value of
!! <code>X</code> being printed to a specified file (<code>datafile</code> in the <code>run</code> subroutine) every
!! <code>output_file_period</code> sweeps (at the beginning of the sweep) - where the file is an argument of the <code>run</code>
!! subroutine. However, if <code>output_file_period=0</code>, then instead the output is <i>after every move</i>; and if
!! <code>output_file_period<0</code> then there is no output to the file. Similar applies to <code>output_stdout_X</code>, but for
!! stdout during a simulation.
!! </p>
!! <p>
!! The variable <code>checkpoint_period</code> determines how often (in sweeps) the simulation is checkpointed, i.e., how often the
!! state of the simulation is exported to a specified file (<code>statefile</code> in the <code>run</code> subroutine). Note that
!! the simulation is automatically checkpointed at the end of every simulation. If <code>checkpoint_period<=0</code>, then there is
!! no output to the file <code>statefile</code>, i.e., it is empty. If one wants the simulation to be checkpointed only at the end
!! of a simulation, and not during it, then set <code>checkpoint_period>stop_sweeps</code>.
!! </p>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! </tr>
!! <tr> <td> <code>output_file_period</code> </td> <td> <code>integer(ik)</code> </td> </tr>
!! <tr> <td> <code>output_file_Lx </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_Ly </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_Lz </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_V </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_R_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_R_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_u </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_lattice </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_E </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_M </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_macro </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_eta </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_moves_lattice </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_accepted_moves_lattice </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_moves_part </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_accepted_moves_part </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_moves_vol </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_accepted_moves_vol </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_rejected_moves_M_OOB </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_M_OOB_high </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_barrier_macro_low </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_barrier_macro_high </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_rejected_moves_M_barrier </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_moves_since_lock </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_melts </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_DeltaFs </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_DeltaFs </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_H_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_H_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_H_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_H_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_V_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_V_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_V_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_V_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_umsd_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_umsd_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_umsd_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_umsd_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_L_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_L_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_equil_L_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_file_sigma_equil_L_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_period</code> </td> <td> <code>integer(ik)</code> </td> </tr>
!! <tr> <td> <code>output_stdout_Lx </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_Ly </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_Lz </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_V </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_R_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_R_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_u </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_lattice </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_E </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_M </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_macro </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_eta </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_accepted_moves_lattice </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_accepted_moves_part </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_accepted_moves_vol </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_rejected_moves_M_OOB </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_M_OOB_high </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_barrier_macro_low </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_barrier_macro_high </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_rejected_moves_M_barrier </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_moves_since_lock </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_DeltaFs </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_DeltaFs </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_H_1s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_H_1s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_H_2s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_H_2s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_V_1s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_V_1s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_V_2s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_V_2s </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_melts </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_umsd_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_umsd_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_umsd_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_umsd_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_L_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_L_1 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_equil_L_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>output_stdout_sigma_equil_L_2 </code> </td> <td> <code>logical</code> </td> </tr>
!! <tr> <td> <code>> <code>checkpoint_period </code> </td> <td> <code>integer(ik)</code> </td> </tr> </code> </td> <td> <code>logical</code> </td> </tr>
!! </table>
integer(ik) :: output_file_period
logical :: output_file_Lx
logical :: output_file_Ly
logical :: output_file_Lz
logical :: output_file_V
logical :: output_file_R_1
logical :: output_file_R_2
logical :: output_file_u
logical :: output_file_lattice
logical :: output_file_E
logical :: output_file_M
logical :: output_file_macro
logical :: output_file_eta
logical :: output_file_moves_lattice
logical :: output_file_accepted_moves_lattice
logical :: output_file_moves_part
logical :: output_file_accepted_moves_part
logical :: output_file_moves_vol
logical :: output_file_accepted_moves_vol
logical :: output_file_rejected_moves_M_OOB
logical :: output_file_M_OOB_high
logical :: output_file_M_OOB_low
logical :: output_file_barrier_macro_low
logical :: output_file_barrier_macro_high
logical :: output_file_rejected_moves_M_barrier
logical :: output_file_moves_since_lock
logical :: output_file_melts
logical :: output_file_equil_DeltaF
logical :: output_file_sigma_equil_DeltaF
logical :: output_file_equil_H_1
logical :: output_file_sigma_equil_H_1
logical :: output_file_equil_H_2
logical :: output_file_sigma_equil_H_2
logical :: output_file_equil_V_1
logical :: output_file_sigma_equil_V_1
logical :: output_file_equil_V_2
logical :: output_file_sigma_equil_V_2
logical :: output_file_equil_umsd_1
logical :: output_file_sigma_equil_umsd_1
logical :: output_file_equil_umsd_2
logical :: output_file_sigma_equil_umsd_2
logical :: output_file_equil_L_1
logical :: output_file_sigma_equil_L_1
logical :: output_file_equil_L_2
logical :: output_file_sigma_equil_L_2
integer(ik) :: output_stdout_period
logical :: output_stdout_Lx
logical :: output_stdout_Ly
logical :: output_stdout_Lz
logical :: output_stdout_V
logical :: output_stdout_R_1
logical :: output_stdout_R_2
logical :: output_stdout_u
logical :: output_stdout_lattice
logical :: output_stdout_E
logical :: output_stdout_M
logical :: output_stdout_macro
logical :: output_stdout_eta
logical :: output_stdout_moves_lattice
logical :: output_stdout_accepted_moves_lattice
logical :: output_stdout_moves_part
logical :: output_stdout_accepted_moves_part
logical :: output_stdout_moves_vol
logical :: output_stdout_accepted_moves_vol
logical :: output_stdout_rejected_moves_M_OOB
logical :: output_stdout_M_OOB_high
logical :: output_stdout_M_OOB_low
logical :: output_stdout_barrier_macro_low
logical :: output_stdout_barrier_macro_high
logical :: output_stdout_rejected_moves_M_barrier
logical :: output_stdout_moves_since_lock
logical :: output_stdout_melts
logical :: output_stdout_equil_DeltaF
logical :: output_stdout_sigma_equil_DeltaF
logical :: output_stdout_equil_H_1
logical :: output_stdout_sigma_equil_H_1
logical :: output_stdout_equil_H_2
logical :: output_stdout_sigma_equil_H_2
logical :: output_stdout_equil_V_1
logical :: output_stdout_sigma_equil_V_1
logical :: output_stdout_equil_V_2
logical :: output_stdout_sigma_equil_V_2
logical :: output_stdout_equil_umsd_1
logical :: output_stdout_sigma_equil_umsd_1
logical :: output_stdout_equil_umsd_2
logical :: output_stdout_sigma_equil_umsd_2
logical :: output_stdout_equil_L_1
logical :: output_stdout_sigma_equil_L_1
logical :: output_stdout_equil_L_2
logical :: output_stdout_sigma_equil_L_2
integer(ik) :: checkpoint_period
!! <h3> Variables defining macrostates </h3>
!! <p>
!! In a simulation we need to define macrostates, each of which consists of a range of order parameters.
!! <code>M_grid</code> is an array containing <code>M_grid_size</code> equidistant
!! values. We emphasise that <code>M_grid</code> must have this form.
!! Furthermore, <code>M_grid_size</code> must be greater than 1.
!! <code>M_grid</code> defines macrostates: the 'i'th macrostate is defined as consisting of microstates
!! whose order parameter is between <code>M_grid(i)</code> (inclusive) and <code>M_grid(i)+(M_grid(2)-M_grid(1))</code>
!! (exclusive) - with <code>(M_grid(2)-M_grid(1))</code> being the spacing between successive elements of <code>M_grid</code>. Moves
!! which take us to states with <code>M>=M_grid(M_grid_size)+(M_grid(2)-M_grid(1))</code> or <code>M<M_grid(1)</code>
!! are always rejected.
!! </p>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! </tr>
!! <tr>
!! <td> <code>M_grid_size</code> </td>
!! <td> <code>integer(ik)</code> </td>
!! </tr>
!! <tr>
!! <td> <code>M_grid</code> </td>
!! <td> <code> real(rk), dimension(:), allocatable</code> </td>
!! </tr>
!! </table>
integer(ik) :: M_grid_size
real(rk), dimension(:), allocatable :: M_grid
!! <h3> Variables describing the state of the system </h3>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! <td> <b> Description </b> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>n_part</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red"> Total number of particles in the system. </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>Lx</code> </font> </td>
!! <td> <font color="red"> <code> real(rk), dimension(2)</code> </font> </td>
!! <td> <font color="red"> Dimension of supercell in x direction for lattice types 1 and 2. </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>Ly</code> </font> </td>
!! <td> <font color="red"> <code> real(rk), dimension(2)</code> </font> </td>
!! <td> <font color="red"> Dimension of supercell in y direction for lattice types 1 and 2. </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>Lz</code> </font> </td>
!! <td> <font color="red"> <code> real(rk), dimension(2)</code> </font> </td>
!! <td> <font color="red"> Dimension of supercell in z direction for lattice types 1 and 2. </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>V</code> </font> </td>
!! <td> <font color="red"> <code> real(rk)</code> </font> </td>
!! <td> <font color="red"> Current volume of the system (for NPT simulations) </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>spec_1</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik), dimension(:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! <code>spec_1(n)</code> is the species of the nth particle for lattice type 1.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>spec_2</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik), dimension(:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! <code>spec_2(n)</code> is the species of the nth particle for lattice type 2.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>pos_1</code> </font> </td>
!! <td> <font color="red"> <code>real(rk), dimension(:,:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! (<code>pos_1(1,n)</code>,<code>pos_1(2,n)</code>,<code>pos_1(3,n)</code>) is the position
!! of the nth particle for lattice type 1.
!! This array should have <code>n_part</code> elements in its first dimension,
!! and 3 elements in its second.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>pos_2</code> </font> </td>
!! <td> <font color="red"> <code>real(rk), dimension(:,:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! (<code>pos_2(1,n)</code>,<code>pos_2(2,n)</code>,<code>pos_2(3,n)</code>) is the position
!! of the nth particle for lattice type 2.
!! This array should have <code>n_part</code> elements in its first dimension,
!! and 3 elements in its second.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>R_1</code> </font> </td>
!! <td> <font color="red"> <code>real(rk), dimension(:,:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! (<code>R_1(1,n)</code>,<code>R_1(2,n)</code>,<code>R_1(3,n)</code>) is the position
!! of the nth particle with zero displacement for lattice type 1.
!! This array should have <code>n_part</code> elements in its first dimension,
!! and 3 elements in its second.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>R_2</code> </font> </td>
!! <td> <font color="red"> <code>real(rk), dimension(:,:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! (<code>R_2(1,n)</code>,<code>R_2(2,n)</code>,<code>R_2(3,n)</code>) is the position
!! of the nth particle with zero displacement for lattice type 2.
!! This array should have <code>n_part</code> elements in its first dimension,
!! and 3 elements in its second.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>u</code> </font> </td>
!! <td> <font color="red"> <code>real(rk), dimension(:,:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! (<code>u(1,n)</code>,<code>u(2,n)</code>,<code>u(3,n)</code>) is the current displacement vector
!! for the nth particle in the system, i.e., if the current lattice
!! is lattice type 1, then <code>u(:,n)+R_1(:,n)</code> are the current positions
!! of the particles in the system. This array should have <code>n_part</code>
!! elements in its first dimension, and 3 elements in its second.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>lattice</code> </font> </td>
!! <td> <font color="red"> <code>integer</code> </font> </td>
!! <td> <font color="red">
!! The current lattice type (1 or 2).
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>E_1</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The energy associated with lattice type 1 for the current set of
!! displacements and volume.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>E_2</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The energy associated with lattice type 2 for the current set of
!! displacements and volume.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>E</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The current energy. (<code>E=E_1</code> if <code>lattice=1</code> and <code>E=E_2</code> if <code>lattice=2</code>).
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>M</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The current order parameter.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>macro</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">
!! The current macrostate.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>eta</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The value of the weight function for the current state.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>M_const</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The value of <code>log(V_1/V_2)/beta</code> determined at the beginning of the simulation, where
!! <code>V_1</code> is the volume of box 1. Note that this is a constant throughout a simulation since <code>V_1/V_2</code>
!! is preserved by the lattice switch. This quantity is stored in memory for optimisation purposes -
!! to avoid calling the log function every move.
!! </font> </td>
!! </tr>
!! </table>
integer(ik) :: n_part
real(rk), dimension(2) :: Lx
real(rk), dimension(2) :: Ly
real(rk), dimension(2) :: Lz
real(rk) :: V
integer(ik), dimension(:), allocatable :: spec_1
integer(ik), dimension(:), allocatable :: spec_2
real(rk), dimension(:,:), allocatable :: pos_1
real(rk), dimension(:,:), allocatable :: pos_2
real(rk), dimension(:,:), allocatable :: R_1
real(rk), dimension(:,:), allocatable :: R_2
real(rk), dimension(:,:), allocatable :: u
integer :: lattice
real(rk) :: E_1
real(rk) :: E_2
real(rk) :: E
real(rk) :: M
integer(ik) :: macro
real(rk) :: eta
real(rk) :: M_const
!! <h3> Variables describing the lattice switch </h3>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! <td> <b> Description </b> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>E</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The current energy. (<code>E=E_1</code> if <code>lattice=1</code> and <code>E=E_2</code> if <code>lattice=2</code>).
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>switchscalex</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The scalefactor of the supercell in the x-dimension to take the system from lattice 1 to lattice 2
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>switchscaley</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The scalefactor of the supercell in the y-dimension to take the system from lattice 1 to lattice 2
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>switchscalez</code> </font> </td>
!! <td> <font color="red"> <code>real(rk)</code> </font> </td>
!! <td> <font color="red">
!! The scalefactor of the supercell in the z-dimension to take the system from lattice 1 to lattice 2
!! </font> </td>
!! </tr>
real(rk) :: switchscalex
real(rk) :: switchscaley
real(rk) :: switchscalez
!! <h3> Variables keeping track of move numbers, acceptance rates, etc. </h3>
!! <table border="1">
!! <tr>
!! <td> <b> Variable </b> </td>
!! <td> <b> Type </b> </td>
!! <td> <b> Description </b> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>M_counts_1</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik), dimension(:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! <code>M_counts_1(i)</code> is the number of microstates visited which belong to
!! macrostate 'i' and lattice 1. This array must have <code>M_grid_size</code> elements.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>M_counts_2</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik), dimension(:), allocatable</code> </font> </td>
!! <td> <font color="red">
!! <code>M_counts_2(i)</code> is the number of microstates visited which belong to
!! macrostate 'i' and lattice 2. This array must have <code>M_grid_size</code> elements.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>sweeps</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">
!! Total number of sweeps which have been made.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>moves</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">
!! Total number of moves which have been made.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>moves_lattice</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">
!! Total number of lattice moves which have been made.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>accepted_moves_lattice</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">
!! Total number of accepted lattice switch moves.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>moves_part</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">
!! Total number of particle moves which have been made.
!! </font> </td>
!! </tr>
!! <tr>
!! <td> <font color="red"> <code>accepted_moves_part</code> </font> </td>
!! <td> <font color="red"> <code>integer(ik)</code> </font> </td>
!! <td> <font color="red">