Subroutine showing examples of code's functionalities.
44
45
46 integer, intent(in) :: my_rank
47
48
49 type(run_params) :: setup
50
51
52 type(metropolis_params) :: metropolis
53
54
55 integer :: i,j, div_steps, accept, n_save
56
57
58 real(real64) :: beta, temp, sim_temp, current_energy, step_E, &
59 step_Esq, C, acceptance
60
61
62 character(len=36) :: xyz_file
63
64
65 character(len=43) :: diagnostics_file
66
67
68 character(len=37) :: radial_file
69
70
71 real(real64), allocatable, dimension(:,:,:) :: r_densities
72
73
74
75 call lattice_shells(setup, shells, config)
76
77
78
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
87
88
89
90
91
92
93
94
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
104 call display_grid(config)
105
106
107
108
109
110
111 j = 1
112
113
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
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
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
158 accept = setup%mc_step(config, beta)
159
160 acceptance = acceptance + accept
161
162
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
172 energies_of_t(j) = step_e/n_save/setup%n_atoms
173
174
175 c = (step_esq/n_save - (step_e/n_save)**2)/(sim_temp*temp)/setup%n_atoms
176
177
178 acceptance_of_t(j) = acceptance/real(metropolis%n_mc_steps)
179
180
181 c_of_t(j) = c
182
183 if (my_rank ==0) then
184
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
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
203 call xyz_writer(xyz_file, config, setup)
204
205
206
207
208
209
210 call lattice_shells(setup, shells, config)
211
212
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
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