Fortran Developer Documentation Main Page

This is the homepage for the Developer Documentation of the Fortran backend of Dice. On this site, you will find definitions for all of the functions and subroutines that we use to run VQMC simulations on our three currently implemented problems (QHO, H2plus and H2), as well as instructions for integrating support for a new problem into the code (see below).

Links

Link	Destination
GitHub	Here
Dev Documentation Main	Here
Fortran Dev Documentation	You are here
Python Dev Documentation	Here

---

Adding a new problem to Dice Quantum Monte Carlo (QMC)

This section presents what the user would need to modify in the Fortran files (.f90) to add a new problem to Dice QMC.

The first step is to add the derived type of the new problem to shared_data.f90, under the section "Derived Types". Below is a sample format of how the derived type should look like.

---

! Sample format of new derived type for new problem
type :: new_problem_type
        integer(int32) :: steps=100000_int32         ! Total number of steps for each MMC chain
        integer(int32) :: burn_step=1000_int32       ! Number of steps to be burnt from the start of the MMC chain
        integer(int32) :: thin_step=10_int32         ! Interval of steps by which the MMC chain is thinned
        real(real64) :: rcut=3.0_real64              ! Maximum radius of search around the atoms             
        real(real64) :: sigma=1.8_real64             ! Standard deviation in MMC transition probability
        logical :: auto_params = .true.              ! True runs the automatic gradient search and false runs the grid search
        real(real64) :: parameter_for_new_problem = 0.4_real64      ! One can add more than one parameter in the same format
        real(real64), dimension(3) :: parameter_grid = (/0,1,10/)   ! Grid of parameter values, format: [min, max, grid points]
        real(real64), dimension(3) :: bond = (/1,3,10/)             ! Grid of bondlength values, format: [min, max, grid points]
    end type H2_type

If the user wants to add any other numerical constants then they can be added to the constants module in shared_data.f90.

The next step is to modify quantum_solvers.f90. Regardless of the new problem the user will need to implement a function that calculates the transition probability of the current step and the proposed step. Below is a sample format of how the transition probability function should look like.

---

! Sample format of transition probability function
function trans_prob(parameters_of_new_problem, &
current_step, next_step, other_constants) result(Prob)
 
        implicit none
 
        real(real64), intent(in) :: parameters_of new_problem 
        real(real64), dimension(3), intent(in) :: current_step, next_step
        real(real64), dimension(:), intent(in) :: other_constants 
 
        real(real64) :: Prob
 
        prob = ...
 
    end function Trans_Prob

The other function that is also necessary to be implemented regardless of the problem is one that calculates the local energy of the problem at a given position. Below is a sample format of how the local energy function should look like.

---

! Sample format of local energy function 
function local_energy(parameters_of_new_problem, &
position) result(E_loc)
 
        implicit none
 
        real(real64), intent(in), dimension(:) :: parameters
        real(real64), intent(in), dimension(:) :: position 
 
        real(real64) :: E_loc
 
        e_loc = ...
 
    end function Local_Energy

Most problems will require the user to implement a wave function corresponding to the problem, in order to calculate the transitional probabilities and local energies. The user is free to implement this however they want, as long as the local energy and transitional probability functions have the formats presented above. Another remark for the user is that when calculating the energy, unless the derivatives can be found analytically the user will need to use some numerical method to obtain the derivatives, which they will have to implement themselves. All of the modifications to quantum_solvers.f90 should be done after the section H_2 Solvers.

The next step is to modify the vqmc.f90 module. The user will need to implement a new VQMC solver for the new problem in this module. The user should use one of the "VQMC_QHO", "VQMC_H_2_plus" or "VQMC_H_2" sub-routines as a template for their VQMC solver. This is because they would only need to modify the transitional probability and local energy functions in the Metropolis algorithm. Below is a sample template of how the Metropolis algorithm should look like.

---

! Metropolis Algorithm without burning and thinning for new problem
do i = 2, n            
    ...
 
    ! Modify here the function for transitional probability
    log_prob_n = log(trans_prob(parameters, x(i-1), x_n(1)))
            
    ...
            
    ! Compare the drawn sample from the uniform distribution with the acceptance ratio 
    if (a_prob_trans >= unif_num) then 
        ...
        ! Modify here the function for local energy  
        energy_chain(i) = local_energy(parameters, x_n(1))
    else 
        ... 
        
        ! Modify here the function for local energy
        energy_chain(i) = local_energy(parameters, x(i-1))
    end if 
end do

If the user uses one of the "VQMC_QHO", "VQMC_H_2_plus" or "VQMC_H_2" sub-routines as a template then the burning and thinning features can be reused for the new problem. It is important to stick to one of the VQMC sub-routines so that the outputs are the same.

The next step is to add the automatic gradient descent-based search (GDBS) to vqmc.f90. The user should use "Auto_H_2", "Auto_H_2_plus" sub-routines as a template and should only modify the VQMC solver being called (see below). This modification should be place after the new VQMC solver.

---

subroutine auto_h_2(...)
        ...
        ! Main iteration loop
        do k = cauto_start, maxiter
 
            ! Calculate gradient at current beta
            ! This should be modified to the VQMC
            ! solver that user created in the 
            ! last step.
            call vqmc_h_2(...)
        ...
        end do 
        
end subroutine Auto_H_2

If the user does not want to do a GDBS then they can ignore the previous step. Instead they can add a grid search to vqmc.f90. The user should use "Grid_H_2", "Grid_H_2_plus" sub-routines as a template and should only modify the VQMC solver being called (see below). This sub-routine should be placed after the new VQMC solver.

---

subroutine grid_h_2_plus(...)
        ...
    ! Loop over c values
    do i = cgrid_start, size(c_array)
 
        ! Calculate energy at current value of c
        ! This should be modified to the VQMC
        ! solver that user created in the 
        ! last step.
        call vqmc_h_2_plus(...)
        ...
    end do
    ...
end subroutine Grid_H_2_plus

The user can of course add both of the search methods mentioned above, but it is important that they at least add one. The final item the user has to add to vqmc.f90 is the equilibration run subroutine. This can be based off "Equil_QHO" or "Equil_H_2" and the user should only modify the VQMC solver being used like in the previous step.

The final step is to modify driver.f90. The user should add a new case with the corresponding problem to the "VQMC Loops" section. The user should use the following "Quantum Harmonic Oscillator" subsection, as a template to write their code. The code should only be modified so that the GDBS or grid search, and equilibration subroutines for the new problem are called (see below).

---

case('New_problem')
    ...
    ! Loop over bond lengths
    do b = 1, int(h2plus%bond(3))
        ...
        ! Main VQMC callers
        if (.not. run_equil) then
            if (h2plus%auto_params) then
            ! Search for best param at this bond length
            ! This should be modified to the 
            ! GDBS sub-routine for this problem
            call auto_h_2_plus(...)
            ...
            else
 
            ! Find best param in grid at this bond length
            ! This should be modified to the 
            ! Grid search routine for this problem 
            call grid_h_2_plus(...)
 
            end if
        ...
        end if 
        ...
        if 
            ...
        else
            ! Run equilibrations serially on rank 0
            if (my_rank==0) then
                ! This shoud be modified to 
                ! the equilibration sub_routine 
                ! for this new problem
                call equil_h_2_plus(...)
            ...
        end if
    ...
    end do
...        

After all the modifications the user should be able to run the code the same way as it is explained in the user manual.

This documentation was generated from commit ece8004