Dice Fortran Backend Documentation
Functions/Subroutines
driver.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

program main
 Main program, contains the driver subroutine that runs all other code. More...
 
subroutine driver ()
 Subroutine for running the main code (needed for Doxygen compatability). More...
 

Function/Subroutine Documentation

◆ main()

program main

Main program, contains the driver subroutine that runs all other code.

Packaged as a separate subroutine within this program so that it remains compatible with Doxygen, which otherwise ignores any further documentation.

Definition at line 14 of file driver.f90.

+ Here is the call graph for this function:

◆ driver()

subroutine main::driver

Subroutine for running the main code (needed for Doxygen compatability).

The main driver code for the VQMC program. Links all other modules to run VQMC calculations on the Quantum Harmonic Oscillator (QHO), the Hydrogen Molecular Ion (H2plus) and the Hydrogen Molecule (H2).

Reads inputs from a provided input file, either generated by a Python frontend script or handwritten in the Fortran Namelist I/O format. Determines the problem to be solved, then runs VQMC routines. These can be run over a set parameter space, or in the case of the H2plus and H2 problems, can automatically determine the best variational parameter for the current bond length. Outputs results in NetCDF format, to be read by a Python frontend visualisation notebook for analysis.

Also includes extensions to start separate equilibration runs for determining optimal burning/thinning parameters for the Markov chain, and to write the energy traces for each parameter / parameter combination. Also capable of logging the full Markov chain to file (slow, should only be used for debug).

Definition at line 53 of file driver.f90.

54 
55 ! ---------------------------------------------------------------------------------------
56 ! SECTION: Variable Declarations
57 ! ---------------------------------------------------------------------------------------
58 
59  real(real64), dimension(:), allocatable :: p_arr
60  real(real64), dimension(:), allocatable :: bonds
61  real(real64), dimension(:), allocatable :: H2plus_param_c_array
62  real(real64), dimension(:), allocatable :: H2_param_a_array
63  real(real64), dimension(:), allocatable :: H2_param_beta_array
64  real(real64), dimension(:), allocatable :: energies
65  real(real64), dimension(:), allocatable :: variances
66  real(real64), dimension(:), allocatable :: accept_rates
68  real(real64) :: H2plus_param_c_guess
69  real(real64) :: H2plus_param_c_opt
70  real(real64) :: H2_param_beta_guess
71  real(real64) :: H2_param_beta_opt
72  real(real64) :: H2_param_a
73  real(real64) :: energy
74  real(real64) :: variance
75  real(real64) :: accept_rate
76 
77  integer :: b, i
78  integer :: logiostat
79  integer :: ierr
80  character(len=1) :: H2plus_c_converged
81  character(len=1) :: H2_beta_converged
82 
83 ! ---------------------------------------------------------------------------------------
84 ! SECTION: Initial Setup
85 ! ---------------------------------------------------------------------------------------
86  ! Initialise MPI
87  ! Initialize MPI
88  call mpi_init(mpierr)
89 
90  ! Find number of processors and current rank
91  call mpi_comm_size(mpi_comm_world,p,mpierr)
92  call mpi_comm_rank(mpi_comm_world,my_rank,mpierr)
93 
94  ! Read the input parameters.
95  call read_inputs()
96 
97  ! File management from rank 0
98  if (my_rank==0) then
99 
100  ! Notify user of the system
101  write(*, "('System: ', A)") p_system
102  write(*, "('')")
103 
104  ! Make sure the equilibration output directory exists if needed.
105  ! There is no simple way to check directory existence that works
106  ! with both ifort and gfortran, so mkdir is just run with stderr
107  ! suppressed instead.
108  if (run_equil) then
109  call system('mkdir ./equil_out 2>/dev/null')
110  write(*, "('Dice QMC started in Equilibration Mode')")
111  ! Check step length
112  if(run_steps>max_chain_len) then
113  write(*,"('')")
114  write(*, "(A, I7, A)") 'WARNING: The program will be slow at writing equilibration files, &
115  &as the requested number of steps exceeds the recommended &
116  &upper limit of ', max_chain_len, ' steps. Refer to the user manual for details.'
117  write(*,"('')")
118  endif
119  end if
120 
121  ! Same directory check as above for write_chains.
122  if (write_chains) then
123  if (.not. run_equil) then
124  call system('mkdir ./chains_out 2>/dev/null')
125  write(*, "('Dice QMC started in full VQMC mode with additional Markov chain outputting')")
126  ! Check step length
127  if(run_steps>max_chain_len) then
128  write(*,*) 'NOTE: For efficiency, only the first ', max_chain_len,' entries of the Markov chain &
129  &(after burning and thinning) will be written out to the chains_out folder.'
130  endif
131  else
132  error stop "Dice QMC cannot be started with both run_equil and write_chains set &
133  &to TRUE. Check your input params.txt and try again."
134  end if
135  end if
136 
137  ! Initialise the log file if required.
138  if (write_log) then
139 
140  ! Open the logfile
141  open(logid, file=logfile, status='unknown', access='sequential', &
142  form='formatted', iostat=logiostat)
143  if (logiostat .ne. 0) write(0, *) 'Error opening log file ', logfile
144 
145  ! Write program details to logfile
146  write(logid, "('*** Dice QMC log file ***')") ! include date/time?
147  write(logid, "('-------------------------------')")
148  write(logid, *) 'System: ', p_system
149  write(logid, "('')")
150 
151  ! Update log-writer tag for rank 0
152  i_log = .true.
153 
154  endif
155 
156  endif
157 
158  ! If p_system='test', run test module.
159  if (trim(p_system)=='test') then
160  ! Tests need to run on a single processor
161  if (my_rank==0) then
162  call unit_tests()
163  endif
164  ! Exit the program
165  call mpi_finalize(mpierr)
166  call exit()
167  endif
168 
169  ! If not a test job, setup MPI worksharing.
170  ! Equilibrations will be run serially on rank 0.
171  ! For non-equilibration runs:
172  if (.not. run_equil) then
173 
174  ! Split workload for VQMC runs
175  my_steps = run_steps / p ! Integer division
176  if (my_rank < (run_steps-(my_steps*p))) then
177  my_steps = my_steps+1
178  endif
179 
180  ! Pre-calculate frequently used factors.
181  my_size = my_steps/tstep ! Integer division
182  call mpi_reduce(my_size,tot_size,1,mpi_integer,mpi_sum,0,mpi_comm_world,mpierr)
183  if (my_rank==0) then
184  ! Factor representing division by total chain size after burning & thinning
185  per_step = 1.0_real64 / real(tot_size,kind=real64)
186  ! Factor representing division by number of processors
187  per_proc = 1.0_real64 / real(p,kind=real64)
188  ! Factor representing division by square of number of processors
189  per_proc2 = per_proc**2
190  endif
191  ! The factor 'per_proc' is required on all ranks for auto_grad calculations.
192  ! To ensure floating point agreement, broadcaste it from rank 0, rather than re-calculating on every rank.
193  call mpi_bcast(per_proc,1,mpi_double_precision,0,mpi_comm_world,mpierr)
194 
195  endif
196 
197  ! Initialise the random number generator
198  if (run_equil) then
199  ! Run calculations on one processor (rank 0)
200  p = 1
201  if(my_rank==0) then
202  ! Initialise the seed on rank 0
203  call init_ran1()
204  ! Update logfile if required
205  if (write_log) then
206  write(logid, "('Initial seed from /dev/urandom: ', I20)") seed
207  write(logid, "('')")
208  endif
209  ! Ensure auto_params is disabled
210  if (trim(p_system) .eq. 'H2plus' .and. h2plus%auto_params) then
211  h2plus%auto_params = .false.
212  write(*, "('WARNING: H2plus%auto_params should not be .TRUE. when in Equilibration Mode, disabling...')")
213  else if (trim(p_system) .eq. 'H2' .and. h2%auto_params) then
214  h2%auto_params = .false.
215  write(*, "('WARNING: H2plus%auto_params should not be .TRUE. when in Equilibration Mode, disabling...')")
216  end if
217  endif
218  else
219  ! For non-equil VQMC runs,
220  ! Initialise a seed on each MPI rank
221  call init_ran1()
222  ! Update logfile if required
223  if (write_log) then
224  if (my_rank>0) then
225  ! Send the seeds to rank 0
226  call mpi_send(seed,1,mpi_integer8,0,tag_seed,mpi_comm_world,mpierr)
227  else
228  ! Write seed of rank=0 to logfile
229  write(logid, "('Rank ',I2,'. Initial seed from /dev/urandom: ', I20)") my_rank, seed
230  do ip = 1, p-1
231  ! Receive seed value from rank=ip
232  call mpi_recv(ip_seed,1,mpi_integer8,ip,tag_seed,mpi_comm_world,status,mpierr)
233  ! Write seed of rank=ip to logfile
234  write(logid, "('Rank ',I2,'. Initial seed from /dev/urandom: ', I20)") ip, ip_seed
235  end do
236  write(logid, "('')")
237  endif
238  endif
239 
240  endif
241 
242 ! ---------------------------------------------------------------------------------------
243 ! SECTION: VQMC Loops
244 ! ---------------------------------------------------------------------------------------
245 
246  ! Execute problem-dependent routines.
247  select case(trim(p_system))
248 
249  ! -----------------------------------------------------------------------------------
250  ! SUBSECTION: Quantum Harmonic Oscillator
251  ! -----------------------------------------------------------------------------------
252 
253  case('QHO')
254 
255  ! Print out system info.
256  if (my_rank==0) then
257  write(*, "('Doing VQMC runs across ', I3, ' parameter(s).')") int(qho%alpha(3))
258  write(*, "('Steps: ', I10, ', Burn Steps: ', I7, ', Thin Every ', I7, ' Steps.')") &
259  qho%steps, qho%burn_step, qho%thin_step
260  write(*, "('')")
261  endif
262 
263  ! Set up alpha array.
264  p_arr = linspace(qho%alpha(1), qho%alpha(2), qho%alpha(3))
265 
266  if(my_rank==0) then
267  write(*, "('Grid search across linearly spaced parameters enabled.')")
268  write(*, "('')")
269  if (write_log) then
270  write(logid, "('Grid search across linearly spaced parameters enabled.')")
271  write(logid, "('Array of alpha values: ')")
272  do i = 1, int(qho%alpha(3))
273  write(logid, "(F12.8)") p_arr(i)
274  end do
275  write(logid, "('')")
276  end if
277  endif
278 
279  ! Main VQMC callers.
280  if (run_equil) then
281  ! For restarts, input_parser would have already set alpha_start
282  if (my_rank==0) then
283  call equil_qho(p_arr)
284  endif
285  else
286  ! Initialise arrays
287  allocate(energies(int(qho%alpha(3)))); energies(:) = 0.0_real64
288  allocate(variances(int(qho%alpha(3)))); variances(:) = 0.0_real64
289  allocate(accept_rates(int(qho%alpha(3)))); accept_rates(:) = 0.0_real64
290 
291  ! Sort out restarting
292  if (run_restart) then
293  ! Insert optimised parameters, energies, accept_rates, and variances into arrays
294  call read_energies(resfile_etot, energies, variances, accept_rates)
295  endif
296 
297  ! Run the grid evaluation
298  call grid_qho(p_arr, energies, accept_rates, variances)
299  end if
300 
301  ! -----------------------------------------------------------------------------------
302  ! SUBSECTION: Hydrogen Molecular Ion
303  ! -----------------------------------------------------------------------------------
304 
305  case('H2plus')
306 
307  ! Initialise arrays
308  allocate(energies(int(h2plus%bond(3))))
309  allocate(variances(int(h2plus%bond(3))))
310  allocate(accept_rates(int(h2plus%bond(3))))
311  allocate(h2plus_param_c_array(int(h2plus%bond(3))))
312 
313  ! Print out system info.
314  if (my_rank==0) then
315  write(*, "('Doing VQMC runs across ', I7, ' bond lengths.')") int(h2plus%bond(3))
316  write(*, "('Steps: ', I10, ', Burn Steps: ', I7, ', Thin Every ', &
317  & I7, ' Steps.')") h2plus%steps, h2plus%burn_step, h2plus%thin_step
318  write(*, "('')")
319  endif
320 
321  ! Set up bond length array
322  bonds = linspace(h2plus%bond(1), h2plus%bond(2), h2plus%bond(3))
323  if (i_log) then
324  write(logid, "('Array of bond lengths: ')")
325  do i = 1, int(h2plus%bond(3))
326  write(logid, "('Bond length ', I4, ': ', F7.4)") i, bonds(i)
327  end do
328  write(logid, "('')")
329  end if
330 
331  ! Set up parameter array, if doing a grid search
332  if (h2plus%auto_params .eqv. .false.) then
333  p_arr = linspace(h2plus%c_grid(1), h2plus%c_grid(2), h2plus%c_grid(3))
334  endif
335 
336  ! Update user from rank 0
337  if (my_rank==0) then
338  if (h2plus%auto_params .eqv. .false.) then
339  write(*, "('Grid search across linearly spaced parameters enabled.')")
340  write(*, "('')")
341  if (write_log) then
342  write(logid, "('Grid search across linearly spaced parameters enabled.')")
343  write(logid, "('Array of c values: ')")
344  do i = 1, int(h2plus%c_grid(3))
345  write(logid, "(F12.8)") p_arr(i)
346  end do
347  write(logid, "('')")
348  end if
349  else
350  write(*, "('Automatic search across parameters enabled.')")
351  write(*, "('')")
352  if (write_log) then
353  write(logid, "('Automatic search across parameters enabled.')")
354  write(logid, "('')")
355  end if
356  endif
357  endif
358 
359  ! Sort out restarting
360  if (run_restart .and. (.not. run_equil)) then
361  ! Insert optimised parameters, energies, accept_rates, and variances into arrays
362  call read_energies_h2plus(resfile_etot, energies, variances, accept_rates,h2plus_param_c_array)
363  endif
364 
365  ! Set starting guess value
366  if ((.not. run_equil) .and. h2plus%auto_params) then
367  if ( run_restart) then
368  if (.not. res%auto) error stop 'Error: Inconsistency in auto_params setting between params.txt and res_nml.txt.'
369  h2plus_param_c_guess = res%guess
370  else
371  h2plus_param_c_guess = h2plus%c
372  endif
373  endif
374 
375  ! Set up printed output for the loop.
376  if (my_rank==0) then
377  write(*, "('------------------ Begin VQMC Runs ------------------')")
378  write(*, "('Bond Length Converged? Optimised c Energy')")
379  if (write_log) then
380  write(logid, "('------------------ Begin VQMC Runs ------------------')")
381  end if
382  endif
383 
384  ! Loop over bond lengths
385  do b = bond_start, int(h2plus%bond(3))
386 
387  ! Store bond index for restart writing
388  i_bond = b
389 
390  if (i_log) then
391  write(logid, "(' ------------- Bond Length = ', F7.4, ' ------------- ')") &
392  bonds(b)
393  end if
394 
395  ! Main VQMC callers
396  if (.not. run_equil) then
397 
398  if (h2plus%auto_params) then
399 
400  ! Search for best param at this bond length
401  call auto_h_2_plus(h2plus_param_c_guess, bonds(b), h2plus_param_c_opt, &
402  energy, variance, accept_rate, h2plus_c_converged)
403 
404  ! Reset start counter
405  auto_start=1
406 
407  ! Use current optimised value of c as guess for next bond length's c,
408  ! provided calculations converged this loop
409  if (h2plus_c_converged .eq. 'Y') h2plus_param_c_guess = h2plus_param_c_opt
410 
411  else
412 
413  ! Find best param in grid at this bond length
414  call grid_h_2_plus(p_arr, bonds(b), h2plus_param_c_opt, energy, &
415  variance, accept_rate)
416 
417  ! Reset start counter
418  grid_start=1
419 
420  end if
421 
422  if(my_rank==0) then
423 
424  ! Append energy and parameters to arrays, on rank 0
425  h2plus_param_c_array(b) = h2plus_param_c_opt
426  energies(b) = energy
427  variances(b) = variance
428  accept_rates(b) = accept_rate
429 
430  ! Write restart files, if required
431  if (write_restart) then
432  ! Store the arrays
433  call write_h2plus_main(h2plus_param_c_array, bonds, energies, variances, &
434  resfile_etot, ierr, h2plus, accept_rates)
435  ! Write res_nml.txt
436  if (h2plus%auto_params) then
437  call write_res_nml(b+1,1,h2plus%auto_params,h2plus_param_c_guess)
438  else
439  call write_res_nml(b+1,1,h2plus%auto_params)
440  endif
441  endif
442 
443  endif
444 
445  else
446  ! Run equilibrations serially on rank 0
447  if (my_rank==0) then
448  call equil_h_2_plus(p_arr, bonds(b), b)
449  ! Write restart files, if required
450  if (write_restart) then
451  call write_res_nml(b+1,1,h2plus%auto_params)
452  endif
453 
454  endif
455  end if
456 
457  end do
458 
459  if (my_rank==0) then
460 
461  write(*, "('------------------- End VQMC Runs -------------------')")
462  write(*, "('')")
463 
464  ! If not in equilibration mode, write results
465  if (.not. run_equil) then
466  write(*, "(2A)") 'Calculations finished, writing results to NetCDF output file ', &
467  outfile
468  if (write_log) then
469  write(logid, "(2A)") 'Calculations finished, writing results to NetCDF &
470  &output file ', outfile
471  end if
472 
473  call write_h2plus_main(h2plus_param_c_array, bonds, energies, variances, &
474  outfile, ierr, h2plus, accept_rates)
475  end if
476 
477  endif
478 
479 
480  ! -----------------------------------------------------------------------------------
481  ! SUBSECTION: Hydrogen Molecule
482  ! -----------------------------------------------------------------------------------
483 
484  case('H2')
485 
486  ! Initialise arrays
487  allocate(energies(int(h2%bond(3))))
488  allocate(variances(int(h2%bond(3))))
489  allocate(accept_rates(int(h2%bond(3))))
490  allocate(h2_param_a_array(int(h2%bond(3))))
491  allocate(h2_param_beta_array(int(h2%bond(3))))
492 
493  ! Print out system info.
494  if (my_rank==0) then
495  write(*, "('Doing VQMC runs across ', I7, ' bond lengths.')") int(h2%bond(3))
496  write(*, "('Steps: ', I10, ', Burn Steps: ', I7, ', Thin Every ', I7, ' Steps.')") &
497  h2%steps, h2%burn_step, h2%thin_step
498  write(*, "('')")
499  endif
500 
501  ! Set up bond length array
502  bonds = linspace(h2%bond(1), h2%bond(2), h2%bond(3))
503  if (i_log) then
504  write(logid, "('Array of bond lengths: ')")
505  do i = 1, int(h2%bond(3))
506  write(logid, "('Bond length ', I4, ': ', F7.4)") i, bonds(i)
507  end do
508  write(logid, "('')")
509  end if
510 
511  ! Set up parameter array, if doing a grid search
512  if (h2%auto_params .eqv. .false.) then
513  p_arr = linspace(h2%beta_grid(1), h2%beta_grid(2), h2%beta_grid(3))
514  endif
515 
516  ! Update user from rank 0
517  if (my_rank==0) then
518  if (h2%auto_params .eqv. .false.) then
519  write(*, "('Grid search across linearly spaced parameters enabled.')")
520  write(*, "('')")
521  if (write_log) then
522  write(logid, "('Grid search across linearly spaced parameters enabled.')")
523  write(logid, "('Array of beta values: ')")
524  do i = 1, int(h2%beta_grid(3))
525  write(logid, "(F12.8)") p_arr(i)
526  end do
527  write(logid, "('')")
528  end if
529  else
530  write(*, "('Automatic search across parameters enabled.')")
531  write(*, "('')")
532  if (write_log) then
533  write(logid, "('Automatic search across parameters enabled.')")
534  write(logid, "('')")
535  end if
536  endif
537  end if
538 
539  ! Sort out restarting
540  if (run_restart .and. (.not. run_equil)) then
541  ! Insert optimised parameters, energies, accept_rates, and variances into arrays
542  call read_energies_h2(resfile_etot, energies, variances, accept_rates,h2_param_a_array,h2_param_beta_array)
543  endif
544 
545  ! Set starting guess value
546  if ((.not. run_equil) .and. h2%auto_params) then
547  if (run_restart) then
548  if (.not. res%auto) error stop 'Error: Inconsistency in auto_params setting between params.txt and res_nml.txt.'
549  h2_param_beta_guess = res%guess
550  else
551  h2_param_beta_guess = h2%beta
552  endif
553  endif
554 
555  ! Set up printed output for the loop.
556  if (my_rank==0) then
557  write(*, "('------------------ Begin VQMC Runs ------------------')")
558  write(*, "('Bond Length Converged? Optimised Beta Energy')")
559  if (i_log) then
560  write(logid, "('------------------ Begin VQMC Runs ------------------')")
561  end if
562  endif
563 
564  ! Loop over bond lengths
565  do b = bond_start, int(h2%bond(3))
566 
567  ! Store bond index for restart writing
568  i_bond = b
569 
570  if (i_log) then
571  write(logid, "(' ------------- Bond Length = ', F7.4, ' ------------- ')") &
572  bonds(b)
573  end if
574 
575  ! Calculate parameter a at this bond length by Newton-Raphson iteration.
576  h2_param_a = h_2_cusp(bonds(b))
577 
578  ! Main VQMC callers
579  if (.not. run_equil) then
580 
581  ! Either iterate through beta values, or automatically search for best beta.
582  if (h2%auto_params) then ! ie. if auto searching
583 
584  ! Solve value of beta at this bond length.
585  call auto_h_2(h2_param_a, h2_param_beta_guess, bonds(b), h2_param_beta_opt, &
586  energy, variance, accept_rate, h2_beta_converged)
587 
588  ! Use current optimised value of beta as guess for next bond length's beta,
589  ! provided calculations converged this loop
590  if (h2_beta_converged .eq. 'Y') h2_param_beta_guess = h2_param_beta_opt
591 
592  ! Reset start counter
593  auto_start=1
594 
595  else ! ie. if doing a grid search
596 
597  ! Get best value of beta at this bond length
598  call grid_h_2(h2_param_a, p_arr, bonds(b), h2_param_beta_opt, energy, &
599  variance, accept_rate)
600 
601  ! Reset start counter
602  grid_start=1
603 
604  ! Write restart files, if required
605  if (write_restart .and. my_rank==0) then
606  call write_h2_main(h2_param_a_array, h2_param_beta_array, bonds, &
607  energies, variances, resfile_etot, ierr, h2, accept_rates)
608  call write_res_nml(b+1,1,h2%auto_params)
609  endif
610 
611  end if
612 
613  if (my_rank==0) then
614 
615  ! Append energy and parameters to arrays
616  h2_param_a_array(b) = h2_param_a
617  h2_param_beta_array(b) = h2_param_beta_opt
618  energies(b) = energy
619  variances(b) = variance
620  accept_rates(b) = accept_rate
621 
622  ! Write restart files, if required
623  if (write_restart) then
624  ! Store the arrays
625  call write_h2_main(h2_param_a_array, h2_param_beta_array, bonds, &
626  energies, variances, resfile_etot, ierr, h2, accept_rates)
627  if (h2%auto_params) then
628  call write_res_nml(b+1,1,h2%auto_params,h2_param_beta_guess)
629  else
630  call write_res_nml(b+1,1,h2%auto_params)
631  endif
632  endif
633 
634  endif
635 
636  else
637 
638  ! Run equilibrations
639  if (my_rank==0) then
640  call equil_h_2(h2_param_a, p_arr, bonds(b), b)
641  ! Write restart files, if required
642  if (write_restart) then
643  call write_res_nml(b+1,1,h2%auto_params)
644  endif
645  endif
646 
647  end if
648 
649  end do
650 
651  if (my_rank==0) then
652 
653  write(*, "('------------------- End VQMC Runs -------------------')")
654  write(*, "('')")
655 
656  ! If not in equilibration mode, write results
657  if (.not. run_equil) then
658 
659  write(*, "(2A)") 'Calculations finished, writing results to NetCDF output file ', &
660  outfile
661 
662  if (write_log) then
663  write(logid, "(2A)") 'Calculations finished, writing results to NetCDF &
664  &output file ', outfile
665  end if
666 
667  call write_h2_main(h2_param_a_array, h2_param_beta_array, bonds, &
668  energies, variances, outfile, ierr, h2, accept_rates)
669 
670  end if
671 
672  endif
673 
674  ! Default case, for if none of the cases match. Should be invalidated by checks in input_parser.
675  case default
676  error stop 'Unknown case. How did you get here?'
677  end select
678 
679  ! Close logfile from rank 0
680  if (i_log) close(logid)
681 
682  ! Remove restart files if complete
683  call system('rm res_etot.nc 2>/dev/null')
684  call system('rm res_nml.txt 2>/dev/null')
685  call system('rm res_epar.nc 2>/dev/null')
686 
687  ! End MPI parallelism
688  call mpi_finalize(mpierr)
689 
+ Here is the call graph for this function:
+ Here is the caller graph for this function: