// WOLFF.C
// Wolff cluster algorithm
// © 2021 Peter J. Meyer
#include "iss.h"
#define STACK_POINTER_CHECK true
#define DEBUG false
unsigned int num_spins_visited_this_run; // set to 0 in EQUIL.C at start of each run
char wolff_map_filepath[128]; // pathname of map file
short int w_i, w_j, w_k, w_l, w_r, r0;
// not static, since accessed in DISP_LAT.C
static int w_spin_value, w_new_spin_value;
static int perform_wolff_update(void);
static void grow_and_flip_wolff_cluster(void);
static void flip_wolff_spin(void);
static int possible_new_wolff_spin(void);
// This implementation is slow for T considerably larger than
// the critical temperature (e.g. T = 4.0), but is about the
// same speed as the Swendsen-Wang algorithm at the critical temperature.
/*--------------------------*/
void perform_wolff_sweep(void)
{
unsigned int target_num, target_num_spins_visited_this_run;
if ( model_is_ising )
target_num_spins_visited_this_run = ( mcs_num+1 )*((unsigned int)num_spins/2);
else // model is q-state Potts
target_num_spins_visited_this_run = ( mcs_num+1 )*(((unsigned int)num_spins*(q-1))/q);
if ( target_num_spins_visited_this_run < num_spins_visited_this_run )
target_num = 0;
else
target_num = target_num_spins_visited_this_run - num_spins_visited_this_run;
do {
num_spins_visited_this_sweep += perform_wolff_update();
num_clusters_traced_this_sweep++;
if ( map_file_requested )
{
strcpy(wolff_map_filepath,map_filepath);
display_map(wolff_map_filepath,4);
}
} while ( num_spins_visited_this_sweep < target_num );
num_spins_visited_this_run += num_spins_visited_this_sweep;
}
/*---------------------------------*/
static int perform_wolff_update(void)
{
mark_all_sites_unvisited();
switch ( dimensionality )
{
case 2:
// Select spin at random.
do {
w_i = (int)(PRNG1()*size);
w_j = (int)(PRNG1()*size);
} while ( !SITE_IS_OCCUPIED_2(w_i,w_j) );
w_spin_value = spin_at_occupied_site(w_i,w_j);
break;
case 3:
do {
w_i = (int)(PRNG1()*size);
w_j = (int)(PRNG1()*size);
w_k = (int)(PRNG1()*size);
} while ( !SITE_IS_OCCUPIED_3(w_i,w_j,w_k) );
w_spin_value = spin_at_occupied_site(w_i,w_j,w_k);
break;
case 4:
do {
w_i = (int)(PRNG1()*size);
w_j = (int)(PRNG1()*size);
w_k = (int)(PRNG1()*size);
w_l = (int)(PRNG1()*size);
} while ( !SITE_IS_OCCUPIED_4(w_i,w_j,w_k,w_l) );
w_spin_value = spin_at_occupied_site(w_i,w_j,w_k,w_l);
break;
}
if ( model_is_q_state_potts )
{
// Pick a new q-state Potts spin value
// randomly selected from values other than the old.
do {
w_new_spin_value = (int)(PRNG1()*q) + 1;
} while ( w_new_spin_value == w_spin_value );
}
clear_stack();
ct_size_of_cluster = 0;
grow_and_flip_wolff_cluster();
if ( ct_size_of_cluster > max_cluster_size )
max_cluster_size = ct_size_of_cluster;
#if STACK_POINTER_CHECK
check_stack_is_empty(__FILE__);
#endif
return ( ct_size_of_cluster );
}
/*
Grow and flip a Wolff cluster (call with 1st spin at i,j):
Flip spin.
Put r0 = -1
Loop
While there is a direction r > r0 such that there is a nn spin in direction r
which has the same value as the seed spin and is unvisited:
If PRNG < probability
Push i,j,r on stack
Put i,j = location of new spin
Put r0 = -1
Flip spin
End If
End While
If stack is empty
finish
End If
Pop i,j,r0 from stack
End Loop
Flip spin:
Flip the spin
Mark spin as visited
Increase cluster count
*/
/*-----------------------------------------*/
static void grow_and_flip_wolff_cluster(void)
{
flip_wolff_spin();
r0 = -1;
switch ( dimensionality )
{
case 2:
if ( map_file_requested )
// Erase v_bonds2
memset(&v_bonds2[0][0],0,num_sites);
loop
{
while ( possible_new_wolff_spin() )
{
// direction of possible new spin is w_r
// and location of possible new spin is (nn_i,nn_j).
if ( PRNG1() < v_bond_probability )
{
if ( push_onto_stack_with_r(w_r,w_i,w_j) )
{
w_i = nn_i;
w_j = nn_j;
r0 = -1;
flip_wolff_spin();
#if DEBUG
if ( map_file_requested )
{
strcpy(wolff_map_filepath,map_filepath);
display_map(wolff_map_filepath,4);
}
#endif
}
}
else
r0 = w_r;
}
if ( stack_empty() )
return;
pop_from_stack_with_r(&r0,&w_i,&w_j);
}
break;
case 3:
if ( map_file_requested )
// Erase v_bonds3
memset(&v_bonds3[0][0][0],0,num_sites);
loop
{
while ( possible_new_wolff_spin() )
{
// direction of possible new spin is w_r
// and location of possible new spin is (nn_i,nn_j,nn_k).
if ( PRNG1() < v_bond_probability )
{
if ( push_onto_stack_with_r(w_r,w_i,w_j,w_k) )
{
w_i = nn_i;
w_j = nn_j;
w_k = nn_k;
r0 = -1;
flip_wolff_spin();
#if DEBUG
if ( map_file_requested )
{
strcpy(wolff_map_filepath,map_filepath);
display_map(wolff_map_filepath,4);
}
#endif
}
}
else
r0 = w_r;
}
if ( stack_empty() )
return;
pop_from_stack_with_r(&r0,&w_i,&w_j,&w_k);
}
break;
case 4:
if ( map_file_requested )
// Erase v_bonds3
memset(&v_bonds4[0][0][0][0],0,num_sites);
loop
{
while ( possible_new_wolff_spin() )
{
// direction of possible new spin is w_r
// and location of possible new spin is (nn_i,nn_j,nn_k,nn_l).
if ( PRNG1() < v_bond_probability )
{
if ( push_onto_stack_with_r(w_r,w_i,w_j,w_k,w_l) )
{
w_i = nn_i;
w_j = nn_j;
w_k = nn_k;
w_l = nn_l;
r0 = -1;
flip_wolff_spin();
#if DEBUG
if ( map_file_requested )
{
strcpy(wolff_map_filepath,map_filepath);
display_map(wolff_map_filepath,4);
}
#endif
}
}
else
r0 = w_r;
}
if ( stack_empty() )
return;
pop_from_stack_with_r(&r0,&w_i,&w_j,&w_k,&w_l);
}
break;
}
}
/*-----------------------------*/
static void flip_wolff_spin(void)
{
ct_size_of_cluster++;
switch ( dimensionality )
{
case 2:
// Flip the spin.
if ( model_is_ising )
do_ising_spin_flip(w_spin_value,w_i,w_j);
else // model is q-state Potts
do_q_state_potts_spin_flip(w_spin_value,w_new_spin_value,w_i,w_j);
MARK_SITE_AS_VISITED_2(w_i,w_j);
break;
case 3:
if ( model_is_ising )
do_ising_spin_flip(w_spin_value,w_i,w_j,w_k);
else // model is q-state Potts
do_q_state_potts_spin_flip(w_spin_value,w_new_spin_value,w_i,w_j,w_k);
MARK_SITE_AS_VISITED_3(w_i,w_j,w_k);
break;
case 4:
if ( model_is_ising )
do_ising_spin_flip(w_spin_value,w_i,w_j,w_k,w_l);
else // model is q-state Potts
do_q_state_potts_spin_flip(w_spin_value,w_new_spin_value,w_i,w_j,w_k,w_l);
MARK_SITE_AS_VISITED_4(w_i,w_j,w_k,w_l);
break;
}
}
// This returns true if there is a direction r > r0 such that
// there is a nn spin in direction r which has the same spin value
// as the seed spin and which is unvisited.
// In this case the direction of the possible new spin on return is in w_r
// and the location of the possible new spin is (nn_i,nn_j,...).
/*------------------------------------*/
static int possible_new_wolff_spin(void)
{
if ( r0 == coord_num - 1 )
return ( false );
switch ( dimensionality )
{
case 2:
if ( !bonds2[w_i][w_j] )
return ( false );
// There is at least one open bond at this site.
get_nn_sites(w_i,w_j);
for ( w_r=r0+1; w_r<coord_num; w_r++ )
{
if ( bonds2[w_i][w_j] & bit[w_r] )
// if an open bond exists in direction w_r
{
nn_i = nn_sites[w_r][0];
nn_j = nn_sites[w_r][1];
if ( !SITE_IS_VISITED_2(nn_i,nn_j) )
// if the nn spin has not been visited (=flipped)
{
// This spin has not been visited, so it has not been flipped
// so it has the same value as when we began to grow the cluster.
if ( w_spin_value == spin_at_occupied_site(nn_i,nn_j) )
// if this spin has the same spin value as the seed spin
{
if ( map_file_requested )
{
// Open the virtual bond to show that
// this bond has been considered.
v_bonds2[w_i][w_j] |= bit[w_r];
// Mark the virtual bond as open at the other site.
v_bonds2[nn_i][nn_j]
|= bit[inverse_direction((w_i+w_j)%2,w_r)];
}
return ( true );
}
}
}
}
break;
case 3:
if ( !bonds3[w_i][w_j][w_k] )
return ( false );
// There is at least one open bond at this site.
get_nn_sites(w_i,w_j,w_k);
for ( w_r=r0+1; w_r<coord_num; w_r++ )
{
if ( bonds3[w_i][w_j][w_k] & bit[w_r] )
// if an open bond exists in direction w_r
{
nn_i = nn_sites[w_r][0];
nn_j = nn_sites[w_r][1];
nn_k = nn_sites[w_r][2];
if ( !SITE_IS_VISITED_3(nn_i,nn_j,nn_k) )
// if the nn spin has not been visited (=flipped)
{
// This spin has not been visited, so it has not been flipped
// so it has the same value as when we began to grow the cluster.
if ( w_spin_value == spin_at_occupied_site(nn_i,nn_j,nn_k) )
// if this spin has the same spin value as the seed spin
{
if ( map_file_requested )
{
// Open the virtual bond to show that
// this bond has been considered.
v_bonds3[w_i][w_j][w_k] |= bit[w_r];
// Mark the virtual bond as open at the other site.
v_bonds3[nn_i][nn_j][nn_k]
|= bit[inverse_direction((w_i+w_j+w_k)%2,w_r)];
}
return ( true );
}
}
}
}
break;
case 4:
if ( !bonds4[w_i][w_j][w_k][w_l] )
return ( false );
// There is at least one open bond at this site.
get_nn_sites(w_i,w_j,w_k,w_l);
for ( w_r=r0+1; w_r<coord_num; w_r++ )
{
if ( bonds4[w_i][w_j][w_k][w_l] & bit[w_r] )
// if an open bond exists in direction w_r
{
nn_i = nn_sites[w_r][0];
nn_j = nn_sites[w_r][1];
nn_k = nn_sites[w_r][2];
nn_l = nn_sites[w_r][3];
if ( !SITE_IS_VISITED_4(nn_i,nn_j,nn_k,nn_l) )
// if the nn spin has not been visited (=flipped)
{
// This spin has not been visited, so it has not been flipped
// so it has the same value as when we began to grow the cluster.
if ( w_spin_value == spin_at_occupied_site(nn_i,nn_j,nn_k,nn_l) )
// if this spin has the same spin value as the seed spin
{
if ( map_file_requested )
{
// Open the virtual bond to show that
// this bond has been considered.
v_bonds4[w_i][w_j][w_k][w_l] |= bit[w_r];
// Mark the virtual bond as open at the other site.
v_bonds4[nn_i][nn_j][nn_k][nn_l]
|= bit[inverse_direction((w_i+w_j+w_k+w_l)%2,w_r)];
}
return ( true );
}
}
}
}
break;
}
return ( false );
}