vqmc Module Reference

Contains routines for handline random number generation, the main VQMC algorithms and routines for determining ideal parameters. More...

Functions/Subroutines
subroutine	init_ran1 ()
	Initialiser for the random number generator. More...
subroutine	ran1 (sample, dseed)
	Implementation of Press et al's random number generator. More...
real(real64) function, dimension(:), allocatable	linspace (start, end, num)
	Adaptation of numpy's linspace function. More...
subroutine	random_normal (dim, x)
	Draws a sample from a normal distribution using the Box-Muller transform. More...
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. More...
subroutine	grid_qho (alpha_array, energies, acc_rates, variances)
	Performs VQMC on an array of values for the variational parameter \( \alpha \). More...
subroutine	equil_qho (alpha_array)
	Performs VQMC on an array of values for the variational parameter \( \alpha \), outputting energy traces for each. More...
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 \( H_{2}^{+} \) system. More...
subroutine	auto_h_2_plus (c_in, bond_length, c_out, energy_out, variance_out, accept_out, converged)
	Automatically finds the best variational parameter \( c \) for the \( H_{2}^{+} \) problem at given bond length. More...
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 \( c \) in the \( H_{2}^{+} \) problem. More...
subroutine	equil_h_2_plus (c_array, bond_length, bonds_pos)
	Performs VQMC on an array of values for the variational parameter \( c \), outputting energy traces for each. More...
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 \( H_{2}^{+} \) system. More...
subroutine	auto_h_2 (a, beta_in, bond_length, beta_out, energy_out, variance_out, accept_out, converged)
	Automatically finds the best variational parameter \( \beta \) for the \( H_{2}^{+} \) problem at a given bond length. More...
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 \( \beta \) in the \( H_{2} \) problem. More...
subroutine	equil_h_2 (a, beta_array, bond_length, bonds_pos)
	Performs VQMC on an array of values for the variational parameter \( \beta \), outputting energy traces for each. More...

Variables
integer, parameter	ntab = 32
	Variable controlling the RNG sequence. More...
integer(int64), dimension(ntab), save	iv = 0
	Variable controlling the RNG sequence. More...
integer(int64), save	iy = 0
	Variable controlling the RNG sequence. More...
integer	logfrac = 10
	1/logfrac is the fraction of samples to report in the log, if shared_data::write_log is TRUE. More...

Detailed Description

Contains routines for handline random number generation, the main VQMC algorithms and routines for determining ideal parameters.

This module contains the VQMC solvers for the quantum harmonic oscillator (QHO), hydrogen ion (H_2_plus) and hydrogen molecule (H_2) problems. It also has the subroutine random normal which generates a sample from a normal distribution using the Box-Muller transfrom, which is needed for the Metropolis algorithm used in VQMC. This uses a random number generator which is initially seeded by environmental noise through /dev/urandom, and then propagated by an algorithm from Numerical Recipes by Press et al. (2007)

The VQMC algorithms use the solvers module to obtain the transition probabilities needed for the Metropolis algorithm and to calculate the local energies given the position of the random walker(s).

The module also contains routines for running these VQMC algorithms across a range of parameters. The Grid routines are used to do a search for the best variational parameter on a user-defined grid, while the Auto routines utilise the gradients calculated by the local energy routines to automatically determine the best variational parameter from a suitable guess. Equil routines are also included to run VQMC over a parameter grid, but to output the position and energy traces from VQMC rather than a total energy. These traces can then be used to determine optimal burning/thinning parameters for the Markov chains, using tools in the Python frontend.

Function/Subroutine Documentation

init_ran1()

subroutine vqmc::init_ran1

Initialiser for the random number generator.

Initialises the random number generator by drawing a seed from the environmental noise in /dev/urandom. Cuts down this seed to a smaller integer to prevent overflows in ran1 when the seed is converted to a double-precision float.

Definition at line 68 of file vqmc.f90.

 
        implicit none
 
        integer :: fu       ! File unit
        integer :: rc       ! I/O status of read from /dev/urandom
        integer :: l_seed   ! Length of unshortened seed
        character(len=30) :: c_seed     ! Character value of seed
 
        ! Read pseudorandom noise from /dev/urandom
        open (access='stream', action='read', file='/dev/urandom', &
        form='unformatted', iostat=rc, newunit=fu)
 
        if (rc == 0) read (fu) seed
 
        close (fu)
 
        ! Convert long seed to characters, shorten and convert back into an integer
        write(c_seed, "(I30)") seed
        l_seed = len(c_seed)
        c_seed = adjustl(c_seed)
        c_seed = c_seed(7:l_seed)
        read(c_seed, *) seed
 
        ! Ensure seed is negative for compatability with ran1
        if (seed .gt. 0) then
            seed = -seed
        end if
 

Here is the caller graph for this function:

ran1()

subroutine vqmc::ran1	(	real(real64), intent(out)	sample,
		integer(int64), intent(inout)	dseed
	)

Implementation of Press et al's random number generator.

Update to David Quigley's implementation of Press et al's random number generator (originally from Numerical Methods ), as seem in PX425 Assignment 3, to be Fortran 2008 compliant.

Remarks

NOTE: The variable dseed must be a negative integer when this routine is first called, otherwise the generator will not initialise properly.

Parameters

[in,out]	dseed	Current value of seed. Must be a negative integer on first call of this routine.
[out]	sample	Uniform random number.

Definition at line 118 of file vqmc.f90.

 
        implicit none
 
        integer(int64), intent(inout) :: dseed
        real(real64), intent(out) :: sample
 
        ! Constants defined in Press et al
        integer(int64), parameter :: ia = 16807, im = 2147483647, iq = 127773
        integer(int64), parameter :: ir = 2836, &
                                     ndiv = 1 + floor((real(im, kind=real64)-1.0)/ntab)
        real(real64), parameter :: am = 1.0/real(im, kind=real64), eps = 1.2e-7, &
                                   rnmx = 1.0-eps
 
        ! Loop counters
        integer(int64) :: j, k
 
        ! The following is transcribed directly from Press et al
        if (dseed<=0 .or. iy==0) then ! initialization
                                      ! NB this may need to be modified to allow for restarting as iy will be 0
            dseed = max(-dseed,1)
 
            do j = ntab+8, 1, -1
                k = floor(real(dseed, kind=real64)/real(iq, kind=real64))
                dseed = ia * (dseed-k*iq) - ir * k
                if (dseed < 0) dseed = dseed + im
                if (j <= ntab) iv(j) = dseed
            end do
 
            iy = iv(1)
        end if
 
        k = floor(real(dseed, kind=real64)/real(iq, kind=real64))
        dseed = ia * (dseed-k*iq) - ir * k
        if (dseed < 0) dseed = dseed + im
        j = 1_int64 + iy / ndiv
        iy = iv(j)
        iv(j) = dseed
        sample = min(am*iy, rnmx)
 

Here is the caller graph for this function:

linspace()

real(real64) function, dimension(:), allocatable vqmc::linspace	(	real(real64), intent(in)	start,
		real(real64), intent(in)	end,
		real(real64), intent(in)	num
	)

Adaptation of numpy's linspace function.

Uses a REAL value for num for compatability with parameter arrays, then converts to int internally. Returns a length 1 array just containing the start value if num is equal to 1.

Parameters

[in]	start	Starting value of the sequence
[in]	end	End value of the sequence
[in]	num	Number of samples to generate between start and end value

Return values

ls_array

An array with evenly specified numbers over an interval

Definition at line 182 of file vqmc.f90.

 
        implicit none
 
        real(real64), intent(in) :: start, end, num
 
        real(real64), dimension(:), allocatable :: ls_array
 
        integer :: inum                                         ! The integer type of num
        integer :: i                                            ! Loop variable 
        real(real64) :: step                                    ! Spacing between two numbers in the array
 
        ! Handle real value for num
        inum = int(num)
        ! Handle single values
        if (inum .eq. 1) then
            allocate(ls_array(inum))
            ls_array(1) = start
        ! Handle multiple values
        else
            ! Generate array
            allocate(ls_array(inum))
            ! Find step size
            step = (end - start) / real(inum-1, kind=real64)
            ! Fill array
            do i = 1, inum
                ls_array(i) = start + (i-1)*step
            end do
        end if
 

Here is the caller graph for this function:

random_normal()

subroutine vqmc::random_normal	(	integer, intent(in)	dim,
		real(real64), dimension(:), intent(out)	x
	)

Draws a sample from a normal distribution using the Box-Muller transform.

Uses the uniform random number generator ran1 to generate an array of normally distributed numbers, by means of the Box-Muller transform.

Parameters

[in]	dim	The dimension of the random vector to be generated.
[out]	x	The vector sample from the normal distribution.

Definition at line 228 of file vqmc.f90.

        
        implicit none 
 
        integer, intent(in) :: dim
        real(real64), intent(out), dimension(:) :: x 
        
        real(real64) :: u_1, u_2                        ! Samples from uniform distribution 
        integer :: i                                    ! Loop variable 
 
        ! Generate vector containg samples from normal distributions 
        do i = 1, dim
            ! This can be done the Fortran intrinsic random number generator but this 
            ! is not a good random number generator 
            ! call random_number(u_1)
            ! call random_number(u_2)
 
            ! Draws from uniform distribution
            call ran1(u_1, seed)
            call ran1(u_2, seed)
 
            ! Calculate the sample from the normal distribution using the Box-Muller transform 
            x(i) = sqrt(-2*log(u_1))*cos(2*pi*u_2)
        end do 
 

Here is the call graph for this function:

Here is the caller graph for this function:

vqmc_qho()

subroutine vqmc::vqmc_qho	(	real(real64), intent(in)	sigma,
		integer, intent(in)	N,
		integer, intent(in)	burn,
		integer, intent(in)	thin,
		real(real64), intent(in)	alpha_const,
		real(real64), dimension(:), intent(out), allocatable	X_burnt,
		real(real64), intent(out)	accept_rate,
		real(real64), dimension(:), intent(out), allocatable	energy_chain,
		real(real64), intent(out)	energy_mean,
		real(real64), intent(out)	energy_variance,
		real(real64), optional	x_start
	)

Performs VQMC for the QHO system.

Computes the total energy of a trial wavefunction by averaging over a Markov chain of local energies generated by QHO_Energy, with random move acceptance/rejection probability determined by QHO_Prob.

Allows for modification of total Markov chain length, as well as the number of steps to be 'burned' while the chain is equilibrating and the interval to 'thin' the chain over to reduce correlation of samples. The size of the proposed random moves can be controlled by modifying the sigma parameter to get an ideal move acceptance rate.

Parameters

[in]	sigma	Step size
[in]	N	Number of steps to take
[in]	burn	Number of sample to discard at the beginning
[in]	thin	The frequency at which we discard the samples from the MCMC chain after burning
[in]	alpha_const	Parameter
[in]	x_start	Starting point of MCMC chain (optional, default is randomly generated)
[out]	X_burnt	The burnt and thinned MCMC chain for the random walker
[out]	accept_rate	The acceptance rate of the MCMC chain
[out]	energy_chain	An array with the local energies of the each points in the MCMC chain
[out]	energy_mean	The average of the local energies recorded in energy_chain
[out]	energy_variance	The variance of the local energies recorded in energy_chain

Definition at line 289 of file vqmc.f90.

 
        implicit none 
 
        real(real64), intent(in) :: sigma, alpha_const
        integer, intent(in) :: N, burn, thin 
        real(real64), optional :: x_start 
         
        real(real64), intent(out) :: accept_rate, energy_mean, energy_variance
        real(real64), intent(out), dimension(:), allocatable :: X_burnt, energy_chain 
        
        real(real64) :: log_prob_n                      ! The natural logarithm of the transition probability
        real(real64) :: a_prob_trans                    ! Acceptance ratio of the transition probability 
        real(real64) :: unif_num                        ! A random sample from the uniform distribution 
        real(real64), dimension(1) :: start             ! The starting point of the MCMC chain 
        real(real64), dimension(1) :: X_n               ! The proposed step 
        real(real64), dimension(1) :: sample_normal     ! A sample from the normal distribution 
        real(real64), dimension(:), allocatable :: X    ! Non-burnt or thinned MCMC chain
        integer :: thin_counter                         ! Counter to keep track of the thinning process
        integer :: reset_counter                        ! Counts how many times thin counter was reset
        integer :: accept                               ! Number of accepted steps 
        integer :: i                                    ! Loop variable 
        integer :: i_bt                                 ! Burned/thinned value of i
        integer :: re_size                              ! The size of the burnt and thinned array
        character(len=1) :: acceptance                  ! Yes/No signpost for acceptance
 
        ! Adjust the output array to be of the burnt and thinned array
        re_size = (n - burn) / thin
        
        allocate(x(1: n+1)); x(:) = 0.0_real64
        allocate(x_burnt(1: re_size)); x_burnt(:) = 0.0_real64
        allocate(energy_chain(1: re_size)); energy_chain(:) = 0.0_real64
        
        ! Allow the starting position to take the defualt value 
        if (present(x_start)) then
            start = x_start
            x(1) = start(1)
            if (i_log) then
                write(logid, "('Starting from provided Markov chain position', F12.8)") x(1)
            end if
        else 
            ! Creates the initial position of the random walker that is within a certain radius  
            call random_normal(1, sample_normal)
            do while(norm2(sample_normal) >= qho%rcut)
                call random_normal(1, sample_normal)
            end do 
            start = sample_normal
            x(1) = start(1)
            if (i_log) then
                write(logid, "('Generated random Markov chain start position', F12.8)") x(1)
            end if
        end if 
 
        ! Start the counter for the number of accepted moves
        accept = 0 
        ! Start the counter for the thinning process
        thin_counter = 0
        ! Start the reset counter 
        reset_counter = 0
 
        !Write initial conditions to the output file
        if (i_log) then
            write(logid, "('Markov chain for alpha = ', F7.4)") alpha_const
            write(logid, "('Sample     Position       Local Energy     Accepted?')")
            write(logid, "(I8, '   ', F12.8, '   ', F12.6)") 1, x(1), energy_chain(1)
        end if
 
        ! In case burn = 0 to start calculating the chain from the beginning 
        if (burn == 0) then
            energy_chain(1) = qho_energy(alpha_const, x(1))
            x_burnt(1) = x(1)
 
            do i = 2, n 
                ! Add one to the thin counter 
                thin_counter = 1 + thin_counter           
                ! Propose the next step
                call random_normal(1, sample_normal)
                x_n = x(i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between proposed step and current step
                log_prob_n = log(qho_prob(alpha_const, x(i-1), x_n(1)))
            
                ! Calculate the acceptance ratio
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
 
                ! Generate sample from uniform distribution 
                call ran1(unif_num, seed)
            
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio 
                if (a_prob_trans >= unif_num) then 
                    ! If accepted add proposed position to the MCMC chain
                    x(i) = x_n(1) 
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain
                    if (thin_counter == thin) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(i_bt) = x_n(1)
                        ! Calculate the local energy at this point 
                        energy_chain(i_bt) = qho_energy(alpha_const, x_n(1))
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'Y'
                        ! Increment the number of accepted steps 
                        accept = 1 + accept 
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.6, '     ', A1)") &
                                    (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
                            end if
                        end if
                    end if 
                else 
                    ! If rejected the proposed step is the current step
                    x(i) = x(i-1)
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain 
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(i_bt) = x(i-1)
                        ! Calculate the local energy at this point
                        energy_chain(i_bt) = qho_energy(alpha_const, x(i-1))
                        ! Reseting the thin counter 
                        thin_counter = 0 
 
                        acceptance = 'N'
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.6, '     ', A1)") &
                                    (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
                            end if
                        end if
 
                    end if 
                end if 
 
            end do 
        ! Loop to advance the chain to the final burn step and not save any of the outputs
        else if (burn >= 2) then
            do i = 2, burn
                ! Propose the next step
                call random_normal(1, sample_normal)
                x_n = x(i-1) + sigma * sample_normal
                ! Calculate the transition probability between proposed step and current step
                log_prob_n = log(qho_prob(alpha_const, x(i-1), x_n(1)))
                ! Calculate the acceptance ratio
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
                ! Generate sample from uniform distribution 
                call ran1(unif_num, seed)
            
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio 
                if (a_prob_trans >= unif_num) then 
                    ! If accepted add proposed position to the MCMC chain
                    x(i) = x_n(1) 
                else 
                    ! If rejected the proposed step is the current step
                    x(i) = x(i-1)
                end if
            end do 
        end if
       
        ! Metropolis Algorithm to calculate the burnt and thinned MCMC chain 
        if (burn >=1) then 
            do i = burn + 1, n+1 
                ! Add one to the thin counter 
                thin_counter = 1 + thin_counter           
                ! Propose the next step
                call random_normal(1, sample_normal)
                x_n = x(i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between proposed step and current step
                log_prob_n = log(qho_prob(alpha_const, x(i-1), x_n(1)))
            
                ! Calculate the acceptance ratio
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
 
                ! Generate sample from uniform distribution 
                call ran1(unif_num, seed)
            
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio 
                if (a_prob_trans >= unif_num) then 
                    ! If accepted add proposed position to the MCMC chain
                    x(i) = x_n(1) 
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(i_bt) = x_n(1)
                        ! Calculate the local energy at this point 
                        energy_chain(i_bt) = qho_energy(alpha_const, x_n(1))
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'Y'
                        ! Increment the number of accepted steps 
                        accept = 1 + accept 
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.6, '     ', A1)") &
                                    (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
                            end if
                        end if
 
                    end if 
 
                else 
                    ! If rejected the proposed step is the current step
                    x(i) = x(i-1)
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain 
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(i_bt) = x(i-1)
                        ! Calculate the local energy at this point
                        energy_chain(i_bt) = qho_energy(alpha_const, x(i-1))
                        ! Reseting the thin counter 
                        thin_counter = 0 
 
                        acceptance = 'N'
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.6, '     ', A1)") &
                                    (i_bt), x_burnt(i_bt), energy_chain(i_bt), acceptance
                            end if
                        end if
 
                    end if 
                end if 
            end do
        end if
 
        ! Calculate the acceptance ratio 
        accept_rate = (real(accept, kind=real64) / real(re_size, kind=real64))
 
        ! Calculate the mean energy and the variance in the energy
        energy_mean = sum(energy_chain) / real(re_size, kind=real64)
        energy_variance = 0.0_real64
        do i = 1, re_size
            energy_variance = energy_variance + ((1.0_real64 / real(re_size - 1, kind=real64)) * &
            (energy_chain(i) - energy_mean) * (energy_chain(i) - energy_mean)) 
        end do
        energy_variance = (1.0_real64 / real(re_size, kind=real64)) * energy_variance
 
        if (i_log) then
            write(logid, "('Local energy = ', F12.8)") energy_mean
            write(logid, "('Variance = ', F12.8)") energy_variance
            write(logid, "('Acceptance rate = ', F12.8)") accept_rate
            write(logid, "('')")
        end if
 

Here is the call graph for this function:

Here is the caller graph for this function:

grid_qho()

subroutine vqmc::grid_qho	(	real(real64), dimension(:), intent(in)	alpha_array,
		real(real64), dimension(:), intent(inout)	energies,
		real(real64), dimension(:), intent(inout)	acc_rates,
		real(real64), dimension(:), intent(inout)	variances
	)

Performs VQMC on an array of values for the variational parameter \( \alpha \).

Evaluates energy of the QHO trial wavefunction at each value of the parameter \( \alpha \) in alpha_array. Also evaluates variances in these energies as a measure of uncertainty. Writes these results to a NetCDF file.

Parameters

[in]	alpha_array	The array of alpha values to obtain energies for
[in,out]	energies	Array of size(alpha_array) to contain total energies
[in,out]	acc_rates	Array of size(alpha_array) to contain acceptance rates
[in,out]	variances	Array of size(alpha_array) to contain variances

Definition at line 579 of file vqmc.f90.

 
        implicit none
 
        real(real64), dimension(:), intent(in) :: alpha_array
        real(real64), dimension(:), intent(inout) :: energies
        real(real64), dimension(:), intent(inout) :: acc_rates
        real(real64), dimension(:), intent(inout) :: variances
 
        real(real64), dimension(:), allocatable :: energy_chain, trimmed_energy_chain
        real(real64), dimension(:), allocatable :: pos_chain, trimmed_pos_chain
        real(real64) :: accept_rate
        real(real64) :: energy_mean
        real(real64) :: energy_variance
        integer :: i, ierr
        character(len=3) :: str_i
        character(len=40) :: chainfile
 
        ! MPI variables
        real(real64)    :: my_ensum
 
        ! Set up printed output for the loop, from rank 0.
        if (my_rank==0) then
            write(*, "('-------------- Begin VQMC Runs --------------')")
            write(*, "('Alpha       Energy            Acceptance Rate')")
            if (write_log) then
                write(logid, "('-------------- Begin VQMC Runs --------------')")
                if (.not. run_equil) then
                    write(logid, "('')")
                    write(logid, "('Logging the chains from Rank 0.')")
                    write(logid, "('')")
                endif 
            end if
        endif 
 
        ! Loop over values of alpha
        do i = alpha_start, size(alpha_array)
 
            ! Calculate energy at current value of alpha
            call vqmc_qho(qho%sigma, my_steps+qho%burn_step, qho%burn_step, qho%thin_step, &
                alpha_array(i), pos_chain, accept_rate, energy_chain, energy_mean, &
                energy_variance)
            
            ! Save energy - serial version
            !energies(i) = energy_mean
            !variances(i) = energy_variance
            !acc_rates(i) = accept_rate
            !write(*, "(F7.4, '     ', F12.8, '      ', F12.8)") alpha_array(i), &
            !    energy_mean, accept_rate
 
            ! Sum energy sums across ranks and save to rank 0 
            my_ensum = sum(energy_chain)
            call mpi_reduce(my_ensum,energies(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Sum accept rates across ranks, to find average in rank 0
            call mpi_reduce(accept_rate,acc_rates(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Sum variances acros ranks, to find average in rank 0
            call mpi_reduce(energy_variance,variances(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Print results from all ranks for checking
            !write(*, "(F7.4, '     ', F12.8, '      ', F12.8, '      ', F12.8, '      ', I4)") alpha_array(i), &
            !    energy_mean, accept_rate, energy_variance, my_rank
 
            ! Compute and print combined results in rank 0
            if (my_rank==0) then
                !print*, 'tot_ensum=',tot_ensum,' ,tot_size=',tot_size
                energies(i) = energies(i) * per_step
                acc_rates(i) = acc_rates(i) * per_proc
                variances(i) = variances(i) * per_proc ! CHECK FORMULA!
                write(*, "(F7.4, '     ', F12.8, '      ', F12.8)") alpha_array(i), &
                    energies(i), acc_rates(i)
 
                ! Write energy chain from rank 0, if required
                if (write_chains) then
 
                    ! Setup filename
                    write(str_i, "(I3)") i
                    write(chainfile, "(3A)") './chains_out/alpha_', trim(adjustl(str_i)), '.nc'
 
                    ! Check chain length
                    if (size(energy_chain) < max_chain_len) then
                        write(*, "(2A)") 'Writing Markov chain to from rank 0 to', chainfile
                        call write_qho_equilibration(pos_chain, energy_chain, chainfile, ierr, &
                            qho, accept_rate, alpha_array(i))
                    else
                        ! Adjust chain length for writing.
                        ! Allocate temporary arrays
                        allocate(trimmed_pos_chain(max_chain_len))
                        allocate(trimmed_energy_chain(max_chain_len))
                        ! Fill in the temporary arrays
                        trimmed_pos_chain(:) = pos_chain(1:max_chain_len)
                        trimmed_energy_chain(:) = energy_chain(1:max_chain_len)
                        ! Write the trimmed temporary arrays
                        write(*, "(2A)") 'Writing trimmed Markov chain to from rank 0 to', chainfile
                        call write_qho_equilibration(trimmed_pos_chain, trimmed_energy_chain, chainfile, ierr, &
                            qho, accept_rate, alpha_array(i))
                        ! Deallocate the temporary arrays
                        deallocate(trimmed_pos_chain)
                        deallocate(trimmed_energy_chain)
                    endif
                end if
 
                ! Write restart files, if required
                if (write_restart) then
                    if (mod(i,restart_num)==0) then
                        ! Write res_nml.txt
                        call write_res_nml(i+1)
                        ! Store energy, variance & accept_rate arrays
                        call write_qho_main(alpha_array, energies, variances, resfile_etot, ierr, qho, acc_rates)
                    endif
                endif
 
            endif 
 
        end do
 
        if (my_rank==0) then
 
            write(*, "('--------------  End VQMC Runs  --------------')")
            if (write_log) then
                write(logid, "('--------------  End VQMC Runs  --------------')")
                write(logid, "('Results:')")
                write(logid, "('alpha            Energy')")
                do i = 1, size(alpha_array)
                    write(logid, "(F12.8, '    ', F12.8)") alpha_array(i), energies(i)
                end do
                write(logid, "('')")
                write(logid, "('Min Energy = ', F12.8, ' at c = ', F12.8)") &
                    minval(energies), alpha_array(minloc(energies, 1))
                write(logid, "('')")
            end if
 
            ! Write main results file.
            write(*, *) 'Calculations finished, writing results to NetCDF output file ', &
                        outfile
            if (write_log) then
                write(logid, *) 'Calculations finished, writing results to NetCDF &
                    &output file ', outfile
            end if
            call write_qho_main(alpha_array, energies, variances, outfile, ierr, qho, acc_rates)
 
        endif 
 

Here is the call graph for this function:

Here is the caller graph for this function:

equil_qho()

subroutine vqmc::equil_qho

(

real(real64), dimension(:), intent(in)

alpha_array

)

Performs VQMC on an array of values for the variational parameter \( \alpha \), outputting energy traces for each.

Generates VQMC position and energy traces of the QHO trial wavefunction at each value of the parameter \( \alpha \) in alpha_array. Outputs each of these traces to a separate NetCDF file, to be used in determining optimal burning/thinning paramters for an accurate full VQMC run.

Parameters

[in]	alpha_array	The array of alpha values to obtain energies for
[in,out]	energies	The energies

Definition at line 742 of file vqmc.f90.

 
        implicit none
 
        real(real64), dimension(:), intent(in)    :: alpha_array
 
        real(real64), dimension(:), allocatable :: energy_chain
        real(real64), dimension(:), allocatable :: pos_chain
        real(real64) :: energy
        real(real64) :: accept_rate
        real(real64) :: variance
        integer :: i, ierr
        character(len=3) :: str_i
        character(len=40) :: equilfile
 
        ! Set up printed output for the loop.
        write(*, "('-------------- Begin VQMC Runs --------------')")
        write(*, "('Alpha       Energy            Acceptance Rate')")
        if (write_log) then
            write(logid, "('-------------- Begin VQMC Runs --------------')")
        end if 
 
        ! Loop over values of alpha
        do i = alpha_start, size(alpha_array)
 
            ! Calculate energy at current value of alpha
            call vqmc_qho(qho%sigma,qho%steps, qho%burn_step, qho%thin_step, &
                alpha_array(i), pos_chain, accept_rate, energy_chain, energy, &
                variance)
 
            ! Report the energy 
            write(*, "(F7.4, '     ', F12.8, '      ', F12.8)") alpha_array(i), &
            energy, accept_rate
            if (write_log) then
                write(logid, "('')")
                write(logid, "('Results:')")
                write(logid, "('alpha = ', F12.8, ', Energy = ', F12.8)") alpha_array(i), energy
                write(logid, "('')")
            end if
        
            ! Save the energy trace for each value of alpha.
            write(str_i, "(I3)") i
            write(equilfile, "(3A)") './equil_out/alpha_', trim(adjustl(str_i)),'.nc'
            write(*, *) 'Writing equilibration outfile to ', equilfile
            call write_qho_equilibration(pos_chain, energy_chain, equilfile, ierr, qho, &
                accept_rate, alpha_array(i))
 
            ! Write restart files, if required
            if (write_restart) then
                if (mod(i,restart_num)==0) then
                    ! Write res_nml.txt
                    call write_res_nml(i+1)
                endif
            endif
 
        end do
 
        ! Closing print statements
 
        write(*, "('--------------  End VQMC Runs  --------------')")
        if (write_log) then
            write(logid, "('--------------  End VQMC Runs  --------------')")
        end if
 
        !if (write_log) then
        !    write(logid, "('Results:')")
        !     write(logid, "('alpha            Energy')")
        !    do i = 1, size(alpha_array)
        !        write(logid, "(F12.8, '    ', F12.8)") alpha_array(i), energies(i)
        !    end do
        !    write(logid, "('')")
        !end if
    

Here is the call graph for this function:

Here is the caller graph for this function:

vqmc_h_2_plus()

subroutine vqmc::vqmc_h_2_plus	(	real(real64), intent(in)	sigma,
		integer, intent(in)	N,
		integer, intent(in)	burn,
		integer, intent(in)	thin,
		real(real64), intent(in)	c_const,
		real(real64), dimension(3), intent(in)	R_a,
		real(real64), dimension(3), intent(in)	R_b,
		real(real64), dimension(:, :), intent(out), allocatable	X_burnt,
		real(real64), intent(out)	accept_rate,
		real(real64), dimension(:), intent(out), allocatable	energy_chain,
		real(real64), intent(out)	energy_mean,
		real(real64), intent(out)	energy_variance,
		real(real64), dimension(:), intent(out), allocatable	DphiDc_chain,
		real(real64), dimension(3), optional	x_start
	)

Performs VQMC for the \( H_{2}^{+} \) system.

Computes the total energy of a trial wavefunction by averaging over a Markov chain of local energies generated by H_2_plus_Energy, with random move acceptance/rejection probability determined by H_2_plus_Prob. The H-H bond length can be changed by modifying the hydrogen atom positions R_a and R_b. Also outputs the gradient with respect to the variational parameter \( c \) at each step in the chain, to be used in Auto_H_2_plus.

Allows for modification of total Markov chain length, as well as the number of steps to be 'burned' while the chain is equilibrating and the interval to 'thin' the chain over to reduce correlation of samples. The size of the proposed random moves can be controlled by modifying the sigma parameter to get an ideal move acceptance rate.

Parameters

[in]	sigma	Step size
[in]	N	Number of steps to take
[in]	burn	Number of sample to discard at the beginning
[in]	thin	The frequency at which we discard the samples from the MCMC chain after burning
[in]	c_const	Parameter(s) of trial wavefunction
[in]	R_a	Position of hydrogen atom A
[in]	R_b	Position of hydrogen atom B
[in]	x_start	Starting point of MCMC chain (optional, default is randomly generated)
[out]	X_burnt	The burnt and thinned MCMC chain for the random walker
[out]	accept_rate	The acceptance rate of the MCMC chain
[out]	energy_chain	An array with the local energies of the each points in the MCMC chain
[out]	energy_mean	The average of the local energies recorded in energy_chain
[out]	energy_variance	The variance of the local energies recorded in energy_chain
[out]	DphiDc_chain	Array of the gradients with respect to \( c \) at each point in the MCMC chain

Definition at line 857 of file vqmc.f90.

 
        implicit none 
 
        real(real64), intent(in), dimension(3) :: R_a, R_b
        real(real64), intent(in) :: sigma
        real(real64), intent(in) :: c_const
        real(real64), dimension(3), optional :: x_start
        integer, intent(in) :: N, burn, thin 
        
         
        real(real64), intent(out) :: accept_rate, energy_mean, energy_variance
        real(real64), intent(out), dimension(:, :), allocatable :: X_burnt
        real(real64), intent(out), dimension(:), allocatable :: energy_chain
        real(real64), intent(out), dimension(:), allocatable :: DphiDc_chain
 
        real(real64) :: E_loc                               ! Local energy at current MCMC walker position
        real(real64) :: DphiDc                              ! Gradient wrt c at current MCMC walker position
        real(real64) :: log_prob_n                          ! The natural logarithm of the transition proability
        real(real64) :: a_prob_trans                        ! Acceptance ratio of the transitional probability 
        real(real64) :: unif_num                            ! A random sample from the uniform distribution
        real(real64), dimension(3) :: start                 ! Starting point of MCMC chain
        real(real64), dimension(:, :), allocatable :: X     ! Unburnt and thinned MCMC chain 
        real(real64), dimension(3) :: X_n                   ! The proposed step 
        real(real64), dimension(3) :: sample_normal         ! A sample from the normal distribution 
        integer :: accept                                   ! Number of accepted steps 
        integer :: thin_counter                             ! Counter to keep track of the thinning process
        integer :: reset_counter                            ! Counts how many times thin counter was reset
        integer :: re_size                                  ! The size of the burnt and thinned array
        integer :: i_bt                                     ! Burned/thinned value of i
        integer :: i                                        ! Loop variable 
        character(len=1) :: acceptance                      ! Yes/No signpost for acceptance
 
        if (i_log) then
            write(logid, "('Markov chain for c = ', F7.4)") c_const
            if (.not. run_equil) write(logid, "('Logging the chain from Rank 0.')")
        endif 
 
        ! Adjust the output array to be of the burnt and thinned array
        re_size = (n - burn) / thin
 
        allocate(x(3, 1:n+1)); x(:, :) = 0.0_real64
        allocate(x_burnt(3, 1:re_size)); x_burnt(:, :) = 0.0_real64
        allocate(energy_chain(1:re_size)); energy_chain(:) = 0.0_real64
        allocate(dphidc_chain(1:re_size)); dphidc_chain(:) = 0.0_real64
 
        ! Allow the starting position to take on the default value
        if (present(x_start)) then 
            start = x_start
            x(:,1) = start
            if (i_log) then
                write(logid, "('Starting from provided Markov chain position:')")
                write(logid, "('    x = ', F12.8, ' y = ', F12.8, ' z = ', F12.8)") &
                    x(1, 1), x(2, 1), x(3, 1)
            end if
        else 
            ! Create the initial position of the random walker that is within a certain radius
            call random_normal(3, sample_normal)
            do while(norm2(sample_normal) >= h2plus%rcut)
                call random_normal(3, sample_normal)
            end do 
            start = sample_normal
            x(:,1) = start 
            if (i_log) then
                write(logid, "('Generated random Markov chain start position:')")
                write(logid, "('X = ', F12.8, ', Y = ', F12.8, ', Z = ', F12.8)") &
                    x(1, 1), x(2, 1), x(3, 1)
            end if
        end if
        
        ! Start the counter for the number of accepted moves
        accept = 0  
        ! Start the counter for the thinning process 
        thin_counter = 0
        ! Start the reset_counter 
        reset_counter = 0
        
        if (i_log) then
            write(logid, "('Sample      X               Y              Z              &
                &Local Energy     Accepted?')")
            write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8, ' ', F12.6)")&
                1, x(1, 1), x(2, 1), x(3, 1), energy_chain(1)
        end if
        
        ! Calculate the local energy of the initial step if not burnt
        if (burn == 0) then 
            call h_2_plus_energy(c_const, x(:, 1), r_a, r_b, e_loc, dphidc)
            energy_chain(1) = e_loc
            dphidc_chain(1) = dphidc
            x_burnt(:, 1) = x(:, 1)
 
            do i = 2, n
                ! Add one to the thin counter 
                thin_counter = 1 + thin_counter
                ! Propose the next step 
                call random_normal(3, sample_normal)
                x_n = x(:, i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between proposed step and current step 
                log_prob_n = log(h_2_plus_prob(c_const, x(:, i-1), x_n, r_a, r_b))
           
                ! Calculate the acceptance ratio 
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
            
                ! Generate sample from uniform distribution
                call ran1(unif_num, seed)
 
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio
                if (a_prob_trans >= unif_num) then 
                    ! If accepted add proposed position to the MCMC chain
                    x(:, i) = x_n 
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(:, i_bt) = x_n
                        ! Calculate the local energy at this point 
                        call h_2_plus_energy(c_const, x_n, r_a, r_b, e_loc, dphidc)
                        energy_chain(i_bt) = e_loc
                        dphidc_chain(i_bt) = dphidc
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'Y'
                        ! Increment the number of accepted steps 
                        accept = 1 + accept 
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    &' ', F12.6, '         ', A1)") (i_bt), x_burnt(1, i_bt), &
                                    x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
                                    acceptance
                            end if
                        end if
                    end if 
                else
                    ! If rejected the proposed step is the current step 
                    x(:, i) = x(:, i-1)
 
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(:, i_bt) = x(:, i-1)
                        ! Calculate the local energy at this point 
                        call h_2_plus_energy(c_const, x(:, i-1), r_a, r_b, e_loc, dphidc)
                        energy_chain(i_bt) = e_loc
                        dphidc_chain(i_bt) = dphidc
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'N'
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    &' ', F12.6, '         ', A1)") (i_bt), x_burnt(1, i_bt), &
                                    x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
                                    acceptance
                            end if
                        end if
                    end if 
                end if 
            end do 
        ! Loop to advance the chain to the final burn step and not save any of the outputs
        else if (burn >= 2) then
            do i = 2, burn
                ! Propose the next step 
                call random_normal(3, sample_normal)
                x_n = x(:, i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between proposed step and current step 
                log_prob_n = log(h_2_plus_prob(c_const, x(:, i-1), x_n, r_a, r_b))
           
                ! Calculate the acceptance ratio 
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
            
                ! Generate sample from uniform distribution
                call ran1(unif_num, seed)
 
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio
                if (a_prob_trans >= unif_num) then 
                    ! If accepted add proposed position to the MCMC chain
                    x(:, i) = x_n 
                else 
                    ! If rejected the proposed step is the current step 
                    x(:, i) = x(:, i-1)
                end if 
            end do 
        end if
 
        ! Metropolis Algorithm to calculate the burnt and thinned MCMC chain 
        if (burn >= 1) then 
            ! Metropolis algorithm 
            do i = burn + 1, n+1 
                ! Add one to the thin counter 
                thin_counter = 1 + thin_counter
                ! Propose the next step 
                call random_normal(3, sample_normal)
                x_n = x(:, i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between proposed step and current step 
                log_prob_n = log(h_2_plus_prob(c_const, x(:, i-1), x_n, r_a, r_b))
                
                ! Calculate the acceptance ratio 
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
            
                ! Generate sample from uniform distribution
                call ran1(unif_num, seed)
 
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio
                if (a_prob_trans >= unif_num) then 
                    ! If accepted add proposed position to the MCMC chain
                    x(:, i) = x_n 
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(:, i_bt) = x_n
                        ! Calculate the local energy at this point 
                        call h_2_plus_energy(c_const, x_n, r_a, r_b, e_loc, dphidc)
                        energy_chain(i_bt) = e_loc
                        dphidc_chain(i_bt) = dphidc
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'Y'
                        ! Increment the number of accepted steps 
                        accept = 1 + accept 
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    &' ', F12.6, '         ', A1)") (i_bt), x_burnt(1, i_bt), &
                                    x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
                                    acceptance
                            end if
                        end if
                    end if 
                 
                else
                    ! If rejected the proposed step is the current step 
                    x(:, i) = x(:, i-1)
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain                   
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_burnt(:, i_bt) = x(:, i-1)
                        ! Calculate the local energy at this point 
                        call h_2_plus_energy(c_const, x(:, i-1), r_a, r_b, e_loc, dphidc)
                        energy_chain(i_bt) = e_loc
                        dphidc_chain(i_bt) = dphidc
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'N'
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    &' ', F12.6, '         ', A1)") (i_bt), x_burnt(1, i_bt), &
                                    x_burnt(2, i_bt), x_burnt(2, i_bt), energy_chain(i_bt), &
                                    acceptance
                            end if
                        end if
 
                    end if 
                end if 
            end do
        end if
 
        ! Calculate the acceptance ratio
        accept_rate = (real(accept, kind=real64) / real(re_size, kind=real64))
 
        ! Calculate the mean energy and the variance in the energy
        energy_mean = sum(energy_chain) / real(re_size, kind=real64)
        energy_variance = 0.0_real64
        do i = 1, re_size
            energy_variance = energy_variance + ((1.0_real64 / real(re_size - 1, kind=real64)) * &
            (energy_chain(i) - energy_mean) * (energy_chain(i) - energy_mean)) 
        end do
        energy_variance = (1.0_real64 / real(re_size, kind=real64)) * energy_variance
 
        if (i_log) then
            write(logid, "('Local energy = ', F12.8)") energy_mean
            write(logid, "('Variance = ', F12.8)") energy_variance
            write(logid, "('Acceptance rate = ', F12.8)") accept_rate
            write(logid, "('')")
        end if
 

Here is the call graph for this function:

Here is the caller graph for this function:

auto_h_2_plus()

subroutine vqmc::auto_h_2_plus	(	real(real64), intent(in)	c_in,
		real(real64), intent(in)	bond_length,
		real(real64), intent(out)	c_out,
		real(real64), intent(out)	energy_out,
		real(real64), intent(out)	variance_out,
		real(real64), intent(out)	accept_out,
		character(len=1), intent(out)	converged
	)

Automatically finds the best variational parameter \( c \) for the \( H_{2}^{+} \) problem at given bond length.

Determines the optimum value of \( c \) at a given bond length from a guess value by utilising the gradient with respect to \( c \) (calculated by first order central difference) in a dampened steepest descent optimisation.

Parameters

[in]	c_in	Guess value for c
[in]	bond_length	H-H bond length
[out]	c_out	Optimised value for c at this bond length
[out]	energy_out	Energy of wavefunction at optimal c
[out]	variance_out	Variance of energy at optimal c
[out]	accept_out	Acceptance rate of MCMC chain at optimal c
[out]	converged	Y/N convergence status

Definition at line 1188 of file vqmc.f90.

 
        implicit none
 
        real(real64), intent(in) :: c_in
        real(real64), intent(in) :: bond_length
        real(real64), intent(out) :: c_out
        real(real64), intent(out) :: energy_out
        real(real64), intent(out) :: variance_out
        real(real64), intent(out) :: accept_out
        character(len=1), intent(out) :: converged
 
        real(real64), dimension(3) :: R_a, R_b
        real(real64) :: c_current, c_new
        real(real64) :: c_grad
        real(real64), dimension(:, :), allocatable :: pos_chain
        real(real64), dimension(:), allocatable :: energy_chain
        real(real64), dimension(:), allocatable :: DphiDc_chain
        integer :: viz_steps
        integer :: k
        integer :: ierr
        integer :: maxiter = 200
        real(real64) :: tol = 1e-4_real64
        character(len=6) :: str_b
        character(len=40) :: chainfile
 
        ! MPI Variables
        real(real64) :: my_c_grad, my_ensum
        real(real64) :: sum_var, sum_acc
 
        ! Initialise
        c_current = c_in
        r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
        r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
 
        if (i_log) then
            write(logid, "('Starting parameter search for c...')")
        end if
 
        ! Main iteration loop
        do k = auto_start, maxiter
 
            ! Calculate gradient at current c
            call vqmc_h_2_plus(h2plus%sigma, my_steps+h2plus%burn_step, h2plus%burn_step, h2plus%thin_step, &
                    c_current, r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, & 
                    variance_out, dphidc_chain)
            
            ! Find gradient     
            my_c_grad = 2.0_real64 * (sum(energy_chain * dphidc_chain) / size(energy_chain)) &
                - (energy_out * (sum(dphidc_chain) / size(dphidc_chain)))
 
            ! Find average c_grad across ranks, on all ranks.
            call mpi_allreduce(my_c_grad,c_grad,1,mpi_double_precision,mpi_sum,mpi_comm_world,mpierr)
            c_grad = c_grad * per_proc
 
            ! Update the log
            if (i_log) then
                write(logid, "('c            Gradient')")
                write(logid, "(F13.7, '       ', F13.7)") c_current, c_grad
            end if
 
            ! Check if gradient is within tolerance.
            if (abs(c_grad) .le. tol) then
 
                ! Update c_out, and note convergence, on all ranks.
                c_out = c_current
                converged = 'Y'
 
                ! Sum the chain energy-sums and store in rank 0
                my_ensum = sum(energy_chain)
                call mpi_reduce(my_ensum,energy_out,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
                ! Sum the chain variances and store in rank 0
                call mpi_reduce(variance_out, sum_var,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
                ! Sum the chain acceptance rates and store in rank 0
                call mpi_reduce(accept_out, sum_acc,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
                ! Calculate averages on rank 0
                if (my_rank==0) then
 
                    energy_out = energy_out * per_step
                    variance_out = sum_var * per_proc2 
                    accept_out = sum_acc * per_proc
 
                    ! Update the log
                    if (write_log) then
                        write(logid, "('Parameter search converged!')")
                        write(logid, "('Results:')")
                        write(logid, "('c = ', F13.7, ', Energy = ', F12.8, ', Variance = ', F12.8)") &
                            c_current, energy_out, variance_out
                        write(logid, "('')")
                    end if
 
                end if 
 
                ! Exit subroutine, as convergence is reached.
                ! Note: Only rank 0 has the correct array outputs. All ranks have correct c_out and 'converged'.
                exit
 
            end if
 
            ! If not, update c by steepest descent, on all ranks.
            call h_2_plus_update_c(c_current, c_grad, c_new)
            c_current = c_new
 
            ! Write restart files, if required
            if (write_restart .and. my_rank==0) then
                if (mod(k,restart_num)==0) then
                    call write_res_nml(i_bond,k+1,h2plus%auto_params,c_current)
                endif
            endif 
 
        end do
 
        ! Deal with non-convergence
        if (k .ge. maxiter) then
            converged = 'N'
            c_out = c_current
            if (i_log) then 
                write(logid, "('Parameter search NOT converged! Consider changing&
                    & the starting value of c, or the damping in the H_2_update_c&
                    & function.')")
            end if
        end if
 
        ! Report results
        if (my_rank==0) then
 
            write(*, "(F12.8, '   ', A, '             ', F13.7, '      ', F12.8)") &
                bond_length, converged, c_out, energy_out
 
            ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
            if (write_chains) then
                write(str_b, "(F6.3)") bond_length
                write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
                write(*, "('Running VQMC at c = ', F12.7, ' to produce single Markov chain...')") c_out
                viz_steps = min(h2plus%steps,max_chain_len)
                call vqmc_h_2_plus(h2plus%sigma, viz_steps, h2plus%burn_step, h2plus%thin_step, &
                    c_out, r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, & 
                    variance_out, dphidc_chain)
                write(*, "(2A)") 'Writing Markov chain to ', chainfile
                call write_h2plus_equilibration(pos_chain, energy_chain, chainfile, ierr, h2plus, &
                    accept_out, c_out, bond_length)
            end if
        end if
        

Here is the call graph for this function:

Here is the caller graph for this function:

grid_h_2_plus()

subroutine vqmc::grid_h_2_plus	(	real(real64), dimension(:), intent(in)	c_array,
		real(real64), intent(in)	bond_length,
		real(real64), intent(out)	c_out,
		real(real64), intent(out)	energy_out,
		real(real64), intent(out)	variance_out,
		real(real64), intent(out)	accept_out
	)

Performs a grid search to find the optimal value for the variational parameter \( c \) in the \( H_{2}^{+} \) problem.

Evaluates energy of the \( H_{2}^{+} \) trial wavefunction at each value of the parameter \( c \) in c_array. Determines the best energy, and reports its corresponding \( c \) as the optimal value at this bond length. Also returns the variance of this energy as a measure of uncertainty, as well as the acceptance rate of the MCMC chain.

Parameters

[in]	c_array	Array of c values to be tested
[in]	bond_length	H-H bond length
[out]	c_out	Optimised value for c at this bond length
[out]	energy_out	Energy of wavefunction at optimal c
[out]	variance_out	Variance of energy at optimal c
[out]	accept_out	Acceptance rate of MCMC chain at optimal c

Definition at line 1359 of file vqmc.f90.

    
        implicit none
 
        real(real64), dimension(:), intent(in) :: c_array
        real(real64), intent(in) :: bond_length
        real(real64), intent(out) :: c_out
        real(real64), intent(out) :: energy_out
        real(real64), intent(out) :: variance_out
        real(real64), intent(out) :: accept_out
 
        real(real64), dimension(3) :: R_a, R_b
        real(real64), dimension(:, :), allocatable :: pos_chain
        real(real64), dimension(:), allocatable :: energy_chain, DphiDc_chain
        real(real64), dimension(:), allocatable :: param_energies, param_accepts, param_variances
        integer :: viz_steps
        integer :: i
        integer :: ierr
        character(len=6) :: str_b
        character(len=40) :: chainfile
 
        ! MPI variables
        real(real64)  :: my_ensum
 
        ! Initialise
        r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
        r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
        allocate(param_energies(size(c_array))); param_energies(:)= 0.0_real64
        allocate(param_accepts(size(c_array))); param_accepts(:) = 0.0_real64
        allocate(param_variances(size(c_array))); param_variances(:) = 0.0_real64
 
        ! Sort out restarting
        if (run_restart) then
            ! Read in array values
            call read_energies(resfile_epar, param_energies, param_variances, param_accepts)
            ! Reset the tag
            run_restart = .false.
        endif 
 
        if (i_log) then
            write(logid, "('Starting parameter search for c...')")
        end if
 
        ! Loop over c values
        do i = grid_start, size(c_array)
 
            ! Calculate energy at current value of c
            call vqmc_h_2_plus(h2plus%sigma, my_steps+h2plus%burn_step, h2plus%burn_step, h2plus%thin_step, &
                c_array(i), r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, &
                variance_out, dphidc_chain)
 
            ! Save energy - serial version
            !param_energies(i) = energy
            !param_variances(i) = variance
            !param_accepts(i) = accept_rate
 
            ! Sum the chain energy-sums and store in rank 0
            my_ensum = sum(energy_chain)
            call mpi_reduce(my_ensum,param_energies(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Sum the chain variances and store in rank 0
            call mpi_reduce(variance_out,param_variances(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Sum the chain acceptance rates and store in rank 0
            call mpi_reduce(accept_out,param_accepts(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Calculate averages on rank 0
            if (my_rank==0) then
                param_energies(i) = param_energies(i) * per_step
                param_variances(i) = param_variances(i) * per_proc2 
                param_accepts(i) = param_accepts(i) * per_proc
 
                ! Write restart files, if required
                if (write_restart) then
                    if (mod(i,restart_num)==0) then
                        ! Write res_nml.txt
                        call write_res_nml(i_bond,i+1,h2plus%auto_params)
                        ! Write param energies, variances and acceptance rates
                        call write_qho_main(c_array, param_energies, param_variances, resfile_epar, ierr, qho, param_accepts)
                    endif
                endif 
            endif 
 
        end do
 
        ! Compile final results of grid search on rank 0
        ! Note: Only rank 0 has the correct _out values.
        if (my_rank==0) then
 
            ! Select lowest energy and corresponding c value
            energy_out = minval(param_energies)
            c_out = c_array(minloc(param_energies, 1))
            variance_out = param_variances(minloc(param_energies, 1))
            accept_out = param_accepts(minloc(param_energies, 1))
 
            write(*, "(F12.8, '   -             ', F13.7, '      ', F12.8)") &
                bond_length, c_out, energy_out
 
            if (write_log) then
                write(logid, "('Results:')")
                write(logid, "('c                Energy')")
                do i = 1, size(c_array)
                    write(logid, "(F12.8, '    ', F12.8)") c_array(i), param_energies(i)
                end do
                write(logid, "('')")
                write(logid, "('Min Energy = ', F12.8, ' at c = ', F12.8)") energy_out, c_out
                write(logid, "('')")
            end if
 
            ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
            if (write_chains) then
                write(str_b, "(F6.3)") bond_length
                write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
                write(*, "('Running VQMC at c = ', F12.7, ' to produce single Markov chain...')") c_out
                viz_steps = min(h2plus%steps,max_chain_len)
                call vqmc_h_2_plus(h2plus%sigma, viz_steps, h2plus%burn_step, h2plus%thin_step, &
                    c_out, r_a, r_b, pos_chain, accept_out, energy_chain, energy_out, & 
                    variance_out, dphidc_chain)
                write(*, "(2A)") 'Writing Markov chain to ', chainfile
                call write_h2plus_equilibration(pos_chain, energy_chain, chainfile, ierr, h2plus, &
                    accept_out, c_out, bond_length)
            end if
 
        endif 
    

Here is the call graph for this function:

Here is the caller graph for this function:

equil_h_2_plus()

subroutine vqmc::equil_h_2_plus	(	real(real64), dimension(:), intent(in)	c_array,
		real(real64), intent(in)	bond_length,
		integer, intent(in)	bonds_pos
	)

Performs VQMC on an array of values for the variational parameter \( c \), outputting energy traces for each.

Generates VQMC position and energy traces of the \( H_{2}^{+} \) trial wavefunction at each value of the parameter \( c \) in c_array. Outputs each of these traces to a separate NetCDF file, to be used in determining optimal burning/thinning paramters for an accurate full VQMC run.

Parameters

[in]	c_array	The array of c values to obtain energy traces for
[in]	bond_length	The current H-H bond length
[in]	bonds_pos	The position in the array of bond lengths, to be used in designating output files

Definition at line 1505 of file vqmc.f90.

 
        implicit none
 
        real(real64), dimension(:), intent(in) :: c_array
        real(real64), intent(in) :: bond_length
        integer, intent(in) :: bonds_pos
 
        real(real64), dimension(3) :: R_a, R_b
        real(real64), dimension(:), allocatable :: variances
        real(real64), dimension(:), allocatable :: energies
        real(real64), dimension(:), allocatable :: energy_chain
        real(real64), dimension(:, :), allocatable :: pos_chain
        real(real64), dimension(:), allocatable :: DphiDc_chain
        real(real64) :: accept_rate
        real(real64) :: energy_variance
        integer :: i, ierr
        character(len=3) :: str_i, str_b
        character(len=40) :: equilfile
 
        ! Initialise
        r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
        r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
        allocate(energies(size(c_array))); energies(:) = 0.0_real64
        allocate(variances(size(c_array)))
 
        ! Loop over values of c
        do i = grid_start, size(c_array)
 
            ! Calculate energy at current value of c
            call vqmc_h_2_plus(h2plus%sigma, h2plus%steps, h2plus%burn_step, h2plus%thin_step, &
                c_array(i), r_a, r_b, pos_chain, accept_rate, energy_chain,energies(i), &
                energy_variance, dphidc_chain)
 
            write(*, "(F12.8, '   -             ', F13.7, '      ', F12.8)") &
                bond_length, c_array(i), energies(i)
 
            ! Save the energy trace for each value of c
            write(str_i, "(I3)") i
            write(str_b, "(I3)") bonds_pos
            write(equilfile, "(5A)") './equil_out/b', trim(adjustl(str_b)), &
                '_c', trim(adjustl(str_i)), '.nc'
            write(*, *) 'Writing equilibration outfile to ', equilfile
            call write_h2plus_equilibration(pos_chain, energy_chain, equilfile, ierr, h2plus, &
                accept_rate, c_array(i), bond_length)
 
            ! Write restart files, if required
            if (write_restart) then
                if (mod(i,restart_num)==0) then
                    call write_res_nml(i_bond,i+1,h2%auto_params)
                endif
            endif
 
        end do
 
        if (write_log) then
            write(logid, "('Results:')")
            write(logid, "('c                Energy')")
            do i = 1, size(c_array)
                write(logid, "(F12.8, '    ', F12.8)") c_array(i), energies(i)
            end do
            write(logid, "('')")
        end if
 
        deallocate(energies)
        deallocate(variances)
    

Here is the call graph for this function:

Here is the caller graph for this function:

vqmc_h_2()

subroutine vqmc::vqmc_h_2	(	real(real64), intent(in)	sigma,
		integer, intent(in)	N,
		integer, intent(in)	burn,
		integer, intent(in)	thin,
		real(real64), intent(in)	a,
		real(real64), intent(in)	beta,
		real(real64), dimension(3), intent(in)	R_a,
		real(real64), dimension(3), intent(in)	R_b,
		real(real64), dimension(:, :), intent(out), allocatable	X_1_burnt,
		real(real64), dimension(:, :), intent(out), allocatable	X_2_burnt,
		real(real64), intent(out)	accept_rate,
		real(real64), dimension(:), intent(out), allocatable	energy_chain,
		real(real64), intent(out)	energy_mean,
		real(real64), intent(out)	energy_variance,
		real(real64), dimension(:), intent(out), allocatable	DphiDbeta_chain,
		real(real64), dimension(3), optional	x_start_1,
		real(real64), dimension(3), optional	x_start_2
	)

Performs VQMC for the \( H_{2}^{+} \) system.

Computes the total energy of a trial wavefunction by averaging over a Markov chain of local energies generated by H_2_Energy, with random move acceptance/rejection probability determined by H_2_Prob. The H-H bond length can be changed by modifying the hydrogen atom positions R_a and R_b. Also outputs the gradient with respect to the variational parameter \( \beta \) at each step in the chain, to be used in Auto_H_2.

Allows for modification of total Markov chain length, as well as the number of steps to be 'burned' while the chain is equilibrating and the interval to 'thin' the chain over to reduce correlation of samples. The size of the proposed random moves can be controlled by modifying the sigma parameter to get an ideal move acceptance rate.

Parameters

[in]	sigma	Step size
[in]	N	Number of steps to take
[in]	burn	Number of sample to discard at the beginning
[in]	thin	The frequency at which we discard the samples from the MCMC chain after burning
[in]	a	Parameter a of trial wavefunction
[in]	beta	Parameter beta of trial wavefunction
[in]	R_a	Position of hydrogen atom A
[in]	R_b	Position of hydrogen atom B
[in]	x_start_1	Starting point of 1st MCMC chain (optional, default is randomly generated)
[in]	x_start_2	Starting point of 2nd MCMC chain (optional, default is randomly generated)
[out]	X_1_burnt	The burnt and thinned MCMC chain for random walker 1
[out]	X_2_burnt	The burnt and thinned MCMC chain for random walker 2
[out]	accept_rate	The acceptance rate of the MCMC chain
[out]	energy_chain	An array with the local energies of the each points in the MCMC chains
[out]	energy_mean	The average of the local energies recorded in energy_chain
[out]	energy_variance	The variance of the local energies recorded in energy_chain
[out]	DphiDbeta_chain	Array of the gradients with respect to \( \beta \) at each point in the MCMC chain

Definition at line 1617 of file vqmc.f90.

 
        implicit none 
 
        real(real64), intent(in), dimension(3) :: R_a, R_b
        real(real64), intent(in) :: sigma
        real(real64), intent(in) :: a, beta
        real(real64), dimension(3), optional :: x_start_1, x_start_2
        integer, intent(in) :: N, burn, thin 
        
         
        real(real64), intent(out) :: accept_rate, energy_mean, energy_variance
        real(real64), intent(out), dimension(:, :), allocatable :: X_1_burnt, X_2_burnt
        real(real64), intent(out), dimension(:), allocatable :: energy_chain
        real(real64), intent(out), dimension(:), allocatable :: DphiDbeta_chain
 
        real(real64) :: E_loc                                    ! Local energy at current MCMC walker position
        real(real64) :: DphiDbeta                                ! Gradient wrt beta at current MCMC walker position
        real(real64) :: log_prob_n                               ! The natural logarithm of the transition proability 
        real(real64) :: a_prob_trans                             ! Acceptance ratio of the transitional probability                   
        real(real64) :: unif_num                                 ! A random sample from the uniform distribution
        real(real64), dimension(3) :: start_1, start_2           ! Starting points for the random walkers 
        real(real64), dimension(3) :: X_n_1, X_n_2               ! The propsed steps for the random walkers 
        real(real64), dimension(3) :: sample_normal              ! Sample from the normal distribution 
        real(real64), dimension(:, :), allocatable :: X_1, X_2   ! Full MCMC chains 
        integer :: accept                                        ! Number of accepted steps 
        integer :: thin_counter                                  ! Counter to keep track of the thinning process
        integer :: reset_counter                                 ! Counts how many times thin counter was reset
        integer :: re_size                                       ! The size of the burnt and thinned array
        integer :: i_bt                                          ! Burned/thinned value of i
        integer :: i                                             ! Loop variable 
        character(len=1) :: acceptance                           ! Yes/No signpost for acceptance
        
        if (i_log) then
            write(logid, "('Markov chain for beta = ', F7.4)") beta
            if (.not. run_equil) write(logid, "('Logging the chain from Rank 0.')")
        endif 
 
        ! Adjust the output array to be of the burnt and thinned array
        re_size = (n - burn) / thin 
 
        allocate(x_1(3, 1:n+1)); x_1(:, :) = 0.0_real64
        allocate(x_2(3, 1:n+1)); x_2(:, :) = 0.0_real64
        allocate(x_1_burnt(3, 1:re_size)); x_1_burnt(:, :) = 0.0_real64
        allocate(x_2_burnt(3, 1:re_size)); x_2_burnt(:, :) = 0.0_real64
        allocate(energy_chain(1:re_size)); energy_chain(:) = 0.0_real64
        allocate(dphidbeta_chain(1:re_size)); dphidbeta_chain(:) = 0.0_real64
 
        ! Allow the starting position to take on the default value 
        if (present(x_start_1) .and. present(x_start_2)) then
            start_1 = x_start_1
            x_1(:,1) = start_1
            
            start_2 = x_start_2
            x_2(:,1) = start_2
 
            if (i_log) then
                write(logid, "('Starting from provided Markov chain position:')")
                write(logid, "('    x1 = ', F12.8, ' y1 = ', F12.8, ' z1 = ', F12.8)") &
                    x_1(1, 1), x_1(2, 1), x_1(3, 1)
                write(logid, "('    x2 = ', F12.8, ' y2 = ', F12.8, ' z2 = ', F12.8)") &
                    x_2(1, 1), x_2(2, 1), x_2(3, 1)
            end if
        else 
            ! Create the initial position of random walker 1 that is within a certain radius
            call random_normal(3, sample_normal)
            do while(norm2(sample_normal) >= h2%rcut)
                call random_normal(3, sample_normal)
            end do 
            start_1 = sample_normal
            x_1(:,1) = start_1
 
            ! Create the initial position of random walker 2 that is within a certain radius
            call random_normal(3, sample_normal)
            do while(norm2(sample_normal) >= h2%rcut)
                call random_normal(3, sample_normal)
            end do 
            start_2 = sample_normal 
            x_2(:,1) = start_2
        end if 
        
        ! Start the counter for the number of accepted moves
        accept = 0  
        ! Start the counter for the thinning process
        thin_counter = 0
        ! Start the reset counter 
        reset_counter = 0
 
        if (i_log) then
            write(logid, "('Sample      X1              Y1             Z1             &
                &X2              Y2             Z2             &
                &Local Energy     Accepted?')")
            write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8, '   ', F12.8, &
                &'   ', F12.8, '   ', F12.8, ' ', F12.6)") 1, x_1(1, 1), x_1(2, 1), x_1(3, 1),&
                x_2(1, 1), x_2(2, 1), x_2(3, 1), energy_chain(1)
        end if
 
        ! Calculate the local energy of the initial step if not burnt
        if (burn == 0) then 
            call h_2_energy(a, beta, x_1(:,1), x_2(:,1), r_a, r_b, e_loc, dphidbeta)
            energy_chain(1) = e_loc
            dphidbeta_chain(1) = dphidbeta
            x_1_burnt(:, 1) = x_1(:, 1)
            x_2_burnt(:, 1) = x_2(:, 1)
 
            do i = 2, n 
                ! Initialize thin counter 
                thin_counter = 1 + thin_counter
 
                ! Propose the next step for random walker 1 
                call random_normal(3, sample_normal)
                x_n_1 = x_1(:, i-1) + sigma * sample_normal
 
                ! Propose the next step for random walker 2 
                call random_normal(3, sample_normal)
                x_n_2 = x_2(:, i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between the proposed steps and current steps
                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))
 
                ! Calculate the acceptance ratio
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
 
                ! Generate sample from uniform distribution
                call ran1(unif_num, seed)
 
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio
                if (a_prob_trans >= unif_num) then 
                    ! If accepted calculate local energy and add proposed position to the MCMC chain
                    x_1(:, i) = x_n_1
                    x_2(:, i) = x_n_2 
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain                                        
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_1_burnt(:, i_bt) = x_n_1
                        x_2_burnt(:, i_bt) = x_n_2
                        ! Calculate the local energy at this point 
                        call h_2_energy(a, beta, x_n_1, x_n_2, r_a, r_b, e_loc, dphidbeta)
                        energy_chain(i_bt) = e_loc
                        dphidbeta_chain(i_bt) = dphidbeta
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        ! Increment the number of accepted steps 
                        accept = 1 + accept 
                        acceptance = 'Y'
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    & '   ', F12.8, '   ', F12.8, '   ', F12.8, ' ', F12.6, '         ', A1)") &
                                    (i-burn), x_1_burnt(1, i-burn), x_1_burnt(2, i-burn), x_1_burnt(3, i-burn),&
                                    x_2_burnt(1, i-burn), x_2_burnt(2, i-burn), x_2_burnt(3, i-burn),&
                                    energy_chain(i-burn), acceptance
                            end if
                        end if
                    end if 
                else 
                    ! If rejected the proposed step is the current step 
                    x_1(:, i) = x_1(:, i-1)
                    x_2(:, i) = x_2(:, i-1) 
                
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_1_burnt(:, i_bt) = x_1(:, i-1)
                        x_2_burnt(:, i_bt) = x_2(:, i-1)
                        ! Calculate the local energy at this point 
                        call h_2_energy(a, beta, x_1(:, i-1), x_2(:, i-1), r_a, r_b, e_loc, dphidbeta)
                        energy_chain(i_bt) = e_loc
                        dphidbeta_chain(i_bt) = dphidbeta
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        acceptance = 'N'
 
                        if (i_log) then
                            if (mod(i_bt, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    & '   ', F12.8, '   ', F12.8, '   ', F12.8, ' ', F12.6, '         ', A1)") &
                                    (i_bt), x_1_burnt(1, i_bt), x_1_burnt(2, i_bt), x_1_burnt(3, i_bt),&
                                    x_2_burnt(1, i_bt), x_2_burnt(2, i_bt), x_2_burnt(3, i_bt),&
                                    energy_chain(i_bt), acceptance
                            end if
                        end if
                    end if 
                end if 
            end do 
        ! Loop to advance the chain to the final burn step and not save any of the outputs
        else if (burn >= 2) then 
            do i = 2, burn 
                ! Propose the next step for random walker 1 
                call random_normal(3, sample_normal)
                x_n_1 = x_1(:, i-1) + sigma * sample_normal
 
                ! Propose the next step for random walker 2 
                call random_normal(3, sample_normal)
                x_n_2 = x_2(:, i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between the proposed steps and current steps
                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))
 
                ! Calculate the acceptance ratio
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
 
                ! Generate sample from uniform distribution
                call ran1(unif_num, seed)
 
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio
                if (a_prob_trans >= unif_num) then 
                    ! If accepted calculate local energy and add proposed position to the MCMC chain
                    x_1(:, i) = x_n_1
                    x_2(:, i) = x_n_2 
                else 
                    ! If rejected the proposed step is the current step 
                    x_1(:, i) = x_1(:, i-1)
                    x_2(:, i) = x_2(:, i-1) 
                end if 
            end do 
 
            if (i_log) write(logid, "('---------- Burning steps completed. ----------')")
 
        end if
 
        ! Metropolis Algorithm to calculate the burnt and thinned MCMC chain 
        if (burn >= 1) then 
            ! Metropolis algorithm 
            do i = (burn + 1) , n+1 
                ! Initialize thin counter 
                thin_counter = 1 + thin_counter
                ! Propose the next step for random walker 1 
                call random_normal(3, sample_normal)
                x_n_1 = x_1(:, i-1) + sigma * sample_normal
 
                ! Propose the next step for random walker 2 
                call random_normal(3, sample_normal)
                x_n_2 = x_2(:, i-1) + sigma * sample_normal
 
                ! Calculate the transition probability between the proposed steps and current steps
                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))
 
                ! Calculate the acceptance ratio
                a_prob_trans = min(1.0_real64, exp(log_prob_n))
 
                ! Generate sample from uniform distribution
                call ran1(unif_num, seed)
 
                ! Compare the drawn sample from the uniform distribution with the acceptance ratio
                if (a_prob_trans >= unif_num) then 
                    ! If accepted calculate local energy and add proposed position to the MCMC chain
                    x_1(:, i) = x_n_1
                    x_2(:, i) = x_n_2 
 
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain                    
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_1_burnt(:, i_bt) = x_n_1
                        x_2_burnt(:, i_bt) = x_n_2
                        ! Calculate the local energy at this point 
                        call h_2_energy(a, beta, x_n_1, x_n_2, r_a, r_b, e_loc, dphidbeta)
                        energy_chain(i_bt) = e_loc
                        dphidbeta_chain(i_bt) = dphidbeta
                        ! Resetting the thin counter 
                        thin_counter = 0
 
                        ! Increment the number of accepted steps 
                        accept = 1 + accept 
                        acceptance = 'Y'
 
                        if (i_log) then
                            if (mod(i-burn, re_size/logfrac) .le. 1e-6_real64) then
                                write(logid, "(I8, '   ', F12.8, '   ', F12.8, '   ', F12.8,&
                                    & '   ', F12.8, '   ', F12.8, '   ', F12.8, ' ', F12.6, '         ', A1)") &
                                    (i_bt), x_1_burnt(1, i_bt), x_1_burnt(2, i_bt), x_1_burnt(3, i_bt),&
                                    x_2_burnt(1, i_bt), x_2_burnt(2, i_bt), x_2_burnt(3, i_bt),&
                                    energy_chain(i_bt), acceptance
                            end if
                        end if
                    end if 
                else 
                    ! If rejected the proposed step is the current step 
                    x_1(:, i) = x_1(:, i-1)
                    x_2(:, i) = x_2(:, i-1) 
                
                    ! To determine whether this point should be added to the burnt and thinned MCMC chain                
                    if ((i > (burn + 1)) .and. (thin_counter == thin)) then
                        ! Increment the reset counter 
                        reset_counter = 1 + reset_counter
                        ! Precalculate i_bt for shorthand
                        i_bt = i - burn - (reset_counter * (thin - 1))
                        ! Add to the burn and thinned MCMC chain
                        x_1_burnt(:, i_bt) = x_1(:, i-1)
                        x_2_burnt(:, i_bt) = x_2(:, i-1)
                        ! Calculate the local energy at this point 
                        call h_2_energy(a, beta, x_1(:, i-1), x_2(:, i-1), r_a, r_b, e_loc, dphidbeta)
                        energy_chain(i_bt) = e_loc
                        dphidbeta_chain(i_bt) = dphidbeta
                        ! Resetting the thin counter 
                        thin_counter = 0
                    end if 
                end if 
            end do 
        end if
 
        ! Calculate the acceptance ratio
        accept_rate = (real(accept, kind=real64) / real(re_size, kind=real64))
 
        ! Calculate the mean energy and the variance in the energy
        energy_mean = sum(energy_chain) / real(re_size, kind=real64)
        energy_variance = 0.0_real64
        do i = 1, re_size
            energy_variance = energy_variance + ((1.0_real64 / real(re_size - 1, kind=real64)) * &
            (energy_chain(i) - energy_mean) * (energy_chain(i) - energy_mean)) 
        end do
        energy_variance = (1.0_real64 / real(re_size, kind=real64)) * energy_variance
 
        if (i_log) then
            write(logid, "('Local energy = ', F12.8)") energy_mean
            write(logid, "('Variance = ', F12.8)") energy_variance
            write(logid, "('Acceptance rate = ', F12.8)") accept_rate
            write(logid, "('')")
        end if
 

Here is the call graph for this function:

Here is the caller graph for this function:

auto_h_2()

subroutine vqmc::auto_h_2	(	real(real64), intent(in)	a,
		real(real64), intent(in)	beta_in,
		real(real64), intent(in)	bond_length,
		real(real64), intent(out)	beta_out,
		real(real64), intent(out)	energy_out,
		real(real64), intent(out)	variance_out,
		real(real64), intent(out)	accept_out,
		character(len=1), intent(out)	converged
	)

Automatically finds the best variational parameter \( \beta \) for the \( H_{2}^{+} \) problem at a given bond length.

Determines the optimum value of \( \beta \) at a given bond length from a guess value by utilising the gradient with respect to \( \beta \) (calculated by first order central difference) in a dampened steepest descent optimisation.

Parameters

[in]	a	Value of parameter a at current bond length, calculated by H_2_cusp
[in]	beta_in	Guess value for beta
[in]	bond_length	H-H bond length
[out]	beta_out	Optimised value for beta at this bond length
[out]	energy_out	Energy of wavefunction at optimal beta
[out]	variance_out	Variance of energy at optimal beta
[out]	accept_out	Acceptance rate of MCMC chain at optimal beta
[out]	converged	Y/N convergence status

Definition at line 1980 of file vqmc.f90.

 
        implicit none
 
        real(real64), intent(in) :: a
        real(real64), intent(in) :: beta_in
        real(real64), intent(in) :: bond_length
        real(real64), intent(out) :: beta_out
        real(real64), intent(out) :: energy_out
        real(real64), intent(out) :: variance_out
        real(real64), intent(out) :: accept_out
        character(len=1), intent(out) :: converged
 
        real(real64), dimension(3) :: R_a, R_b
        real(real64) :: beta_current, beta_new
        real(real64) :: beta_grad
        real(real64), dimension(:, :), allocatable :: e1_pos_chain, e2_pos_chain
        real(real64), dimension(:), allocatable :: energy_chain, DphiDbeta_chain
        integer :: viz_steps
        integer :: k
        integer :: ierr
        integer :: maxiter = 200
        real(real64) :: tol = 1e-4_real64
        character(len=6) :: str_b
        character(len=40) :: chainfile
 
        ! MPI Variables
        real(real64) :: my_beta_grad, my_ensum
        real(real64) :: sum_var, sum_acc
 
        ! Initialise
        beta_current = beta_in
        r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
        r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
 
        if (i_log) then
            write(logid, "('Starting parameter search for beta...')")
        end if
 
        ! Main iteration loop
        do k = auto_start, maxiter
 
            ! Calculate gradient at current beta
            call vqmc_h_2(h2%sigma, my_steps+h2%burn_step, h2%burn_step, h2%thin_step, a, &
                    beta_current, r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
                    energy_chain, energy_out, variance_out, dphidbeta_chain)
 
            ! Find gradient 
            my_beta_grad = 2.0_real64 * ((sum(energy_chain * dphidbeta_chain) / size(energy_chain)) &
                - energy_out * (sum(dphidbeta_chain) / size(energy_chain))) 
 
            ! Find average beta_grad across ranks, on all ranks.
            call mpi_allreduce(my_beta_grad,beta_grad,1,mpi_double_precision,mpi_sum,mpi_comm_world,mpierr)
            beta_grad = beta_grad * per_proc
 
            ! Update the log
            if (i_log) then
                write(logid, "('Beta            Gradient')")
                write(logid, "(F13.7, '    ', F13.7)") beta_current, beta_grad
            end if
               
            ! Check if gradient is within tolerance.
            if (abs(beta_grad) .le. tol) then
 
                ! Update beta_out, and note convergence, on all ranks.
                beta_out = beta_current
                converged = 'Y'
 
                ! Sum the chain energy-sums and store in rank 0
                my_ensum = sum(energy_chain)
                call mpi_reduce(my_ensum,energy_out,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
                ! Sum the chain variances and store in rank 0
                call mpi_reduce(variance_out, sum_var,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
                ! Sum the chain acceptance rates and store in rank 0
                call mpi_reduce(accept_out, sum_acc,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
                ! Calculate averages on rank 0
                if (my_rank==0) then
 
                    energy_out = energy_out * per_step
                    variance_out = sum_var * per_proc2
                    accept_out = sum_acc * per_proc
 
                    if (write_log) then
                        write(logid, "('Parameter search converged!')")
                        write(logid, "('Results:')")
                        write(logid, "('Beta = ', F13.7, ', Energy = ', F12.8, ', Variance = ', F12.8)") &
                            beta_current, energy_out, variance_out
                        write(logid, "('')")
                    end if
 
                endif 
 
                ! Exit subroutine, as convergence is reached.
                ! Note: Only rank 0 has the correct array outputs. All ranks have correct beta_out and 'converged'.
                exit
 
            end if
 
            ! If not, update beta by steepest descent, on all ranks.
            call h_2_update_beta(beta_current, beta_grad, beta_new)
            beta_current = beta_new
 
            ! Write restart files, if required
            if (write_restart .and. my_rank==0) then
                if (mod(k,restart_num)==0) then
                    ! Write res_nml.txt
                    call write_res_nml(i_bond,k+1,h2%auto_params,beta_current)
                endif
            endif 
 
        end do
 
        ! Deal with non-convergence
        if (k .ge. maxiter) then
            converged = 'N'
            beta_out = beta_current
            if (i_log) then 
                write(logid, "('Parameter search NOT converged! Consider changing&
                    & the starting value of beta, or the damping in the H_2_update_beta&
                    & function.')")
            end if
        end if
 
        ! Report results
        if (my_rank==0) then
            write(*, "(F12.8, '   ', A, '             ', F13.7, '      ', F12.8)") &
                bond_length, converged, beta_out, energy_out
        endif
 
        ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
        if (my_rank == 0) then
            if (write_chains) then
                write(str_b, "(F6.3)") bond_length
                write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
                write(*, "('Running VQMC at beta = ', F12.7, ' to produce single Markov chain...')") beta_out
                viz_steps = min(h2plus%steps,max_chain_len)
                call vqmc_h_2(h2%sigma, viz_steps, h2%burn_step, h2%thin_step, a, &
                        beta_out, r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
                        energy_chain, energy_out, variance_out, dphidbeta_chain)
                write(*, "(2A)") 'Writing Markov chain to ', chainfile
                call write_h2_equilibration(e1_pos_chain, energy_chain, chainfile, ierr, &
                    h2, accept_out, beta_out, bond_length)
            end if
        end if
        

Here is the call graph for this function:

Here is the caller graph for this function:

grid_h_2()

subroutine vqmc::grid_h_2	(	real(real64), intent(in)	a,
		real(real64), dimension(:), intent(in)	beta_array,
		real(real64), intent(in)	bond_length,
		real(real64), intent(out)	beta_out,
		real(real64), intent(out)	energy_out,
		real(real64), intent(out)	variance_out,
		real(real64), intent(out)	accept_out
	)

Performs a grid search to find the optimal value for the variational parameter \( \beta \) in the \( H_{2} \) problem.

Evaluates energy of the \( H_{2} \) trial wavefunction at each value of the parameter \( \beta \) in beta_array. Determines the best energy, and reports its corresponding \( \beta \) as the optimal value at this bond length. Also returns the variance of this energy as a measure of uncertainty, as well as the acceptance rate of the MCMC chain.

Parameters

[in]	a	Value of parameter a at current bond length, calculated by H_2_cusp
[in]	beta_array	Array of beta values to be tested
[in]	bond_length	H-H bond length
[out]	beta_out	Optimised value for beta at this bond length
[out]	energy_out	Energy of wavefunction at optimal beta
[out]	variance_out	Variance of energy at optimal beta
[out]	accept_out	Acceptance rate of MCMC chain at optimal beta

Definition at line 2153 of file vqmc.f90.

    
        implicit none
 
        real(real64), intent(in) :: a
        real(real64), dimension(:), intent(in) :: beta_array
        real(real64), intent(in) :: bond_length
        real(real64), intent(out) :: beta_out
        real(real64), intent(out) :: energy_out
        real(real64), intent(out) :: variance_out
        real(real64), intent(out) :: accept_out
 
        real(real64), dimension(3) :: R_a, R_b
        real(real64), dimension(:, :), allocatable :: e1_pos_chain, e2_pos_chain
        real(real64), dimension(:), allocatable :: energy_chain, DphiDbeta_chain
        real(real64), dimension(:), allocatable :: param_energies, param_accepts, param_variances
        integer :: viz_steps
        integer :: i
        integer :: ierr
        character(len=6) :: str_b
        character(len=40) :: chainfile
 
        ! MPI variables
        real(real64) :: my_ensum
 
        ! Initialise
        r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
        r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
        allocate(param_energies(size(beta_array))); param_energies(:)= 0.0_real64
        allocate(param_accepts(size(beta_array))); param_accepts(:) = 0.0_real64
        allocate(param_variances(size(beta_array))) ; param_variances(:) = 0.0_real64
 
        ! Sort out restarting
        if (run_restart) then
            ! Read in array values
            call read_energies(resfile_epar, param_energies, param_variances, param_accepts)
            ! Reset the tag
            run_restart = .false.
        endif 
 
        if (i_log) then
            write(logid, "('Starting parameter search for beta...')")
        end if
 
        ! Loop over beta values
        do i = grid_start, size(beta_array)
 
            ! Calculate energy at current value of beta
            call vqmc_h_2(h2%sigma, my_steps+h2%burn_step, h2%burn_step, h2%thin_step, a, &
                    beta_array(i), r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
                    energy_chain, energy_out, variance_out, dphidbeta_chain)
            
            ! Save energy - serial version
            !param_energies(i) = energy
            !param_variances(i) = variance
            !param_accepts(i) = accept_rate
 
            ! Sum the chain energy-sums and store in rank 0
            my_ensum = sum(energy_chain)
            call mpi_reduce(my_ensum,param_energies(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Sum the chain variances and store in rank 0
            call mpi_reduce(variance_out,param_variances(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Sum the chain acceptance rates and store in rank 0
            call mpi_reduce(accept_out,param_accepts(i),1,mpi_double_precision,mpi_sum,0,mpi_comm_world,mpierr)
 
            ! Calculate averages on rank 0
            if (my_rank==0) then
                param_energies(i) = param_energies(i) * per_step
                param_variances(i) = param_variances(i) * per_proc2
                param_accepts(i) = param_accepts(i) * per_proc
 
                ! Write restart files, if required
                if (write_restart) then
                    ! Write res_nml.txt
                    call write_res_nml(i_bond,i+1,h2%auto_params)
                    ! Write param energies, variances and acceptance rates
                    call write_qho_main(beta_array, param_energies, param_variances, resfile_epar, ierr, qho, param_accepts)
                endif 
            endif 
 
        end do
 
        ! Note: Only rank 0 has the correct _out values.
        if (my_rank==0) then
 
            ! Select lowest energy and corresponding beta value
            energy_out = minval(param_energies)
            beta_out = beta_array(minloc(param_energies, 1))
            variance_out = param_variances(minloc(param_energies, 1))
            accept_out = param_accepts(minloc(param_energies, 1))
 
            ! Report results
            write(*, "(F12.8, '   -             ', F13.7, '      ', F12.8)") &
                bond_length, beta_out, energy_out
 
            if (write_log) then
                write(logid, "('Results:')")
                write(logid, "('beta             Energy')")
                do i = 1, size(beta_array)
                    write(logid, "(F12.8, '    ', F12.8)") beta_array(i), param_energies(i)
                end do
                write(logid, "('')")
                write(logid, "('Min Energy = ', F12.8, ' at c = ', F12.8)") energy_out, beta_out
                write(logid, "('')")
            end if
 
            ! Run a new visualisation chain at the optimised parameter if write_chains is enabled.
            if (write_chains) then
                write(str_b, "(F6.3)") bond_length
                write(chainfile, "(3A)") './chains_out/bond_', trim(adjustl(str_b)), '.nc'
                write(*, "('Running VQMC at beta = ', F12.7, ' to produce single Markov chain...')") beta_out
                viz_steps = min(h2plus%steps,max_chain_len)
                call vqmc_h_2(h2%sigma, viz_steps, h2%burn_step, h2%thin_step, a, &
                        beta_out, r_a, r_b, e1_pos_chain, e2_pos_chain, accept_out, &
                        energy_chain, energy_out, variance_out, dphidbeta_chain)
                write(*, "(2A)") 'Writing Markov chain to ', chainfile
                call write_h2_equilibration(e1_pos_chain, energy_chain, chainfile, ierr, &
                    h2, accept_out, beta_out, bond_length)
            end if
 
        endif  
    

Here is the call graph for this function:

Here is the caller graph for this function:

equil_h_2()

subroutine vqmc::equil_h_2	(	real(real64), intent(in)	a,
		real(real64), dimension(:), intent(in)	beta_array,
		real(real64), intent(in)	bond_length,
		integer, intent(in)	bonds_pos
	)

Performs VQMC on an array of values for the variational parameter \( \beta \), outputting energy traces for each.

Generates VQMC position and energy traces of the \( H_{2} \) trial wavefunction at each value of the parameter \( \beta \) in beta_array. Outputs each of these traces to a separate NetCDF file, to be used in determining optimal burning/thinning paramters for an accurate full VQMC run.

Parameters

[in]	a	Value of parameter a at current bond length, calculated by H_2_cusp
[in]	beta_array	The array of beta values to obtain energy traces for
[in]	bond_length	The current H-H bond length
[in]	bonds_pos	The position in the array of bond lengths, to be used in designating output files

Definition at line 2299 of file vqmc.f90.

 
        implicit none
 
        real(real64), intent(in) :: a
        real(real64), dimension(:), intent(in) :: beta_array
        real(real64), intent(in) :: bond_length
        integer, intent(in) :: bonds_pos
 
        real(real64), dimension(3) :: R_a, R_b
        real(real64), dimension(:), allocatable :: variances
        real(real64), dimension(:), allocatable :: energies
        real(real64), dimension(:, :), allocatable :: e1_pos_chain, e2_pos_chain
        real(real64), dimension(:), allocatable :: energy_chain, DphiDbeta_chain
        real(real64) :: accept_rate
        real(real64) :: energy_mean
        real(real64) :: energy_variance
        integer :: i, ierr
        character(len=3) :: str_i, str_b
        character(len=40) :: equilfile
 
        ! Initialise
        r_a = (/-bond_length/2, 0.0_real64, 0.0_real64/)
        r_b = (/bond_length/2, 0.0_real64, 0.0_real64/)
        allocate(energies(size(beta_array))); energies(:) = 0.0_real64
        allocate(variances(size(beta_array)))
 
        ! Loop over values of beta
        do i = grid_start, size(beta_array)
 
            ! Calculate energy at current value of beta
            call vqmc_h_2(h2%sigma, h2%steps, h2%burn_step, h2%thin_step, a, &
                    beta_array(i), r_a, r_b, e1_pos_chain, e2_pos_chain, accept_rate, &
                    energy_chain, energy_mean, energy_variance, dphidbeta_chain)
            ! Save energy
            energies(i) = energy_mean
 
            write(*, "(F12.8, '   -             ', F13.7, '      ', F12.8)") &
                bond_length, beta_array(i), energies(i)
 
            ! Save the energy trace for each value of c
            write(str_i, "(I3)") i
            write(str_b, "(I3)") bonds_pos
            write(equilfile, "(5A)") './equil_out/b', trim(adjustl(str_b)), &
                '_beta', trim(adjustl(str_i)), '.nc'
            write(*, *) 'Writing equilibration outfile to ', equilfile
            call write_h2_equilibration(e1_pos_chain, energy_chain, equilfile, ierr, &
                h2, accept_rate, beta_array(i), bond_length)
 
            ! Write restart files, if required
            if (write_restart) then
                if (mod(i,restart_num)==0) then
                    call write_res_nml(i_bond,i+1,h2%auto_params)
                endif
            endif
 
        end do
 
        if (write_log) then
            write(logid, "('Results:')")
            write(logid, "('beta             Energy')")
            do i = 1, size(beta_array)
                write(logid, "(F12.8, '    ', F12.8)") beta_array(i), energies(i)
            end do
            write(logid, "('')")
        end if
 
        deallocate(energies)
        deallocate(variances)
    

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

ntab

integer, parameter vqmc::ntab = 32

Variable controlling the RNG sequence.

Definition at line 44 of file vqmc.f90.

    integer, parameter :: ntab = 32

iv

integer(int64), dimension(ntab), save vqmc::iv = 0

Variable controlling the RNG sequence.

Definition at line 45 of file vqmc.f90.

    integer(int64), dimension(ntab), save :: iv = 0

iy

integer(int64), save vqmc::iy = 0

Variable controlling the RNG sequence.

Definition at line 46 of file vqmc.f90.

    integer(int64), save :: iy = 0

logfrac

integer vqmc::logfrac = 10

1/logfrac is the fraction of samples to report in the log, if shared_data::write_log is TRUE.

Definition at line 49 of file vqmc.f90.

    integer :: logfrac = 10