BraWl
Loading...
Searching...
No Matches
howto_examples Module Reference

Functions/Subroutines

subroutine examples (setup, metropolis, my_rank)
 Subroutine showing examples of code's functionalities.
 

Function/Subroutine Documentation

◆ examples()

subroutine howto_examples::examples ( type(run_params) setup,
type(metropolis_params) metropolis,
integer, intent(in) my_rank )

Subroutine showing examples of code's functionalities.

Author
C. D. Woodgate
Date
2019-2025
Parameters
setupDerived type containing simulation parameters
my_rankRank of current MPI process
Returns
None

Definition at line 43 of file howto_examples.f90.

44
45 ! Rank of this processor
46 integer, intent(in) :: my_rank
47
48 ! Arrays for storing data
49 type(run_params) :: setup
50
51 ! Arrays for storing data
52 type(metropolis_params) :: metropolis
53
54 ! Integers used in calculations
55 integer :: i,j, div_steps, accept, n_save
56
57 ! Temperature and temperature steps
58 real(real64) :: beta, temp, sim_temp, current_energy, step_E, &
59 step_Esq, C, acceptance
60
61 ! Name of file for grid state and xyz file at this temperature
62 character(len=36) :: xyz_file
63
64 ! Name of file for writing diagnostics at the end
65 character(len=43) :: diagnostics_file
66
67 ! Name of file for writing radial densities at the end
68 character(len=37) :: radial_file
69
70 ! Radial densities array (to populate)
71 real(real64), allocatable, dimension(:,:,:) :: r_densities
72
73 ! This is for calculating radial density functions later
74 ! (It populates the array "shells".)
75 call lattice_shells(setup, shells, config)
76
77 ! Are we swapping just nearest-neighbours or
78 ! allowing any pair of lattice sites?
79 if (metropolis%nbr_swap) then
80 setup%mc_step => monte_carlo_step_nbr
81 else
82 setup%mc_step => monte_carlo_step_lattice
83 end if
84
85 !----------!
86 ! EXAMPLES !
87 !----------!
88
89 !--------------------------------------------!
90 ! 1. Initialise a random grid and display it !
91 !--------------------------------------------!
92
93 ! Set up the initial state of the grid (random occupancies)
94 ! on this rank. (Each rank works on different grid currently.)
95 call initial_setup(setup, config)
96
97 if(my_rank == 0) then
98 print*, ' '
99 print*, ' This is what the grid initially looks like: '
100 print*, ' '
101 end if
102
103 ! You can see what it looks like (sort of) using this routine
104 call display_grid(config)
105
106 !---------------------------------!
107 ! 2. Perform a Metropolis MC step !
108 !---------------------------------!
109
110 ! Would normally use j to loop over temperatures
111 j = 1
112
113 ! Work out the temperature and corresponding beta
114 temp = metropolis%T + real(j-1, real64)*metropolis%delta_T
115 sim_temp = temp*k_b_in_ry
116 beta = 1.0_real64/sim_temp
117
118 accept=0.0_real64
119
120 do while (accept .lt. 0.5_real64)
121 accept = setup%mc_step(config, beta)
122 end do
123
124 if(my_rank == 0) then
125 print*, ' '
126 print*, ' This is what the grid looks like after one successful step: '
127 print*, ' '
128 end if
129
130 call display_grid(config)
131
132 !---------------------------------------------------!
133 ! 3. Get the energy associated with a configuration !
134 !---------------------------------------------------!
135
136 current_energy = setup%full_energy(config)
137
138 if(my_rank == 0) then
139 print*, ' '
140 print*, ' The energy of that configuration is: ', current_energy, ' Ry'
141 print*, ' or: ', current_energy*13.606*1000.0, ' meV'
142 print*, ' or: ', current_energy*13.606*1000.0/setup%n_atoms, ' meV/atom'
143 print*, ' '
144 end if
145
146 !------------------------------------!
147 ! 4. Run some kind of equillibration !
148 !------------------------------------!
149 n_save=floor(real(metropolis%n_mc_steps)/real(metropolis%n_sample_steps))
150
151 acceptance = 0.0_real64
152 step_e = 0.0_real64
153 step_esq = 0.0_real64
154
155 do i=1, metropolis%n_mc_steps
156
157 ! Make one MC move
158 accept = setup%mc_step(config, beta)
159
160 acceptance = acceptance + accept
161
162 ! Write percentage progress to screen
163 if (mod(i, metropolis%n_sample_steps) .eq. 0) then
164 current_energy = setup%full_energy(config)
165 step_e = step_e + current_energy
166 step_esq = step_esq + current_energy**2
167 end if
168
169 end do
170
171 ! Store the average energy per atom at this temperature
172 energies_of_t(j) = step_e/n_save/setup%n_atoms
173
174 ! Heat capacity (per atom) at this temperature
175 c = (step_esq/n_save - (step_e/n_save)**2)/(sim_temp*temp)/setup%n_atoms
176
177 ! Acceptance rate at this temperature
178 acceptance_of_t(j) = acceptance/real(metropolis%n_mc_steps)
179
180 ! Store the specific heat capacity at this temperature
181 c_of_t(j) = c
182
183 if (my_rank ==0) then
184 ! Write that we have completed a particular temperature
185 write(6,'(a,f7.2,a)',advance='yes') &
186 " Temperature ", temp, " complete on process 0."
187 write(6,'(a,f7.2,a)',advance='yes') &
188 " Internal energy was ", 13.606_real64*1000*energies_of_t(j), " meV/atom"
189 write(6,'(a,f9.4,a)',advance='yes') &
190 " Heat capacity was ", c_of_t(j)/k_b_in_ry, " kB/atom"
191 write(6,'(a,f7.2,a,/)',advance='yes') &
192 " Swap acceptance rate was ", 100.0*acceptance_of_t(j), "%"
193 end if
194
195 !------------------------------------!
196 ! 5. Write the grid to a .xyz file !
197 !------------------------------------!
198
199 write(xyz_file, '(A11 I3.3 A12 I4.4 F2.1 A4)') 'configs/proc_', &
200 my_rank, 'config_at_T_', int(temp), temp-int(temp),'.xyz'
201
202 ! Write xyz file
203 call xyz_writer(xyz_file, config, setup)
204
205 !----------------------------------------------------------------!
206 ! 6. Compute the Warren-Cowley ASRO parameters and write to file !
207 !----------------------------------------------------------------!
208
209 ! Work out how to find each shell
210 call lattice_shells(setup, shells, config)
211
212 ! Compute the radial densities
213 r_densities = radial_densities(setup, config, setup%wc_range, shells)
214
215 write(radial_file, '(A,I3.3,A)') 'asro/proc_', my_rank, '_rho_of_T.nc'
216
217 ! Write the radial densities to file
218 call ncdf_radial_density_writer_once(trim(radial_file), r_densities, shells, setup)
219
220 if(my_rank == 0) then
221 write(6,'(25("-"),x,"Simulation Complete!",x,25("-"))')
222 end if
223
Here is the call graph for this function:
Here is the caller graph for this function: