Dice Fortran Backend Documentation
vqmc.f90
Go to the documentation of this file.
1 ! --------------------------------------------------------------------------------------------------------------------
2 ! Dice Quantum Monte Carlo
3 ! --------------------------------------------------------------------------------------------------------------------
4 ! MODULE: VQMC
5 !
6 ! DESCRIPTION:
10 !
29 ! --------------------------------------------------------------------------------------------------------------------
30 module vqmc
31 
32  use iso_fortran_env
33  use shared_data
34  use constants
35  use solvers
36  use write_netcdf
37  use restart_fns
38 
39  implicit none
40 
41  !public seed ! Allow testing module access so seed can be fixed for testing.
42 
43  ! Variables which control the RNG sequence
44  integer, parameter :: ntab = 32
45  integer(int64), dimension(ntab), save :: iv = 0
46  integer(int64), save :: iy = 0
47  !integer(int64) :: seed !< Seed for the random number generator.
48 
49  integer :: logfrac = 10
50 
51 contains
52 
53 ! --------------------------------------------------------------------------------------------------------------------
54 ! SECTION: Random number generator
55 ! --------------------------------------------------------------------------------------------------------------------
56 
57  ! ----------------------------------------------------------------------------------------------------------------
58  ! ROUTINE: init_ran
59  !
60  ! DESCRIPTION:
63  !
67  ! ----------------------------------------------------------------------------------------------------------------
68  subroutine init_ran1()
69 
70  implicit none
71 
72  integer :: fu ! File unit
73  integer :: rc ! I/O status of read from /dev/urandom
74  integer :: l_seed ! Length of unshortened seed
75  character(len=30) :: c_seed ! Character value of seed
76 
77  ! Read pseudorandom noise from /dev/urandom
78  open (access='stream', action='read', file='/dev/urandom', &
79  form='unformatted', iostat=rc, newunit=fu)
80 
81  if (rc == 0) read (fu) seed
82 
83  close (fu)
84 
85  ! Convert long seed to characters, shorten and convert back into an integer
86  write(c_seed, "(I30)") seed
87  l_seed = len(c_seed)
88  c_seed = adjustl(c_seed)
89  c_seed = c_seed(7:l_seed)
90  read(c_seed, *) seed
91 
92  ! Ensure seed is negative for compatability with ran1
93  if (seed .gt. 0) then
94  seed = -seed
95  end if
96 
97  end subroutine init_ran1
98 
99  ! ----------------------------------------------------------------------------------------------------------------
100  ! ROUTINE: ran1
101  !
102  ! DESCRIPTION:
105  !
113  !
114  ! PARAMETERS:
117  ! ----------------------------------------------------------------------------------------------------------------
118  subroutine ran1(sample, dseed)
119 
120  implicit none
121 
122  integer(int64), intent(inout) :: dseed
123  real(real64), intent(out) :: sample
124 
125  ! Constants defined in Press et al
126  integer(int64), parameter :: ia = 16807, im = 2147483647, iq = 127773
127  integer(int64), parameter :: ir = 2836, &
128  ndiv = 1 + floor((real(im, kind=real64)-1.0)/ntab)
129  real(real64), parameter :: am = 1.0/real(im, kind=real64), eps = 1.2e-7, &
130  rnmx = 1.0-eps
131 
132  ! Loop counters
133  integer(int64) :: j, k
134 
135  ! The following is transcribed directly from Press et al
136  if (dseed<=0 .or. iy==0) then ! initialization
137  ! NB this may need to be modified to allow for restarting as iy will be 0
138  dseed = max(-dseed,1)
139 
140  do j = ntab+8, 1, -1
141  k = floor(real(dseed, kind=real64)/real(iq, kind=real64))
142  dseed = ia * (dseed-k*iq) - ir * k
143  if (dseed < 0) dseed = dseed + im
144  if (j <= ntab) iv(j) = dseed
145  end do
146 
147  iy = iv(1)
148  end if
149 
150  k = floor(real(dseed, kind=real64)/real(iq, kind=real64))
151  dseed = ia * (dseed-k*iq) - ir * k
152  if (dseed < 0) dseed = dseed + im
153  j = 1_int64 + iy / ndiv
154  iy = iv(j)
155  iv(j) = dseed
156  sample = min(am*iy, rnmx)
157 
158  end subroutine ran1
159 
160 ! --------------------------------------------------------------------------------------------------------------------
161 ! SECTION: Useful functions and subroutines
162 ! --------------------------------------------------------------------------------------------------------------------
163 
164  ! ----------------------------------------------------------------------------------------------------------------
165  ! ROUTINE: linspace
166  !
167  ! DESCRIPTION:
170  !
174  !
175  ! PARAMETERS:
179  !
181  ! ----------------------------------------------------------------------------------------------------------------
182  function linspace(start, end, num) result(ls_array)
183 
184  implicit none
185 
186  real(real64), intent(in) :: start, end, num
187 
188  real(real64), dimension(:), allocatable :: ls_array
189 
190  integer :: inum ! The integer type of num
191  integer :: i ! Loop variable
192  real(real64) :: step ! Spacing between two numbers in the array
193 
194  ! Handle real value for num
195  inum = int(num)
196  ! Handle single values
197  if (inum .eq. 1) then
198  allocate(ls_array(inum))
199  ls_array(1) = start
200  ! Handle multiple values
201  else
202  ! Generate array
203  allocate(ls_array(inum))
204  ! Find step size
205  step = (end - start) / real(inum-1, kind=real64)
206  ! Fill array
207  do i = 1, inum
208  ls_array(i) = start + (i-1)*step
209  end do
210  end if
211 
212  end function linspace
213 
214  ! ----------------------------------------------------------------------------------------------------------------
215  ! ROUTINE: random_normal
216  !
217  ! DESCRIPTION:
220  !
223  !
224  ! PARAMETERS:
227  ! ----------------------------------------------------------------------------------------------------------------
228  subroutine random_normal(dim, x)
229 
230  implicit none
231 
232  integer, intent(in) :: dim
233  real(real64), intent(out), dimension(:) :: x
234 
235  real(real64) :: u_1, u_2 ! Samples from uniform distribution
236  integer :: i ! Loop variable
237 
238  ! Generate vector containg samples from normal distributions
239  do i = 1, dim
240  ! This can be done the Fortran intrinsic random number generator but this
241  ! is not a good random number generator
242  ! call random_number(u_1)
243  ! call random_number(u_2)
244 
245  ! Draws from uniform distribution
246  call ran1(u_1, seed)
247  call ran1(u_2, seed)
248 
249  ! Calculate the sample from the normal distribution using the Box-Muller transform
250  x(i) = sqrt(-2*log(u_1))*cos(2*pi*u_2)
251  end do
252 
253  end subroutine random_normal
254 
255 ! --------------------------------------------------------------------------------------------------------------------
256 ! SECTION: VQMC for QHO
257 ! --------------------------------------------------------------------------------------------------------------------
258 
259  ! ----------------------------------------------------------------------------------------------------------------
260  ! ROUTINE: VQMC_QHO
261  !
262  ! DESCRIPTION:
265  !
274  !
275  ! PARAMETERS:
282  !
288  ! ----------------------------------------------------------------------------------------------------------------
289  subroutine vqmc_qho(sigma, N, burn, thin, alpha_const, X_burnt, &
290  accept_rate, energy_chain, energy_mean, energy_variance, x_start)
291 
292  implicit none
293 
294  real(real64), intent(in) :: sigma, alpha_const
295  integer, intent(in) :: N, burn, thin
296  real(real64), optional :: x_start
297 
298  real(real64), intent(out) :: accept_rate, energy_mean, energy_variance
299  real(real64), intent(out), dimension(:), allocatable :: X_burnt, energy_chain
300 
301  real(real64) :: log_prob_n ! The natural logarithm of the transition probability
302  real(real64) :: a_prob_trans ! Acceptance ratio of the transition probability
303  real(real64) :: unif_num ! A random sample from the uniform distribution
304  real(real64), dimension(1) :: start ! The starting point of the MCMC chain
305  real(real64), dimension(1) :: X_n ! The proposed step
306  real(real64), dimension(1) :: sample_normal ! A sample from the normal distribution
307  real(real64), dimension(:), allocatable :: X ! Non-burnt or thinned MCMC chain
308  integer :: thin_counter ! Counter to keep track of the thinning process
309  integer :: reset_counter ! Counts how many times thin counter was reset
310  integer :: accept ! Number of accepted steps
311  integer :: i ! Loop variable
312  integer :: i_bt ! Burned/thinned value of i
313  integer :: re_size ! The size of the burnt and thinned array
314  character(len=1) :: acceptance ! Yes/No signpost for acceptance
315 
316  ! Adjust the output array to be of the burnt and thinned array
317  re_size = (n - burn) / thin
318 
319  allocate(x(1: n+1)); x(:) = 0.0_real64
320  allocate(x_burnt(1: re_size)); x_burnt(:) = 0.0_real64
321  allocate(energy_chain(1: re_size)); energy_chain(:) = 0.0_real64
322 
323  ! Allow the starting position to take the defualt value
324  if (present(x_start)) then
325  start = x_start
326  x(1) = start(1)
327  if (i_log) then
328  write(logid, "('Starting from provided Markov chain position', F12.8)") x(1)
329  end if
330  else
331  ! Creates the initial position of the random walker that is within a certain radius
332  call random_normal(1, sample_normal)
333  do while(norm2(sample_normal) >= qho%rcut)
334  call random_normal(1, sample_normal)
335  end do
336  start = sample_normal
337  x(1) = start(1)
338  if (i_log) then
339  write(logid, "('Generated random Markov chain start position', F12.8)") x(1)
340  end if
341  end if
342 
343  ! Start the counter for the number of accepted moves
344  accept = 0
345  ! Start the counter for the thinning process
346  thin_counter = 0
347  ! Start the reset counter
348  reset_counter = 0
349 
350  !Write initial conditions to the output file
351  if (i_log) then
352  write(logid, "('Markov chain for alpha = ', F7.4)") alpha_const
353  write(logid, "('Sample Position Local Energy Accepted?')")
354  write(logid, "(I8, ' ', F12.8, ' ', F12.6)") 1, x(1), energy_chain(1)
355  end if
356 
357  ! In case burn = 0 to start calculating the chain from the beginning
358  if (burn == 0) then
359  energy_chain(1) = qho_energy(alpha_const, x(1))
360  x_burnt(1) = x(1)
361 
362  do i = 2, n
363  ! Add one to the thin counter
364  thin_counter = 1 + thin_counter
365  ! Propose the next step
366  call random_normal(1, sample_normal)
367  x_n = x(i-1) + sigma * sample_normal
368 
369  ! Calculate the transition probability between proposed step and current step
370  log_prob_n = log(qho_prob(alpha_const, x(i-1), x_n(1)))
371 
372  ! Calculate the acceptance ratio
373  a_prob_trans = min(1.0_real64, exp(log_prob_n))
374 
375  ! Generate sample from uniform distribution
376  call ran1(unif_num, seed)
377 
378  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
379  if (a_prob_trans >= unif_num) then
380  ! If accepted add proposed position to the MCMC chain
381  x(i) = x_n(1)
382 
383  ! To determine whether this point should be added to the burnt and thinned MCMC chain
384  if (thin_counter == thin) then
385  ! Increment the reset counter
386  reset_counter = 1 + reset_counter
387  ! Precalculate i_bt for shorthand
388  i_bt = i - burn - (reset_counter * (thin - 1))
389  ! Add to the burn and thinned MCMC chain
390  x_burnt(i_bt) = x_n(1)
391  ! Calculate the local energy at this point
392  energy_chain(i_bt) = qho_energy(alpha_const, x_n(1))
393  ! Resetting the thin counter
394  thin_counter = 0
395 
396  acceptance = 'Y'
397  ! Increment the number of accepted steps
398  accept = 1 + accept
399 
400  if (i_log) then
401  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
402  write(logid, "(I8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
403  (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
404  end if
405  end if
406  end if
407  else
408  ! If rejected the proposed step is the current step
409  x(i) = x(i-1)
410 
411  ! To determine whether this point should be added to the burnt and thinned MCMC chain
412  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
413  ! Increment the reset counter
414  reset_counter = 1 + reset_counter
415  ! Precalculate i_bt for shorthand
416  i_bt = i - burn - (reset_counter * (thin - 1))
417  ! Add to the burn and thinned MCMC chain
418  x_burnt(i_bt) = x(i-1)
419  ! Calculate the local energy at this point
420  energy_chain(i_bt) = qho_energy(alpha_const, x(i-1))
421  ! Reseting the thin counter
422  thin_counter = 0
423 
424  acceptance = 'N'
425 
426  if (i_log) then
427  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
428  write(logid, "(I8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
429  (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
430  end if
431  end if
432 
433  end if
434  end if
435 
436  end do
437  ! Loop to advance the chain to the final burn step and not save any of the outputs
438  else if (burn >= 2) then
439  do i = 2, burn
440  ! Propose the next step
441  call random_normal(1, sample_normal)
442  x_n = x(i-1) + sigma * sample_normal
443  ! Calculate the transition probability between proposed step and current step
444  log_prob_n = log(qho_prob(alpha_const, x(i-1), x_n(1)))
445  ! Calculate the acceptance ratio
446  a_prob_trans = min(1.0_real64, exp(log_prob_n))
447  ! Generate sample from uniform distribution
448  call ran1(unif_num, seed)
449 
450  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
451  if (a_prob_trans >= unif_num) then
452  ! If accepted add proposed position to the MCMC chain
453  x(i) = x_n(1)
454  else
455  ! If rejected the proposed step is the current step
456  x(i) = x(i-1)
457  end if
458  end do
459  end if
460 
461  ! Metropolis Algorithm to calculate the burnt and thinned MCMC chain
462  if (burn >=1) then
463  do i = burn + 1, n+1
464  ! Add one to the thin counter
465  thin_counter = 1 + thin_counter
466  ! Propose the next step
467  call random_normal(1, sample_normal)
468  x_n = x(i-1) + sigma * sample_normal
469 
470  ! Calculate the transition probability between proposed step and current step
471  log_prob_n = log(qho_prob(alpha_const, x(i-1), x_n(1)))
472 
473  ! Calculate the acceptance ratio
474  a_prob_trans = min(1.0_real64, exp(log_prob_n))
475 
476  ! Generate sample from uniform distribution
477  call ran1(unif_num, seed)
478 
479  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
480  if (a_prob_trans >= unif_num) then
481  ! If accepted add proposed position to the MCMC chain
482  x(i) = x_n(1)
483 
484  ! To determine whether this point should be added to the burnt and thinned MCMC chain
485  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
486  ! Increment the reset counter
487  reset_counter = 1 + reset_counter
488  ! Precalculate i_bt for shorthand
489  i_bt = i - burn - (reset_counter * (thin - 1))
490  ! Add to the burn and thinned MCMC chain
491  x_burnt(i_bt) = x_n(1)
492  ! Calculate the local energy at this point
493  energy_chain(i_bt) = qho_energy(alpha_const, x_n(1))
494  ! Resetting the thin counter
495  thin_counter = 0
496 
497  acceptance = 'Y'
498  ! Increment the number of accepted steps
499  accept = 1 + accept
500 
501  if (i_log) then
502  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
503  write(logid, "(I8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
504  (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
505  end if
506  end if
507 
508  end if
509 
510  else
511  ! If rejected the proposed step is the current step
512  x(i) = x(i-1)
513 
514  ! To determine whether this point should be added to the burnt and thinned MCMC chain
515  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
516  ! Increment the reset counter
517  reset_counter = 1 + reset_counter
518  ! Precalculate i_bt for shorthand
519  i_bt = i - burn - (reset_counter * (thin - 1))
520  ! Add to the burn and thinned MCMC chain
521  x_burnt(i_bt) = x(i-1)
522  ! Calculate the local energy at this point
523  energy_chain(i_bt) = qho_energy(alpha_const, x(i-1))
524  ! Reseting the thin counter
525  thin_counter = 0
526 
527  acceptance = 'N'
528 
529  if (i_log) then
530  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
531  write(logid, "(I8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
532  (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
533  end if
534  end if
535 
536  end if
537  end if
538  end do
539  end if
540 
541  ! Calculate the acceptance ratio
542  accept_rate = (real(accept, kind=real64) / real(re_size, kind=real64))
543 
544  ! Calculate the mean energy and the variance in the energy
545  energy_mean = sum(energy_chain) / real(re_size, kind=real64)
546  energy_variance = 0.0_real64
547  do i = 1, re_size
548  energy_variance = energy_variance + ((1.0_real64 / real(re_size - 1, kind=real64)) * &
549  (energy_chain(i) - energy_mean) * (energy_chain(i) - energy_mean))
550  end do
551  energy_variance = (1.0_real64 / real(re_size, kind=real64)) * energy_variance
552 
553  if (i_log) then
554  write(logid, "('Local energy = ', F12.8)") energy_mean
555  write(logid, "('Variance = ', F12.8)") energy_variance
556  write(logid, "('Acceptance rate = ', F12.8)") accept_rate
557  write(logid, "('')")
558  end if
559 
560  end subroutine vqmc_qho
561 
562  ! ----------------------------------------------------------------------------------------------------------------
563  ! ROUTINE: Grid_QHO
564  !
565  ! DESCRIPTION:
568  !
572  !
573  ! PARAMETERS:
578  ! ----------------------------------------------------------------------------------------------------------------
579  subroutine grid_qho(alpha_array, energies, acc_rates, variances)
580 
581  implicit none
582 
583  real(real64), dimension(:), intent(in) :: alpha_array
584  real(real64), dimension(:), intent(inout) :: energies
585  real(real64), dimension(:), intent(inout) :: acc_rates
586  real(real64), dimension(:), intent(inout) :: variances
587 
588  real(real64), dimension(:), allocatable :: energy_chain, trimmed_energy_chain
589  real(real64), dimension(:), allocatable :: pos_chain, trimmed_pos_chain
590  real(real64) :: accept_rate
591  real(real64) :: energy_mean
592  real(real64) :: energy_variance
593  integer :: i, ierr
594  character(len=3) :: str_i
595  character(len=40) :: chainfile
596 
597  ! MPI variables
598  real(real64) :: my_ensum
599 
600  ! Set up printed output for the loop, from rank 0.
601  if (my_rank==0) then
602  write(*, "('-------------- Begin VQMC Runs --------------')")
603  write(*, "('Alpha Energy Acceptance Rate')")
604  if (write_log) then
605  write(logid, "('-------------- Begin VQMC Runs --------------')")
606  if (.not. run_equil) then
607  write(logid, "('')")
608  write(logid, "('Logging the chains from Rank 0.')")
609  write(logid, "('')")
610  endif
611  end if
612  endif
613 
614  ! Loop over values of alpha
615  do i = alpha_start, size(alpha_array)
616 
617  ! Calculate energy at current value of alpha
618  call vqmc_qho(qho%sigma, my_steps+qho%burn_step, qho%burn_step, qho%thin_step, &
619  alpha_array(i), pos_chain, accept_rate, energy_chain, energy_mean, &
620  energy_variance)
621 
622  ! Save energy - serial version
623  !energies(i) = energy_mean
624  !variances(i) = energy_variance
625  !acc_rates(i) = accept_rate
626  !write(*, "(F7.4, ' ', F12.8, ' ', F12.8)") alpha_array(i), &
627  ! energy_mean, accept_rate
628 
629  ! Sum energy sums across ranks and save to rank 0
630  my_ensum = sum(energy_chain)
631  call mpi_reduce(my_ensum,energies(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
632 
633  ! Sum accept rates across ranks, to find average in rank 0
634  call mpi_reduce(accept_rate,acc_rates(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
635 
636  ! Sum variances acros ranks, to find average in rank 0
637  call mpi_reduce(energy_variance,variances(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
638 
639  ! Print results from all ranks for checking
640  !write(*, "(F7.4, ' ', F12.8, ' ', F12.8, ' ', F12.8, ' ', I4)") alpha_array(i), &
641  ! energy_mean, accept_rate, energy_variance, my_rank
642 
643  ! Compute and print combined results in rank 0
644  if (my_rank==0) then
645  !print*, 'tot_ensum=',tot_ensum,' ,tot_size=',tot_size
646  energies(i) = energies(i) * per_step
647  acc_rates(i) = acc_rates(i) * per_proc
648  variances(i) = variances(i) * per_proc ! CHECK FORMULA!
649  write(*, "(F7.4, ' ', F12.8, ' ', F12.8)") alpha_array(i), &
650  energies(i), acc_rates(i)
651 
652  ! Write energy chain from rank 0, if required
653  if (write_chains) then
654 
655  ! Setup filename
656  write(str_i, "(I3)") i
657  write(chainfile, "(3A)") './chains_out/alpha_', trim(adjustl(str_i)), '.nc'
658 
659  ! Check chain length
660  if (size(energy_chain) < max_chain_len) then
661  write(*, "(2A)") 'Writing Markov chain to from rank 0 to', chainfile
662  call write_qho_equilibration(pos_chain, energy_chain, chainfile, ierr, &
663  qho, accept_rate, alpha_array(i))
664  else
665  ! Adjust chain length for writing.
666  ! Allocate temporary arrays
667  allocate(trimmed_pos_chain(max_chain_len))
668  allocate(trimmed_energy_chain(max_chain_len))
669  ! Fill in the temporary arrays
670  trimmed_pos_chain(:) = pos_chain(1:max_chain_len)
671  trimmed_energy_chain(:) = energy_chain(1:max_chain_len)
672  ! Write the trimmed temporary arrays
673  write(*, "(2A)") 'Writing trimmed Markov chain to from rank 0 to', chainfile
674  call write_qho_equilibration(trimmed_pos_chain, trimmed_energy_chain, chainfile, ierr, &
675  qho, accept_rate, alpha_array(i))
676  ! Deallocate the temporary arrays
677  deallocate(trimmed_pos_chain)
678  deallocate(trimmed_energy_chain)
679  endif
680  end if
681 
682  ! Write restart files, if required
683  if (write_restart) then
684  if (mod(i,restart_num)==0) then
685  ! Write res_nml.txt
686  call write_res_nml(i+1)
687  ! Store energy, variance & accept_rate arrays
688  call write_qho_main(alpha_array, energies, variances, resfile_etot, ierr, qho, acc_rates)
689  endif
690  endif
691 
692  endif
693 
694  end do
695 
696  if (my_rank==0) then
697 
698  write(*, "('-------------- End VQMC Runs --------------')")
699  if (write_log) then
700  write(logid, "('-------------- End VQMC Runs --------------')")
701  write(logid, "('Results:')")
702  write(logid, "('alpha Energy')")
703  do i = 1, size(alpha_array)
704  write(logid, "(F12.8, ' ', F12.8)") alpha_array(i), energies(i)
705  end do
706  write(logid, "('')")
707  write(logid, "('Min Energy = ', F12.8, ' at c = ', F12.8)") &
708  minval(energies), alpha_array(minloc(energies, 1))
709  write(logid, "('')")
710  end if
711 
712  ! Write main results file.
713  write(*, *) 'Calculations finished, writing results to NetCDF output file ', &
714  outfile
715  if (write_log) then
716  write(logid, *) 'Calculations finished, writing results to NetCDF &
717  &output file ', outfile
718  end if
719  call write_qho_main(alpha_array, energies, variances, outfile, ierr, qho, acc_rates)
720 
721  endif
722 
723  end subroutine grid_qho
724 
725  ! ----------------------------------------------------------------------------------------------------------------
726  ! ROUTINE: Equil_QHO
727  !
728  ! DESCRIPTION:
732  !
737  !
738  ! PARAMETERS:
741  ! ----------------------------------------------------------------------------------------------------------------
742  subroutine equil_qho(alpha_array)
743 
744  implicit none
745 
746  real(real64), dimension(:), intent(in) :: alpha_array
747 
748  real(real64), dimension(:), allocatable :: energy_chain
749  real(real64), dimension(:), allocatable :: pos_chain
750  real(real64) :: energy
751  real(real64) :: accept_rate
752  real(real64) :: variance
753  integer :: i, ierr
754  character(len=3) :: str_i
755  character(len=40) :: equilfile
756 
757  ! Set up printed output for the loop.
758  write(*, "('-------------- Begin VQMC Runs --------------')")
759  write(*, "('Alpha Energy Acceptance Rate')")
760  if (write_log) then
761  write(logid, "('-------------- Begin VQMC Runs --------------')")
762  end if
763 
764  ! Loop over values of alpha
765  do i = alpha_start, size(alpha_array)
766 
767  ! Calculate energy at current value of alpha
768  call vqmc_qho(qho%sigma,qho%steps, qho%burn_step, qho%thin_step, &
769  alpha_array(i), pos_chain, accept_rate, energy_chain, energy, &
770  variance)
771 
772  ! Report the energy
773  write(*, "(F7.4, ' ', F12.8, ' ', F12.8)") alpha_array(i), &
774  energy, accept_rate
775  if (write_log) then
776  write(logid, "('')")
777  write(logid, "('Results:')")
778  write(logid, "('alpha = ', F12.8, ', Energy = ', F12.8)") alpha_array(i), energy
779  write(logid, "('')")
780  end if
781 
782  ! Save the energy trace for each value of alpha.
783  write(str_i, "(I3)") i
784  write(equilfile, "(3A)") './equil_out/alpha_', trim(adjustl(str_i)),'.nc'
785  write(*, *) 'Writing equilibration outfile to ', equilfile
786  call write_qho_equilibration(pos_chain, energy_chain, equilfile, ierr, qho, &
787  accept_rate, alpha_array(i))
788 
789  ! Write restart files, if required
790  if (write_restart) then
791  if (mod(i,restart_num)==0) then
792  ! Write res_nml.txt
793  call write_res_nml(i+1)
794  endif
795  endif
796 
797  end do
798 
799  ! Closing print statements
800 
801  write(*, "('-------------- End VQMC Runs --------------')")
802  if (write_log) then
803  write(logid, "('-------------- End VQMC Runs --------------')")
804  end if
805 
806  !if (write_log) then
807  ! write(logid, "('Results:')")
808  ! write(logid, "('alpha Energy')")
809  ! do i = 1, size(alpha_array)
810  ! write(logid, "(F12.8, ' ', F12.8)") alpha_array(i), energies(i)
811  ! end do
812  ! write(logid, "('')")
813  !end if
814 
815  end subroutine equil_qho
816 
817 ! --------------------------------------------------------------------------------------------------------------------
818 ! SECTION: VQMC for H_2_plus
819 ! --------------------------------------------------------------------------------------------------------------------
820 
821  ! ----------------------------------------------------------------------------------------------------------------
822  ! ROUTINE: VQMC_H_2_plus
823  !
824  ! DESCRIPTION:
827  !
839  !
840  ! PARAMETERS:
849  !
856  ! ----------------------------------------------------------------------------------------------------------------
857  subroutine vqmc_h_2_plus(sigma, N, burn, thin, c_const, R_a, R_b, &
858  X_burnt, accept_rate, energy_chain, energy_mean, energy_variance, &
859  DphiDc_chain, x_start)
860 
861  implicit none
862 
863  real(real64), intent(in), dimension(3) :: R_a, R_b
864  real(real64), intent(in) :: sigma
865  real(real64), intent(in) :: c_const
866  real(real64), dimension(3), optional :: x_start
867  integer, intent(in) :: N, burn, thin
868 
869 
870  real(real64), intent(out) :: accept_rate, energy_mean, energy_variance
871  real(real64), intent(out), dimension(:, :), allocatable :: X_burnt
872  real(real64), intent(out), dimension(:), allocatable :: energy_chain
873  real(real64), intent(out), dimension(:), allocatable :: DphiDc_chain
874 
875  real(real64) :: E_loc ! Local energy at current MCMC walker position
876  real(real64) :: DphiDc ! Gradient wrt c at current MCMC walker position
877  real(real64) :: log_prob_n ! The natural logarithm of the transition proability
878  real(real64) :: a_prob_trans ! Acceptance ratio of the transitional probability
879  real(real64) :: unif_num ! A random sample from the uniform distribution
880  real(real64), dimension(3) :: start ! Starting point of MCMC chain
881  real(real64), dimension(:, :), allocatable :: X ! Unburnt and thinned MCMC chain
882  real(real64), dimension(3) :: X_n ! The proposed step
883  real(real64), dimension(3) :: sample_normal ! A sample from the normal distribution
884  integer :: accept ! Number of accepted steps
885  integer :: thin_counter ! Counter to keep track of the thinning process
886  integer :: reset_counter ! Counts how many times thin counter was reset
887  integer :: re_size ! The size of the burnt and thinned array
888  integer :: i_bt ! Burned/thinned value of i
889  integer :: i ! Loop variable
890  character(len=1) :: acceptance ! Yes/No signpost for acceptance
891 
892  if (i_log) then
893  write(logid, "('Markov chain for c = ', F7.4)") c_const
894  if (.not. run_equil) write(logid, "('Logging the chain from Rank 0.')")
895  endif
896 
897  ! Adjust the output array to be of the burnt and thinned array
898  re_size = (n - burn) / thin
899 
900  allocate(x(3, 1:n+1)); x(:, :) = 0.0_real64
901  allocate(x_burnt(3, 1:re_size)); x_burnt(:, :) = 0.0_real64
902  allocate(energy_chain(1:re_size)); energy_chain(:) = 0.0_real64
903  allocate(dphidc_chain(1:re_size)); dphidc_chain(:) = 0.0_real64
904 
905  ! Allow the starting position to take on the default value
906  if (present(x_start)) then
907  start = x_start
908  x(:,1) = start
909  if (i_log) then
910  write(logid, "('Starting from provided Markov chain position:')")
911  write(logid, "(' x = ', F12.8, ' y = ', F12.8, ' z = ', F12.8)") &
912  x(1, 1), x(2, 1), x(3, 1)
913  end if
914  else
915  ! Create the initial position of the random walker that is within a certain radius
916  call random_normal(3, sample_normal)
917  do while(norm2(sample_normal) >= h2plus%rcut)
918  call random_normal(3, sample_normal)
919  end do
920  start = sample_normal
921  x(:,1) = start
922  if (i_log) then
923  write(logid, "('Generated random Markov chain start position:')")
924  write(logid, "('X = ', F12.8, ', Y = ', F12.8, ', Z = ', F12.8)") &
925  x(1, 1), x(2, 1), x(3, 1)
926  end if
927  end if
928 
929  ! Start the counter for the number of accepted moves
930  accept = 0
931  ! Start the counter for the thinning process
932  thin_counter = 0
933  ! Start the reset_counter
934  reset_counter = 0
935 
936  if (i_log) then
937  write(logid, "('Sample X Y Z &
938  &Local Energy Accepted?')")
939  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8, ' ', F12.6)")&
940  1, x(1, 1), x(2, 1), x(3, 1), energy_chain(1)
941  end if
942 
943  ! Calculate the local energy of the initial step if not burnt
944  if (burn == 0) then
945  call h_2_plus_energy(c_const, x(:, 1), r_a, r_b, e_loc, dphidc)
946  energy_chain(1) = e_loc
947  dphidc_chain(1) = dphidc
948  x_burnt(:, 1) = x(:, 1)
949 
950  do i = 2, n
951  ! Add one to the thin counter
952  thin_counter = 1 + thin_counter
953  ! Propose the next step
954  call random_normal(3, sample_normal)
955  x_n = x(:, i-1) + sigma * sample_normal
956 
957  ! Calculate the transition probability between proposed step and current step
958  log_prob_n = log(h_2_plus_prob(c_const, x(:, i-1), x_n, r_a, r_b))
959 
960  ! Calculate the acceptance ratio
961  a_prob_trans = min(1.0_real64, exp(log_prob_n))
962 
963  ! Generate sample from uniform distribution
964  call ran1(unif_num, seed)
965 
966  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
967  if (a_prob_trans >= unif_num) then
968  ! If accepted add proposed position to the MCMC chain
969  x(:, i) = x_n
970 
971  ! To determine whether this point should be added to the burnt and thinned MCMC chain
972  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
973  ! Increment the reset counter
974  reset_counter = 1 + reset_counter
975  ! Precalculate i_bt for shorthand
976  i_bt = i - burn - (reset_counter * (thin - 1))
977  ! Add to the burn and thinned MCMC chain
978  x_burnt(:, i_bt) = x_n
979  ! Calculate the local energy at this point
980  call h_2_plus_energy(c_const, x_n, r_a, r_b, e_loc, dphidc)
981  energy_chain(i_bt) = e_loc
982  dphidc_chain(i_bt) = dphidc
983  ! Resetting the thin counter
984  thin_counter = 0
985 
986  acceptance = 'Y'
987  ! Increment the number of accepted steps
988  accept = 1 + accept
989 
990  if (i_log) then
991  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
992  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
993  &' ', F12.6, ' ', A1)") (i_bt), x_burnt(1, i_bt), &
994  x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
995  acceptance
996  end if
997  end if
998  end if
999  else
1000  ! If rejected the proposed step is the current step
1001  x(:, i) = x(:, i-1)
1002 
1003  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1004  ! Increment the reset counter
1005  reset_counter = 1 + reset_counter
1006  ! Precalculate i_bt for shorthand
1007  i_bt = i - burn - (reset_counter * (thin - 1))
1008  ! Add to the burn and thinned MCMC chain
1009  x_burnt(:, i_bt) = x(:, i-1)
1010  ! Calculate the local energy at this point
1011  call h_2_plus_energy(c_const, x(:, i-1), r_a, r_b, e_loc, dphidc)
1012  energy_chain(i_bt) = e_loc
1013  dphidc_chain(i_bt) = dphidc
1014  ! Resetting the thin counter
1015  thin_counter = 0
1016 
1017  acceptance = 'N'
1018 
1019  if (i_log) then
1020  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
1021  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
1022  &' ', F12.6, ' ', A1)") (i_bt), x_burnt(1, i_bt), &
1023  x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
1024  acceptance
1025  end if
1026  end if
1027  end if
1028  end if
1029  end do
1030  ! Loop to advance the chain to the final burn step and not save any of the outputs
1031  else if (burn >= 2) then
1032  do i = 2, burn
1033  ! Propose the next step
1034  call random_normal(3, sample_normal)
1035  x_n = x(:, i-1) + sigma * sample_normal
1036 
1037  ! Calculate the transition probability between proposed step and current step
1038  log_prob_n = log(h_2_plus_prob(c_const, x(:, i-1), x_n, r_a, r_b))
1039 
1040  ! Calculate the acceptance ratio
1041  a_prob_trans = min(1.0_real64, exp(log_prob_n))
1042 
1043  ! Generate sample from uniform distribution
1044  call ran1(unif_num, seed)
1045 
1046  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
1047  if (a_prob_trans >= unif_num) then
1048  ! If accepted add proposed position to the MCMC chain
1049  x(:, i) = x_n
1050  else
1051  ! If rejected the proposed step is the current step
1052  x(:, i) = x(:, i-1)
1053  end if
1054  end do
1055  end if
1056 
1057  ! Metropolis Algorithm to calculate the burnt and thinned MCMC chain
1058  if (burn >= 1) then
1059  ! Metropolis algorithm
1060  do i = burn + 1, n+1
1061  ! Add one to the thin counter
1062  thin_counter = 1 + thin_counter
1063  ! Propose the next step
1064  call random_normal(3, sample_normal)
1065  x_n = x(:, i-1) + sigma * sample_normal
1066 
1067  ! Calculate the transition probability between proposed step and current step
1068  log_prob_n = log(h_2_plus_prob(c_const, x(:, i-1), x_n, r_a, r_b))
1069 
1070  ! Calculate the acceptance ratio
1071  a_prob_trans = min(1.0_real64, exp(log_prob_n))
1072 
1073  ! Generate sample from uniform distribution
1074  call ran1(unif_num, seed)
1075 
1076  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
1077  if (a_prob_trans >= unif_num) then
1078  ! If accepted add proposed position to the MCMC chain
1079  x(:, i) = x_n
1080 
1081  ! To determine whether this point should be added to the burnt and thinned MCMC chain
1082  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1083  ! Increment the reset counter
1084  reset_counter = 1 + reset_counter
1085  ! Precalculate i_bt for shorthand
1086  i_bt = i - burn - (reset_counter * (thin - 1))
1087  ! Add to the burn and thinned MCMC chain
1088  x_burnt(:, i_bt) = x_n
1089  ! Calculate the local energy at this point
1090  call h_2_plus_energy(c_const, x_n, r_a, r_b, e_loc, dphidc)
1091  energy_chain(i_bt) = e_loc
1092  dphidc_chain(i_bt) = dphidc
1093  ! Resetting the thin counter
1094  thin_counter = 0
1095 
1096  acceptance = 'Y'
1097  ! Increment the number of accepted steps
1098  accept = 1 + accept
1099 
1100  if (i_log) then
1101  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
1102  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
1103  &' ', F12.6, ' ', A1)") (i_bt), x_burnt(1, i_bt), &
1104  x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
1105  acceptance
1106  end if
1107  end if
1108  end if
1109 
1110  else
1111  ! If rejected the proposed step is the current step
1112  x(:, i) = x(:, i-1)
1113 
1114  ! To determine whether this point should be added to the burnt and thinned MCMC chain
1115  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1116  ! Increment the reset counter
1117  reset_counter = 1 + reset_counter
1118  ! Precalculate i_bt for shorthand
1119  i_bt = i - burn - (reset_counter * (thin - 1))
1120  ! Add to the burn and thinned MCMC chain
1121  x_burnt(:, i_bt) = x(:, i-1)
1122  ! Calculate the local energy at this point
1123  call h_2_plus_energy(c_const, x(:, i-1), r_a, r_b, e_loc, dphidc)
1124  energy_chain(i_bt) = e_loc
1125  dphidc_chain(i_bt) = dphidc
1126  ! Resetting the thin counter
1127  thin_counter = 0
1128 
1129  acceptance = 'N'
1130 
1131  if (i_log) then
1132  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
1133  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
1134  &' ', F12.6, ' ', A1)") (i_bt), x_burnt(1, i_bt), &
1135  x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
1136  acceptance
1137  end if
1138  end if
1139 
1140  end if
1141  end if
1142  end do
1143  end if
1144 
1145  ! Calculate the acceptance ratio
1146  accept_rate = (real(accept, kind=real64) / real(re_size, kind=real64))
1147 
1148  ! Calculate the mean energy and the variance in the energy
1149  energy_mean = sum(energy_chain) / real(re_size, kind=real64)
1150  energy_variance = 0.0_real64
1151  do i = 1, re_size
1152  energy_variance = energy_variance + ((1.0_real64 / real(re_size - 1, kind=real64)) * &
1153  (energy_chain(i) - energy_mean) * (energy_chain(i) - energy_mean))
1154  end do
1155  energy_variance = (1.0_real64 / real(re_size, kind=real64)) * energy_variance
1156 
1157  if (i_log) then
1158  write(logid, "('Local energy = ', F12.8)") energy_mean
1159  write(logid, "('Variance = ', F12.8)") energy_variance
1160  write(logid, "('Acceptance rate = ', F12.8)") accept_rate
1161  write(logid, "('')")
1162  end if
1163 
1164  end subroutine vqmc_h_2_plus
1165 
1166  ! ----------------------------------------------------------------------------------------------------------------
1167  ! ROUTINE: Auto_H_2_plus
1168  !
1169  ! DESCRIPTION:
1173  !
1178  !
1179  ! PARAMETERS:
1187  ! ----------------------------------------------------------------------------------------------------------------
1188  subroutine auto_h_2_plus(c_in, bond_length, c_out, energy_out, variance_out, &
1189  accept_out, converged)
1190 
1191  implicit none
1192 
1193  real(real64), intent(in) :: c_in
1194  real(real64), intent(in) :: bond_length
1195  real(real64), intent(out) :: c_out
1196  real(real64), intent(out) :: energy_out
1197  real(real64), intent(out) :: variance_out
1198  real(real64), intent(out) :: accept_out
1199  character(len=1), intent(out) :: converged
1200 
1201  real(real64), dimension(3) :: R_a, R_b
1202  real(real64) :: c_current, c_new
1203  real(real64) :: c_grad
1204  real(real64), dimension(:, :), allocatable :: pos_chain
1205  real(real64), dimension(:), allocatable :: energy_chain
1206  real(real64), dimension(:), allocatable :: DphiDc_chain
1207  integer :: viz_steps
1208  integer :: k
1209  integer :: ierr
1210  integer :: maxiter = 200
1211  real(real64) :: tol = 1e-4_real64
1212  character(len=6) :: str_b
1213  character(len=40) :: chainfile
1214 
1215  ! MPI Variables
1216  real(real64) :: my_c_grad, my_ensum
1217  real(real64) :: sum_var, sum_acc
1218 
1219  ! Initialise
1220  c_current = c_in
1221  r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
1222  r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
1223 
1224  if (i_log) then
1225  write(logid, "('Starting parameter search for c...')")
1226  end if
1227 
1228  ! Main iteration loop
1229  do k = auto_start, maxiter
1230 
1231  ! Calculate gradient at current c
1232  call vqmc_h_2_plus(h2plus%sigma, my_steps+h2plus%burn_step, h2plus%burn_step, h2plus%thin_step, &
1233  c_current, r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, &
1234  variance_out, dphidc_chain)
1235 
1236  ! Find gradient
1237  my_c_grad = 2.0_real64 * (sum(energy_chain * dphidc_chain) / size(energy_chain)) &
1238  - (energy_out * (sum(dphidc_chain) / size(dphidc_chain)))
1239 
1240  ! Find average c_grad across ranks, on all ranks.
1241  call mpi_allreduce(my_c_grad,c_grad,1,mpi_double_precision,mpi_sum,mpi_comm_world,mpierr)
1242  c_grad = c_grad * per_proc
1243 
1244  ! Update the log
1245  if (i_log) then
1246  write(logid, "('c Gradient')")
1247  write(logid, "(F13.7, ' ', F13.7)") c_current, c_grad
1248  end if
1249 
1250  ! Check if gradient is within tolerance.
1251  if (abs(c_grad) .le. tol) then
1252 
1253  ! Update c_out, and note convergence, on all ranks.
1254  c_out = c_current
1255  converged = 'Y'
1256 
1257  ! Sum the chain energy-sums and store in rank 0
1258  my_ensum = sum(energy_chain)
1259  call mpi_reduce(my_ensum,energy_out,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
1260 
1261  ! Sum the chain variances and store in rank 0
1262  call mpi_reduce(variance_out, sum_var,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
1263 
1264  ! Sum the chain acceptance rates and store in rank 0
1265  call mpi_reduce(accept_out, sum_acc,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
1266 
1267  ! Calculate averages on rank 0
1268  if (my_rank==0) then
1269 
1270  energy_out = energy_out * per_step
1271  variance_out = sum_var * per_proc2
1272  accept_out = sum_acc * per_proc
1273 
1274  ! Update the log
1275  if (write_log) then
1276  write(logid, "('Parameter search converged!')")
1277  write(logid, "('Results:')")
1278  write(logid, "('c = ', F13.7, ', Energy = ', F12.8, ', Variance = ', F12.8)") &
1279  c_current, energy_out, variance_out
1280  write(logid, "('')")
1281  end if
1282 
1283  end if
1284 
1285  ! Exit subroutine, as convergence is reached.
1286  ! Note: Only rank 0 has the correct array outputs. All ranks have correct c_out and 'converged'.
1287  exit
1288 
1289  end if
1290 
1291  ! If not, update c by steepest descent, on all ranks.
1292  call h_2_plus_update_c(c_current, c_grad, c_new)
1293  c_current = c_new
1294 
1295  ! Write restart files, if required
1296  if (write_restart .and. my_rank==0) then
1297  if (mod(k,restart_num)==0) then
1298  call write_res_nml(i_bond,k+1,h2plus%auto_params,c_current)
1299  endif
1300  endif
1301 
1302  end do
1303 
1304  ! Deal with non-convergence
1305  if (k .ge. maxiter) then
1306  converged = 'N'
1307  c_out = c_current
1308  if (i_log) then
1309  write(logid, "('Parameter search NOT converged! Consider changing&
1310  & the starting value of c, or the damping in the H_2_update_c&
1311  & function.')")
1312  end if
1313  end if
1314 
1315  ! Report results
1316  if (my_rank==0) then
1317 
1318  write(*, "(F12.8, ' ', A, ' ', F13.7, ' ', F12.8)") &
1319  bond_length, converged, c_out, energy_out
1320 
1321  ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
1322  if (write_chains) then
1323  write(str_b, "(F6.3)") bond_length
1324  write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
1325  write(*, "('Running VQMC at c = ', F12.7, ' to produce single Markov chain...')") c_out
1326  viz_steps = min(h2plus%steps,max_chain_len)
1327  call vqmc_h_2_plus(h2plus%sigma, viz_steps, h2plus%burn_step, h2plus%thin_step, &
1328  c_out, r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, &
1329  variance_out, dphidc_chain)
1330  write(*, "(2A)") 'Writing Markov chain to ', chainfile
1331  call write_h2plus_equilibration(pos_chain, energy_chain, chainfile, ierr, h2plus, &
1332  accept_out, c_out, bond_length)
1333  end if
1334  end if
1335 
1336  end subroutine auto_h_2_plus
1337 
1338  ! ----------------------------------------------------------------------------------------------------------------
1339  ! ROUTINE: Grid_H_2_plus
1340  !
1341  ! DESCRIPTION:
1345  !
1350  !
1351  ! PARAMETERS:
1358  ! ----------------------------------------------------------------------------------------------------------------
1359  subroutine grid_h_2_plus(c_array, bond_length, c_out, energy_out, variance_out, &
1360  accept_out)
1361 
1362  implicit none
1363 
1364  real(real64), dimension(:), intent(in) :: c_array
1365  real(real64), intent(in) :: bond_length
1366  real(real64), intent(out) :: c_out
1367  real(real64), intent(out) :: energy_out
1368  real(real64), intent(out) :: variance_out
1369  real(real64), intent(out) :: accept_out
1370 
1371  real(real64), dimension(3) :: R_a, R_b
1372  real(real64), dimension(:, :), allocatable :: pos_chain
1373  real(real64), dimension(:), allocatable :: energy_chain, DphiDc_chain
1374  real(real64), dimension(:), allocatable :: param_energies, param_accepts, param_variances
1375  integer :: viz_steps
1376  integer :: i
1377  integer :: ierr
1378  character(len=6) :: str_b
1379  character(len=40) :: chainfile
1380 
1381  ! MPI variables
1382  real(real64) :: my_ensum
1383 
1384  ! Initialise
1385  r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
1386  r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
1387  allocate(param_energies(size(c_array))); param_energies(:)= 0.0_real64
1388  allocate(param_accepts(size(c_array))); param_accepts(:) = 0.0_real64
1389  allocate(param_variances(size(c_array))); param_variances(:) = 0.0_real64
1390 
1391  ! Sort out restarting
1392  if (run_restart) then
1393  ! Read in array values
1394  call read_energies(resfile_epar, param_energies, param_variances, param_accepts)
1395  ! Reset the tag
1396  run_restart = .false.
1397  endif
1398 
1399  if (i_log) then
1400  write(logid, "('Starting parameter search for c...')")
1401  end if
1402 
1403  ! Loop over c values
1404  do i = grid_start, size(c_array)
1405 
1406  ! Calculate energy at current value of c
1407  call vqmc_h_2_plus(h2plus%sigma, my_steps+h2plus%burn_step, h2plus%burn_step, h2plus%thin_step, &
1408  c_array(i), r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, &
1409  variance_out, dphidc_chain)
1410 
1411  ! Save energy - serial version
1412  !param_energies(i) = energy
1413  !param_variances(i) = variance
1414  !param_accepts(i) = accept_rate
1415 
1416  ! Sum the chain energy-sums and store in rank 0
1417  my_ensum = sum(energy_chain)
1418  call mpi_reduce(my_ensum,param_energies(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
1419 
1420  ! Sum the chain variances and store in rank 0
1421  call mpi_reduce(variance_out,param_variances(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
1422 
1423  ! Sum the chain acceptance rates and store in rank 0
1424  call mpi_reduce(accept_out,param_accepts(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
1425 
1426  ! Calculate averages on rank 0
1427  if (my_rank==0) then
1428  param_energies(i) = param_energies(i) * per_step
1429  param_variances(i) = param_variances(i) * per_proc2
1430  param_accepts(i) = param_accepts(i) * per_proc
1431 
1432  ! Write restart files, if required
1433  if (write_restart) then
1434  if (mod(i,restart_num)==0) then
1435  ! Write res_nml.txt
1436  call write_res_nml(i_bond,i+1,h2plus%auto_params)
1437  ! Write param energies, variances and acceptance rates
1438  call write_qho_main(c_array, param_energies, param_variances, resfile_epar, ierr, qho, param_accepts)
1439  endif
1440  endif
1441  endif
1442 
1443  end do
1444 
1445  ! Compile final results of grid search on rank 0
1446  ! Note: Only rank 0 has the correct _out values.
1447  if (my_rank==0) then
1448 
1449  ! Select lowest energy and corresponding c value
1450  energy_out = minval(param_energies)
1451  c_out = c_array(minloc(param_energies, 1))
1452  variance_out = param_variances(minloc(param_energies, 1))
1453  accept_out = param_accepts(minloc(param_energies, 1))
1454 
1455  write(*, "(F12.8, ' - ', F13.7, ' ', F12.8)") &
1456  bond_length, c_out, energy_out
1457 
1458  if (write_log) then
1459  write(logid, "('Results:')")
1460  write(logid, "('c Energy')")
1461  do i = 1, size(c_array)
1462  write(logid, "(F12.8, ' ', F12.8)") c_array(i), param_energies(i)
1463  end do
1464  write(logid, "('')")
1465  write(logid, "('Min Energy = ', F12.8, ' at c = ', F12.8)") energy_out, c_out
1466  write(logid, "('')")
1467  end if
1468 
1469  ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
1470  if (write_chains) then
1471  write(str_b, "(F6.3)") bond_length
1472  write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
1473  write(*, "('Running VQMC at c = ', F12.7, ' to produce single Markov chain...')") c_out
1474  viz_steps = min(h2plus%steps,max_chain_len)
1475  call vqmc_h_2_plus(h2plus%sigma, viz_steps, h2plus%burn_step, h2plus%thin_step, &
1476  c_out, r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, &
1477  variance_out, dphidc_chain)
1478  write(*, "(2A)") 'Writing Markov chain to ', chainfile
1479  call write_h2plus_equilibration(pos_chain, energy_chain, chainfile, ierr, h2plus, &
1480  accept_out, c_out, bond_length)
1481  end if
1482 
1483  endif
1484 
1485  end subroutine grid_h_2_plus
1486 
1487  ! ----------------------------------------------------------------------------------------------------------------
1488  ! ROUTINE: Equil_H_2_plus
1489  !
1490  ! DESCRIPTION:
1494  !
1499  !
1500  ! PARAMETERS:
1504  ! ----------------------------------------------------------------------------------------------------------------
1505  subroutine equil_h_2_plus(c_array, bond_length, bonds_pos)
1506 
1507  implicit none
1508 
1509  real(real64), dimension(:), intent(in) :: c_array
1510  real(real64), intent(in) :: bond_length
1511  integer, intent(in) :: bonds_pos
1512 
1513  real(real64), dimension(3) :: R_a, R_b
1514  real(real64), dimension(:), allocatable :: variances
1515  real(real64), dimension(:), allocatable :: energies
1516  real(real64), dimension(:), allocatable :: energy_chain
1517  real(real64), dimension(:, :), allocatable :: pos_chain
1518  real(real64), dimension(:), allocatable :: DphiDc_chain
1519  real(real64) :: accept_rate
1520  real(real64) :: energy_variance
1521  integer :: i, ierr
1522  character(len=3) :: str_i, str_b
1523  character(len=40) :: equilfile
1524 
1525  ! Initialise
1526  r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
1527  r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
1528  allocate(energies(size(c_array))); energies(:) = 0.0_real64
1529  allocate(variances(size(c_array)))
1530 
1531  ! Loop over values of c
1532  do i = grid_start, size(c_array)
1533 
1534  ! Calculate energy at current value of c
1535  call vqmc_h_2_plus(h2plus%sigma, h2plus%steps, h2plus%burn_step, h2plus%thin_step, &
1536  c_array(i), r_a, r_b, pos_chain, accept_rate, energy_chain,energies(i), &
1537  energy_variance, dphidc_chain)
1538 
1539  write(*, "(F12.8, ' - ', F13.7, ' ', F12.8)") &
1540  bond_length, c_array(i), energies(i)
1541 
1542  ! Save the energy trace for each value of c
1543  write(str_i, "(I3)") i
1544  write(str_b, "(I3)") bonds_pos
1545  write(equilfile, "(5A)") './equil_out/b', trim(adjustl(str_b)), &
1546  '_c', trim(adjustl(str_i)), '.nc'
1547  write(*, *) 'Writing equilibration outfile to ', equilfile
1548  call write_h2plus_equilibration(pos_chain, energy_chain, equilfile, ierr, h2plus, &
1549  accept_rate, c_array(i), bond_length)
1550 
1551  ! Write restart files, if required
1552  if (write_restart) then
1553  if (mod(i,restart_num)==0) then
1554  call write_res_nml(i_bond,i+1,h2%auto_params)
1555  endif
1556  endif
1557 
1558  end do
1559 
1560  if (write_log) then
1561  write(logid, "('Results:')")
1562  write(logid, "('c Energy')")
1563  do i = 1, size(c_array)
1564  write(logid, "(F12.8, ' ', F12.8)") c_array(i), energies(i)
1565  end do
1566  write(logid, "('')")
1567  end if
1568 
1569  deallocate(energies)
1570  deallocate(variances)
1571 
1572  end subroutine equil_h_2_plus
1573 
1574 ! --------------------------------------------------------------------------------------------------------------------
1575 ! SECTION: VQMC for H_2
1576 ! --------------------------------------------------------------------------------------------------------------------
1577 
1578  ! ----------------------------------------------------------------------------------------------------------------
1579  ! ROUTINE: VQMC_H_2
1580  !
1581  ! DESCRIPTION:
1584  !
1596  !
1597  ! PARAMETERS:
1608  !
1616  ! ----------------------------------------------------------------------------------------------------------------
1617  subroutine vqmc_h_2(sigma, N, burn, thin, a, beta, R_a, R_b, &
1618  X_1_burnt, X_2_burnt, accept_rate, energy_chain, energy_mean, energy_variance, &
1619  DphiDbeta_chain, x_start_1, x_start_2)
1620 
1621  implicit none
1622 
1623  real(real64), intent(in), dimension(3) :: R_a, R_b
1624  real(real64), intent(in) :: sigma
1625  real(real64), intent(in) :: a, beta
1626  real(real64), dimension(3), optional :: x_start_1, x_start_2
1627  integer, intent(in) :: N, burn, thin
1628 
1629 
1630  real(real64), intent(out) :: accept_rate, energy_mean, energy_variance
1631  real(real64), intent(out), dimension(:, :), allocatable :: X_1_burnt, X_2_burnt
1632  real(real64), intent(out), dimension(:), allocatable :: energy_chain
1633  real(real64), intent(out), dimension(:), allocatable :: DphiDbeta_chain
1634 
1635  real(real64) :: E_loc ! Local energy at current MCMC walker position
1636  real(real64) :: DphiDbeta ! Gradient wrt beta at current MCMC walker position
1637  real(real64) :: log_prob_n ! The natural logarithm of the transition proability
1638  real(real64) :: a_prob_trans ! Acceptance ratio of the transitional probability
1639  real(real64) :: unif_num ! A random sample from the uniform distribution
1640  real(real64), dimension(3) :: start_1, start_2 ! Starting points for the random walkers
1641  real(real64), dimension(3) :: X_n_1, X_n_2 ! The propsed steps for the random walkers
1642  real(real64), dimension(3) :: sample_normal ! Sample from the normal distribution
1643  real(real64), dimension(:, :), allocatable :: X_1, X_2 ! Full MCMC chains
1644  integer :: accept ! Number of accepted steps
1645  integer :: thin_counter ! Counter to keep track of the thinning process
1646  integer :: reset_counter ! Counts how many times thin counter was reset
1647  integer :: re_size ! The size of the burnt and thinned array
1648  integer :: i_bt ! Burned/thinned value of i
1649  integer :: i ! Loop variable
1650  character(len=1) :: acceptance ! Yes/No signpost for acceptance
1651 
1652  if (i_log) then
1653  write(logid, "('Markov chain for beta = ', F7.4)") beta
1654  if (.not. run_equil) write(logid, "('Logging the chain from Rank 0.')")
1655  endif
1656 
1657  ! Adjust the output array to be of the burnt and thinned array
1658  re_size = (n - burn) / thin
1659 
1660  allocate(x_1(3, 1:n+1)); x_1(:, :) = 0.0_real64
1661  allocate(x_2(3, 1:n+1)); x_2(:, :) = 0.0_real64
1662  allocate(x_1_burnt(3, 1:re_size)); x_1_burnt(:, :) = 0.0_real64
1663  allocate(x_2_burnt(3, 1:re_size)); x_2_burnt(:, :) = 0.0_real64
1664  allocate(energy_chain(1:re_size)); energy_chain(:) = 0.0_real64
1665  allocate(dphidbeta_chain(1:re_size)); dphidbeta_chain(:) = 0.0_real64
1666 
1667  ! Allow the starting position to take on the default value
1668  if (present(x_start_1) .and. present(x_start_2)) then
1669  start_1 = x_start_1
1670  x_1(:,1) = start_1
1671 
1672  start_2 = x_start_2
1673  x_2(:,1) = start_2
1674 
1675  if (i_log) then
1676  write(logid, "('Starting from provided Markov chain position:')")
1677  write(logid, "(' x1 = ', F12.8, ' y1 = ', F12.8, ' z1 = ', F12.8)") &
1678  x_1(1, 1), x_1(2, 1), x_1(3, 1)
1679  write(logid, "(' x2 = ', F12.8, ' y2 = ', F12.8, ' z2 = ', F12.8)") &
1680  x_2(1, 1), x_2(2, 1), x_2(3, 1)
1681  end if
1682  else
1683  ! Create the initial position of random walker 1 that is within a certain radius
1684  call random_normal(3, sample_normal)
1685  do while(norm2(sample_normal) >= h2%rcut)
1686  call random_normal(3, sample_normal)
1687  end do
1688  start_1 = sample_normal
1689  x_1(:,1) = start_1
1690 
1691  ! Create the initial position of random walker 2 that is within a certain radius
1692  call random_normal(3, sample_normal)
1693  do while(norm2(sample_normal) >= h2%rcut)
1694  call random_normal(3, sample_normal)
1695  end do
1696  start_2 = sample_normal
1697  x_2(:,1) = start_2
1698  end if
1699 
1700  ! Start the counter for the number of accepted moves
1701  accept = 0
1702  ! Start the counter for the thinning process
1703  thin_counter = 0
1704  ! Start the reset counter
1705  reset_counter = 0
1706 
1707  if (i_log) then
1708  write(logid, "('Sample X1 Y1 Z1 &
1709  &X2 Y2 Z2 &
1710  &Local Energy Accepted?')")
1711  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8, ' ', F12.8, &
1712  &' ', F12.8, ' ', F12.8, ' ', F12.6)") 1, x_1(1, 1), x_1(2, 1), x_1(3, 1),&
1713  x_2(1, 1), x_2(2, 1), x_2(3, 1), energy_chain(1)
1714  end if
1715 
1716  ! Calculate the local energy of the initial step if not burnt
1717  if (burn == 0) then
1718  call h_2_energy(a, beta, x_1(:,1), x_2(:,1), r_a, r_b, e_loc, dphidbeta)
1719  energy_chain(1) = e_loc
1720  dphidbeta_chain(1) = dphidbeta
1721  x_1_burnt(:, 1) = x_1(:, 1)
1722  x_2_burnt(:, 1) = x_2(:, 1)
1723 
1724  do i = 2, n
1725  ! Initialize thin counter
1726  thin_counter = 1 + thin_counter
1727 
1728  ! Propose the next step for random walker 1
1729  call random_normal(3, sample_normal)
1730  x_n_1 = x_1(:, i-1) + sigma * sample_normal
1731 
1732  ! Propose the next step for random walker 2
1733  call random_normal(3, sample_normal)
1734  x_n_2 = x_2(:, i-1) + sigma * sample_normal
1735 
1736  ! Calculate the transition probability between the proposed steps and current steps
1737  log_prob_n = log(h_2_prob(a, beta, x_1(:, i-1), x_2(:, i-1), x_n_1, x_n_2, r_a, r_b))
1738 
1739  ! Calculate the acceptance ratio
1740  a_prob_trans = min(1.0_real64, exp(log_prob_n))
1741 
1742  ! Generate sample from uniform distribution
1743  call ran1(unif_num, seed)
1744 
1745  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
1746  if (a_prob_trans >= unif_num) then
1747  ! If accepted calculate local energy and add proposed position to the MCMC chain
1748  x_1(:, i) = x_n_1
1749  x_2(:, i) = x_n_2
1750 
1751  ! To determine whether this point should be added to the burnt and thinned MCMC chain
1752  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1753  ! Increment the reset counter
1754  reset_counter = 1 + reset_counter
1755  ! Precalculate i_bt for shorthand
1756  i_bt = i - burn - (reset_counter * (thin - 1))
1757  ! Add to the burn and thinned MCMC chain
1758  x_1_burnt(:, i_bt) = x_n_1
1759  x_2_burnt(:, i_bt) = x_n_2
1760  ! Calculate the local energy at this point
1761  call h_2_energy(a, beta, x_n_1, x_n_2, r_a, r_b, e_loc, dphidbeta)
1762  energy_chain(i_bt) = e_loc
1763  dphidbeta_chain(i_bt) = dphidbeta
1764  ! Resetting the thin counter
1765  thin_counter = 0
1766 
1767  ! Increment the number of accepted steps
1768  accept = 1 + accept
1769  acceptance = 'Y'
1770 
1771  if (i_log) then
1772  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
1773  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
1774  & ' ', F12.8, ' ', F12.8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
1775  (i-burn), x_1_burnt(1, i-burn), x_1_burnt(2, i-burn), x_1_burnt(3, i-burn),&
1776  x_2_burnt(1, i-burn), x_2_burnt(2, i-burn), x_2_burnt(3, i-burn),&
1777  energy_chain(i-burn), acceptance
1778  end if
1779  end if
1780  end if
1781  else
1782  ! If rejected the proposed step is the current step
1783  x_1(:, i) = x_1(:, i-1)
1784  x_2(:, i) = x_2(:, i-1)
1785 
1786  ! To determine whether this point should be added to the burnt and thinned MCMC chain
1787  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1788  ! Increment the reset counter
1789  reset_counter = 1 + reset_counter
1790  ! Precalculate i_bt for shorthand
1791  i_bt = i - burn - (reset_counter * (thin - 1))
1792  ! Add to the burn and thinned MCMC chain
1793  x_1_burnt(:, i_bt) = x_1(:, i-1)
1794  x_2_burnt(:, i_bt) = x_2(:, i-1)
1795  ! Calculate the local energy at this point
1796  call h_2_energy(a, beta, x_1(:, i-1), x_2(:, i-1), r_a, r_b, e_loc, dphidbeta)
1797  energy_chain(i_bt) = e_loc
1798  dphidbeta_chain(i_bt) = dphidbeta
1799  ! Resetting the thin counter
1800  thin_counter = 0
1801 
1802  acceptance = 'N'
1803 
1804  if (i_log) then
1805  if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
1806  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
1807  & ' ', F12.8, ' ', F12.8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
1808  (i_bt), x_1_burnt(1, i_bt), x_1_burnt(2, i_bt), x_1_burnt(3, i_bt),&
1809  x_2_burnt(1, i_bt), x_2_burnt(2, i_bt), x_2_burnt(3, i_bt),&
1810  energy_chain(i_bt), acceptance
1811  end if
1812  end if
1813  end if
1814  end if
1815  end do
1816  ! Loop to advance the chain to the final burn step and not save any of the outputs
1817  else if (burn >= 2) then
1818  do i = 2, burn
1819  ! Propose the next step for random walker 1
1820  call random_normal(3, sample_normal)
1821  x_n_1 = x_1(:, i-1) + sigma * sample_normal
1822 
1823  ! Propose the next step for random walker 2
1824  call random_normal(3, sample_normal)
1825  x_n_2 = x_2(:, i-1) + sigma * sample_normal
1826 
1827  ! Calculate the transition probability between the proposed steps and current steps
1828  log_prob_n = log(h_2_prob(a, beta, x_1(:, i-1), x_2(:, i-1), x_n_1, x_n_2, r_a, r_b))
1829 
1830  ! Calculate the acceptance ratio
1831  a_prob_trans = min(1.0_real64, exp(log_prob_n))
1832 
1833  ! Generate sample from uniform distribution
1834  call ran1(unif_num, seed)
1835 
1836  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
1837  if (a_prob_trans >= unif_num) then
1838  ! If accepted calculate local energy and add proposed position to the MCMC chain
1839  x_1(:, i) = x_n_1
1840  x_2(:, i) = x_n_2
1841  else
1842  ! If rejected the proposed step is the current step
1843  x_1(:, i) = x_1(:, i-1)
1844  x_2(:, i) = x_2(:, i-1)
1845  end if
1846  end do
1847 
1848  if (i_log) write(logid, "('---------- Burning steps completed. ----------')")
1849 
1850  end if
1851 
1852  ! Metropolis Algorithm to calculate the burnt and thinned MCMC chain
1853  if (burn >= 1) then
1854  ! Metropolis algorithm
1855  do i = (burn + 1) , n+1
1856  ! Initialize thin counter
1857  thin_counter = 1 + thin_counter
1858  ! Propose the next step for random walker 1
1859  call random_normal(3, sample_normal)
1860  x_n_1 = x_1(:, i-1) + sigma * sample_normal
1861 
1862  ! Propose the next step for random walker 2
1863  call random_normal(3, sample_normal)
1864  x_n_2 = x_2(:, i-1) + sigma * sample_normal
1865 
1866  ! Calculate the transition probability between the proposed steps and current steps
1867  log_prob_n = log(h_2_prob(a, beta, x_1(:, i-1), x_2(:, i-1), x_n_1, x_n_2, r_a, r_b))
1868 
1869  ! Calculate the acceptance ratio
1870  a_prob_trans = min(1.0_real64, exp(log_prob_n))
1871 
1872  ! Generate sample from uniform distribution
1873  call ran1(unif_num, seed)
1874 
1875  ! Compare the drawn sample from the uniform distribution with the acceptance ratio
1876  if (a_prob_trans >= unif_num) then
1877  ! If accepted calculate local energy and add proposed position to the MCMC chain
1878  x_1(:, i) = x_n_1
1879  x_2(:, i) = x_n_2
1880 
1881  ! To determine whether this point should be added to the burnt and thinned MCMC chain
1882  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1883  ! Increment the reset counter
1884  reset_counter = 1 + reset_counter
1885  ! Precalculate i_bt for shorthand
1886  i_bt = i - burn - (reset_counter * (thin - 1))
1887  ! Add to the burn and thinned MCMC chain
1888  x_1_burnt(:, i_bt) = x_n_1
1889  x_2_burnt(:, i_bt) = x_n_2
1890  ! Calculate the local energy at this point
1891  call h_2_energy(a, beta, x_n_1, x_n_2, r_a, r_b, e_loc, dphidbeta)
1892  energy_chain(i_bt) = e_loc
1893  dphidbeta_chain(i_bt) = dphidbeta
1894  ! Resetting the thin counter
1895  thin_counter = 0
1896 
1897  ! Increment the number of accepted steps
1898  accept = 1 + accept
1899  acceptance = 'Y'
1900 
1901  if (i_log) then
1902  if (mod(i-burn, re_size/logfrac) .le. 1e-6_real64) then
1903  write(logid, "(I8, ' ', F12.8, ' ', F12.8, ' ', F12.8,&
1904  & ' ', F12.8, ' ', F12.8, ' ', F12.8, ' ', F12.6, ' ', A1)") &
1905  (i_bt), x_1_burnt(1, i_bt), x_1_burnt(2, i_bt), x_1_burnt(3, i_bt),&
1906  x_2_burnt(1, i_bt), x_2_burnt(2, i_bt), x_2_burnt(3, i_bt),&
1907  energy_chain(i_bt), acceptance
1908  end if
1909  end if
1910  end if
1911  else
1912  ! If rejected the proposed step is the current step
1913  x_1(:, i) = x_1(:, i-1)
1914  x_2(:, i) = x_2(:, i-1)
1915 
1916  ! To determine whether this point should be added to the burnt and thinned MCMC chain
1917  if ((i > (burn + 1)) .and. (thin_counter == thin)) then
1918  ! Increment the reset counter
1919  reset_counter = 1 + reset_counter
1920  ! Precalculate i_bt for shorthand
1921  i_bt = i - burn - (reset_counter * (thin - 1))
1922  ! Add to the burn and thinned MCMC chain
1923  x_1_burnt(:, i_bt) = x_1(:, i-1)
1924  x_2_burnt(:, i_bt) = x_2(:, i-1)
1925  ! Calculate the local energy at this point
1926  call h_2_energy(a, beta, x_1(:, i-1), x_2(:, i-1), r_a, r_b, e_loc, dphidbeta)
1927  energy_chain(i_bt) = e_loc
1928  dphidbeta_chain(i_bt) = dphidbeta
1929  ! Resetting the thin counter
1930  thin_counter = 0
1931  end if
1932  end if
1933  end do
1934  end if
1935 
1936  ! Calculate the acceptance ratio
1937  accept_rate = (real(accept, kind=real64) / real(re_size, kind=real64))
1938 
1939  ! Calculate the mean energy and the variance in the energy
1940  energy_mean = sum(energy_chain) / real(re_size, kind=real64)
1941  energy_variance = 0.0_real64
1942  do i = 1, re_size
1943  energy_variance = energy_variance + ((1.0_real64 / real(re_size - 1, kind=real64)) * &
1944  (energy_chain(i) - energy_mean) * (energy_chain(i) - energy_mean))
1945  end do
1946  energy_variance = (1.0_real64 / real(re_size, kind=real64)) * energy_variance
1947 
1948  if (i_log) then
1949  write(logid, "('Local energy = ', F12.8)") energy_mean
1950  write(logid, "('Variance = ', F12.8)") energy_variance
1951  write(logid, "('Acceptance rate = ', F12.8)") accept_rate
1952  write(logid, "('')")
1953  end if
1954 
1955  end subroutine vqmc_h_2
1956 
1957  ! ----------------------------------------------------------------------------------------------------------------
1958  ! ROUTINE: Auto_H_2
1959  !
1960  ! DESCRIPTION:
1964  !
1969  !
1970  ! PARAMETERS:
1979  ! ----------------------------------------------------------------------------------------------------------------
1980  subroutine auto_h_2(a, beta_in, bond_length, beta_out, energy_out, variance_out, &
1981  accept_out, converged)
1982 
1983  implicit none
1984 
1985  real(real64), intent(in) :: a
1986  real(real64), intent(in) :: beta_in
1987  real(real64), intent(in) :: bond_length
1988  real(real64), intent(out) :: beta_out
1989  real(real64), intent(out) :: energy_out
1990  real(real64), intent(out) :: variance_out
1991  real(real64), intent(out) :: accept_out
1992  character(len=1), intent(out) :: converged
1993 
1994  real(real64), dimension(3) :: R_a, R_b
1995  real(real64) :: beta_current, beta_new
1996  real(real64) :: beta_grad
1997  real(real64), dimension(:, :), allocatable :: e1_pos_chain, e2_pos_chain
1998  real(real64), dimension(:), allocatable :: energy_chain, DphiDbeta_chain
1999  integer :: viz_steps
2000  integer :: k
2001  integer :: ierr
2002  integer :: maxiter = 200
2003  real(real64) :: tol = 1e-4_real64
2004  character(len=6) :: str_b
2005  character(len=40) :: chainfile
2006 
2007  ! MPI Variables
2008  real(real64) :: my_beta_grad, my_ensum
2009  real(real64) :: sum_var, sum_acc
2010 
2011  ! Initialise
2012  beta_current = beta_in
2013  r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
2014  r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
2015 
2016  if (i_log) then
2017  write(logid, "('Starting parameter search for beta...')")
2018  end if
2019 
2020  ! Main iteration loop
2021  do k = auto_start, maxiter
2022 
2023  ! Calculate gradient at current beta
2024  call vqmc_h_2(h2%sigma, my_steps+h2%burn_step, h2%burn_step, h2%thin_step, a, &
2025  beta_current, r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
2026  energy_chain, energy_out, variance_out, dphidbeta_chain)
2027 
2028  ! Find gradient
2029  my_beta_grad = 2.0_real64 * ((sum(energy_chain * dphidbeta_chain) / size(energy_chain)) &
2030  - energy_out * (sum(dphidbeta_chain) / size(energy_chain)))
2031 
2032  ! Find average beta_grad across ranks, on all ranks.
2033  call mpi_allreduce(my_beta_grad,beta_grad,1,mpi_double_precision,mpi_sum,mpi_comm_world,mpierr)
2034  beta_grad = beta_grad * per_proc
2035 
2036  ! Update the log
2037  if (i_log) then
2038  write(logid, "('Beta Gradient')")
2039  write(logid, "(F13.7, ' ', F13.7)") beta_current, beta_grad
2040  end if
2041 
2042  ! Check if gradient is within tolerance.
2043  if (abs(beta_grad) .le. tol) then
2044 
2045  ! Update beta_out, and note convergence, on all ranks.
2046  beta_out = beta_current
2047  converged = 'Y'
2048 
2049  ! Sum the chain energy-sums and store in rank 0
2050  my_ensum = sum(energy_chain)
2051  call mpi_reduce(my_ensum,energy_out,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
2052 
2053  ! Sum the chain variances and store in rank 0
2054  call mpi_reduce(variance_out, sum_var,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
2055 
2056  ! Sum the chain acceptance rates and store in rank 0
2057  call mpi_reduce(accept_out, sum_acc,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
2058 
2059  ! Calculate averages on rank 0
2060  if (my_rank==0) then
2061 
2062  energy_out = energy_out * per_step
2063  variance_out = sum_var * per_proc2
2064  accept_out = sum_acc * per_proc
2065 
2066  if (write_log) then
2067  write(logid, "('Parameter search converged!')")
2068  write(logid, "('Results:')")
2069  write(logid, "('Beta = ', F13.7, ', Energy = ', F12.8, ', Variance = ', F12.8)") &
2070  beta_current, energy_out, variance_out
2071  write(logid, "('')")
2072  end if
2073 
2074  endif
2075 
2076  ! Exit subroutine, as convergence is reached.
2077  ! Note: Only rank 0 has the correct array outputs. All ranks have correct beta_out and 'converged'.
2078  exit
2079 
2080  end if
2081 
2082  ! If not, update beta by steepest descent, on all ranks.
2083  call h_2_update_beta(beta_current, beta_grad, beta_new)
2084  beta_current = beta_new
2085 
2086  ! Write restart files, if required
2087  if (write_restart .and. my_rank==0) then
2088  if (mod(k,restart_num)==0) then
2089  ! Write res_nml.txt
2090  call write_res_nml(i_bond,k+1,h2%auto_params,beta_current)
2091  endif
2092  endif
2093 
2094  end do
2095 
2096  ! Deal with non-convergence
2097  if (k .ge. maxiter) then
2098  converged = 'N'
2099  beta_out = beta_current
2100  if (i_log) then
2101  write(logid, "('Parameter search NOT converged! Consider changing&
2102  & the starting value of beta, or the damping in the H_2_update_beta&
2103  & function.')")
2104  end if
2105  end if
2106 
2107  ! Report results
2108  if (my_rank==0) then
2109  write(*, "(F12.8, ' ', A, ' ', F13.7, ' ', F12.8)") &
2110  bond_length, converged, beta_out, energy_out
2111  endif
2112 
2113  ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
2114  if (my_rank == 0) then
2115  if (write_chains) then
2116  write(str_b, "(F6.3)") bond_length
2117  write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
2118  write(*, "('Running VQMC at beta = ', F12.7, ' to produce single Markov chain...')") beta_out
2119  viz_steps = min(h2plus%steps,max_chain_len)
2120  call vqmc_h_2(h2%sigma, viz_steps, h2%burn_step, h2%thin_step, a, &
2121  beta_out, r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
2122  energy_chain, energy_out, variance_out, dphidbeta_chain)
2123  write(*, "(2A)") 'Writing Markov chain to ', chainfile
2124  call write_h2_equilibration(e1_pos_chain, energy_chain, chainfile, ierr, &
2125  h2, accept_out, beta_out, bond_length)
2126  end if
2127  end if
2128 
2129  end subroutine auto_h_2
2130 
2131  ! ----------------------------------------------------------------------------------------------------------------
2132  ! ROUTINE: Grid_H_2
2133  !
2134  ! DESCRIPTION:
2138  !
2143  !
2144  ! PARAMETERS:
2152  ! ----------------------------------------------------------------------------------------------------------------
2153  subroutine grid_h_2(a, beta_array, bond_length, beta_out, energy_out, variance_out, &
2154  accept_out)
2155 
2156  implicit none
2157 
2158  real(real64), intent(in) :: a
2159  real(real64), dimension(:), intent(in) :: beta_array
2160  real(real64), intent(in) :: bond_length
2161  real(real64), intent(out) :: beta_out
2162  real(real64), intent(out) :: energy_out
2163  real(real64), intent(out) :: variance_out
2164  real(real64), intent(out) :: accept_out
2165 
2166  real(real64), dimension(3) :: R_a, R_b
2167  real(real64), dimension(:, :), allocatable :: e1_pos_chain, e2_pos_chain
2168  real(real64), dimension(:), allocatable :: energy_chain, DphiDbeta_chain
2169  real(real64), dimension(:), allocatable :: param_energies, param_accepts, param_variances
2170  integer :: viz_steps
2171  integer :: i
2172  integer :: ierr
2173  character(len=6) :: str_b
2174  character(len=40) :: chainfile
2175 
2176  ! MPI variables
2177  real(real64) :: my_ensum
2178 
2179  ! Initialise
2180  r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
2181  r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
2182  allocate(param_energies(size(beta_array))); param_energies(:)= 0.0_real64
2183  allocate(param_accepts(size(beta_array))); param_accepts(:) = 0.0_real64
2184  allocate(param_variances(size(beta_array))) ; param_variances(:) = 0.0_real64
2185 
2186  ! Sort out restarting
2187  if (run_restart) then
2188  ! Read in array values
2189  call read_energies(resfile_epar, param_energies, param_variances, param_accepts)
2190  ! Reset the tag
2191  run_restart = .false.
2192  endif
2193 
2194  if (i_log) then
2195  write(logid, "('Starting parameter search for beta...')")
2196  end if
2197 
2198  ! Loop over beta values
2199  do i = grid_start, size(beta_array)
2200 
2201  ! Calculate energy at current value of beta
2202  call vqmc_h_2(h2%sigma, my_steps+h2%burn_step, h2%burn_step, h2%thin_step, a, &
2203  beta_array(i), r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
2204  energy_chain, energy_out, variance_out, dphidbeta_chain)
2205 
2206  ! Save energy - serial version
2207  !param_energies(i) = energy
2208  !param_variances(i) = variance
2209  !param_accepts(i) = accept_rate
2210 
2211  ! Sum the chain energy-sums and store in rank 0
2212  my_ensum = sum(energy_chain)
2213  call mpi_reduce(my_ensum,param_energies(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
2214 
2215  ! Sum the chain variances and store in rank 0
2216  call mpi_reduce(variance_out,param_variances(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
2217 
2218  ! Sum the chain acceptance rates and store in rank 0
2219  call mpi_reduce(accept_out,param_accepts(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
2220 
2221  ! Calculate averages on rank 0
2222  if (my_rank==0) then
2223  param_energies(i) = param_energies(i) * per_step
2224  param_variances(i) = param_variances(i) * per_proc2
2225  param_accepts(i) = param_accepts(i) * per_proc
2226 
2227  ! Write restart files, if required
2228  if (write_restart) then
2229  ! Write res_nml.txt
2230  call write_res_nml(i_bond,i+1,h2%auto_params)
2231  ! Write param energies, variances and acceptance rates
2232  call write_qho_main(beta_array, param_energies, param_variances, resfile_epar, ierr, qho, param_accepts)
2233  endif
2234  endif
2235 
2236  end do
2237 
2238  ! Note: Only rank 0 has the correct _out values.
2239  if (my_rank==0) then
2240 
2241  ! Select lowest energy and corresponding beta value
2242  energy_out = minval(param_energies)
2243  beta_out = beta_array(minloc(param_energies, 1))
2244  variance_out = param_variances(minloc(param_energies, 1))
2245  accept_out = param_accepts(minloc(param_energies, 1))
2246 
2247  ! Report results
2248  write(*, "(F12.8, ' - ', F13.7, ' ', F12.8)") &
2249  bond_length, beta_out, energy_out
2250 
2251  if (write_log) then
2252  write(logid, "('Results:')")
2253  write(logid, "('beta Energy')")
2254  do i = 1, size(beta_array)
2255  write(logid, "(F12.8, ' ', F12.8)") beta_array(i), param_energies(i)
2256  end do
2257  write(logid, "('')")
2258  write(logid, "('Min Energy = ', F12.8, ' at c = ', F12.8)") energy_out, beta_out
2259  write(logid, "('')")
2260  end if
2261 
2262  ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
2263  if (write_chains) then
2264  write(str_b, "(F6.3)") bond_length
2265  write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
2266  write(*, "('Running VQMC at beta = ', F12.7, ' to produce single Markov chain...')") beta_out
2267  viz_steps = min(h2plus%steps,max_chain_len)
2268  call vqmc_h_2(h2%sigma, viz_steps, h2%burn_step, h2%thin_step, a, &
2269  beta_out, r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
2270  energy_chain, energy_out, variance_out, dphidbeta_chain)
2271  write(*, "(2A)") 'Writing Markov chain to ', chainfile
2272  call write_h2_equilibration(e1_pos_chain, energy_chain, chainfile, ierr, &
2273  h2, accept_out, beta_out, bond_length)
2274  end if
2275 
2276  endif
2277 
2278  end subroutine grid_h_2
2279 
2280  ! ----------------------------------------------------------------------------------------------------------------
2281  ! ROUTINE: Equil_H_2
2282  !
2283  ! DESCRIPTION:
2287  !
2292  !
2293  ! PARAMETERS:
2298  ! ----------------------------------------------------------------------------------------------------------------
2299  subroutine equil_h_2(a, beta_array, bond_length, bonds_pos)
2300 
2301  implicit none
2302 
2303  real(real64), intent(in) :: a
2304  real(real64), dimension(:), intent(in) :: beta_array
2305  real(real64), intent(in) :: bond_length
2306  integer, intent(in) :: bonds_pos
2307 
2308  real(real64), dimension(3) :: R_a, R_b
2309  real(real64), dimension(:), allocatable :: variances
2310  real(real64), dimension(:), allocatable :: energies
2311  real(real64), dimension(:, :), allocatable :: e1_pos_chain, e2_pos_chain
2312  real(real64), dimension(:), allocatable :: energy_chain, DphiDbeta_chain
2313  real(real64) :: accept_rate
2314  real(real64) :: energy_mean
2315  real(real64) :: energy_variance
2316  integer :: i, ierr
2317  character(len=3) :: str_i, str_b
2318  character(len=40) :: equilfile
2319 
2320  ! Initialise
2321  r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
2322  r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
2323  allocate(energies(size(beta_array))); energies(:) = 0.0_real64
2324  allocate(variances(size(beta_array)))
2325 
2326  ! Loop over values of beta
2327  do i = grid_start, size(beta_array)
2328 
2329  ! Calculate energy at current value of beta
2330  call vqmc_h_2(h2%sigma, h2%steps, h2%burn_step, h2%thin_step, a, &
2331  beta_array(i), r_a, r_b, e1_pos_chain, e2_pos_chain, accept_rate, &
2332  energy_chain, energy_mean, energy_variance, dphidbeta_chain)
2333  ! Save energy
2334  energies(i) = energy_mean
2335 
2336  write(*, "(F12.8, ' - ', F13.7, ' ', F12.8)") &
2337  bond_length, beta_array(i), energies(i)
2338 
2339  ! Save the energy trace for each value of c
2340  write(str_i, "(I3)") i
2341  write(str_b, "(I3)") bonds_pos
2342  write(equilfile, "(5A)") './equil_out/b', trim(adjustl(str_b)), &
2343  '_beta', trim(adjustl(str_i)), '.nc'
2344  write(*, *) 'Writing equilibration outfile to ', equilfile
2345  call write_h2_equilibration(e1_pos_chain, energy_chain, equilfile, ierr, &
2346  h2, accept_rate, beta_array(i), bond_length)
2347 
2348  ! Write restart files, if required
2349  if (write_restart) then
2350  if (mod(i,restart_num)==0) then
2351  call write_res_nml(i_bond,i+1,h2%auto_params)
2352  endif
2353  endif
2354 
2355  end do
2356 
2357  if (write_log) then
2358  write(logid, "('Results:')")
2359  write(logid, "('beta Energy')")
2360  do i = 1, size(beta_array)
2361  write(logid, "(F12.8, ' ', F12.8)") beta_array(i), energies(i)
2362  end do
2363  write(logid, "('')")
2364  end if
2365 
2366  deallocate(energies)
2367  deallocate(variances)
2368 
2369  end subroutine equil_h_2
2370 
2371 end module VQMC
constants
Contains mathematical constants needed throughout the code, all of which are globally accessible.
Definition: shared_data.f90:16
restart_fns
Contains routines to read and write restart files.
Definition: restart_fns.f90:15
shared_data
Contains derived types, and global variables to store input values of simulation parameters.
Definition: shared_data.f90:49
shared_data::seed
integer(int64) seed
Seed for the random number generator.
Definition: shared_data.f90:152
solvers
Contains local energy and transition probability solvers for available problems.
Definition: quantum_solvers.f90:16
vqmc
Contains routines for handline random number generation, the main VQMC algorithms and routines for de...
Definition: vqmc.f90:30
vqmc::ran1
subroutine ran1(sample, dseed)
Implementation of Press et al's random number generator.
Definition: vqmc.f90:119
vqmc::vqmc_qho
subroutine vqmc_qho(sigma, N, burn, thin, alpha_const, X_burnt, accept_rate, energy_chain, energy_mean, energy_variance, x_start)
Performs VQMC for the QHO system.
Definition: vqmc.f90:291
vqmc::grid_h_2
subroutine grid_h_2(a, beta_array, bond_length, beta_out, energy_out, variance_out, accept_out)
Performs a grid search to find the optimal value for the variational parameter in the problem.
Definition: vqmc.f90:2155
vqmc::auto_h_2_plus
subroutine auto_h_2_plus(c_in, bond_length, c_out, energy_out, variance_out, accept_out, converged)
Automatically finds the best variational parameter for the problem at given bond length.
Definition: vqmc.f90:1190
vqmc::vqmc_h_2
subroutine vqmc_h_2(sigma, N, burn, thin, a, beta, R_a, R_b, X_1_burnt, X_2_burnt, accept_rate, energy_chain, energy_mean, energy_variance, DphiDbeta_chain, x_start_1, x_start_2)
Performs VQMC for the system.
Definition: vqmc.f90:1620
vqmc::iy
integer(int64), save iy
Variable controlling the RNG sequence.
Definition: vqmc.f90:46
vqmc::equil_h_2_plus
subroutine equil_h_2_plus(c_array, bond_length, bonds_pos)
Performs VQMC on an array of values for the variational parameter , outputting energy traces for each...
Definition: vqmc.f90:1506
vqmc::ntab
integer, parameter ntab
Variable controlling the RNG sequence.
Definition: vqmc.f90:44
vqmc::logfrac
integer logfrac
1/logfrac is the fraction of samples to report in the log, if shared_data::write_log is TRUE.
Definition: vqmc.f90:49
vqmc::equil_qho
subroutine equil_qho(alpha_array)
Performs VQMC on an array of values for the variational parameter , outputting energy traces for each...
Definition: vqmc.f90:743
vqmc::random_normal
subroutine random_normal(dim, x)
Draws a sample from a normal distribution using the Box-Muller transform.
Definition: vqmc.f90:229
vqmc::linspace
real(real64) function, dimension(:), allocatable linspace(start, end, num)
Adaptation of numpy's linspace function.
Definition: vqmc.f90:183
vqmc::grid_h_2_plus
subroutine grid_h_2_plus(c_array, bond_length, c_out, energy_out, variance_out, accept_out)
Performs a grid search to find the optimal value for the variational parameter in the problem.
Definition: vqmc.f90:1361
vqmc::auto_h_2
subroutine auto_h_2(a, beta_in, bond_length, beta_out, energy_out, variance_out, accept_out, converged)
Automatically finds the best variational parameter for the problem at a given bond length.
Definition: vqmc.f90:1982
vqmc::grid_qho
subroutine grid_qho(alpha_array, energies, acc_rates, variances)
Performs VQMC on an array of values for the variational parameter .
Definition: vqmc.f90:580
vqmc::vqmc_h_2_plus
subroutine vqmc_h_2_plus(sigma, N, burn, thin, c_const, R_a, R_b, X_burnt, accept_rate, energy_chain, energy_mean, energy_variance, DphiDc_chain, x_start)
Performs VQMC for the system.
Definition: vqmc.f90:860
vqmc::iv
integer(int64), dimension(ntab), save iv
Variable controlling the RNG sequence.
Definition: vqmc.f90:45
vqmc::equil_h_2
subroutine equil_h_2(a, beta_array, bond_length, bonds_pos)
Performs VQMC on an array of values for the variational parameter , outputting energy traces for each...
Definition: vqmc.f90:2300
vqmc::init_ran1
subroutine init_ran1()
Initialiser for the random number generator.
Definition: vqmc.f90:69
write_netcdf
Contains NetCDF output functions for different problem main & equilibration runs.
Definition: netcdf_out.f90:16