// SW-WANG.C
// Equilibration using Swendsen-Wang technique
// © 2021 Peter J. Meyer
#include "iss.h"
#define STACK_POINTER_CHECK false
#if TIMING_DIAGNOSTICS
clock_t create_v_cluster_clock = 0;
clock_t flip_sw_clusters_clock = 0;
clock_t flip_sw_cluster_clock = 0;
clock_t another_sw_spin_clock = 0;
clock_t flip_sw_spin_clock = 0;
#endif
unsigned int num_virtual_bonds_opened = 0;
static short int s_i, s_j, s_k, s_l, s_r;
static int s_spin_value, s_new_spin_value;
char sw_map_filepath[128]; // pathname of map file
void create_virtual_bonds(void);
static void flip_swendsen_wang_clusters(void);
static void flip_swendsen_wang_cluster(void);
static void flip_swendsen_wang_spin(void);
static int another_swendsen_wang_spin(void);
// Swendsen-Wang method has two parts:
// 1. At the beginning of each Monte Carlo step go through the lattice sequentially
// and create virtual bonds between same-spin-value sites with a probability p,
// where p = 1 - exp(-2J/temperature).
// 2. Then for each virtual cluster assign a random spin value to all sites
// within that cluster.
/*----------------------------------*/
void perform_swendsen_wang_sweep(void)
{
int context;
create_virtual_bonds();
if ( map_file_requested )
{
context = 2;
strcpy(sw_map_filepath,map_filepath);
*(strstr(sw_map_filepath,".MAP")+3) = '0'+context;
display_map(sw_map_filepath,2);
}
flip_swendsen_wang_clusters();
if ( map_file_requested )
{
context = 3;
strcpy(sw_map_filepath,map_filepath);
*(strstr(sw_map_filepath,".MAP")+3) = '0'+context;
display_map(sw_map_filepath,3);
}
}
// J #defined in ISS_DEF.H as 1.0.
// According to Coddington
// sw_probability = 1 - exp(-J/temperature);
// For the critical temperature this gives sw_probability = 0.3564.
// According to Gould & Tobochnik
// sw_probability = 1 - exp(-2*J/temperature),
// which gives sw_probability = 0.5858.
// Comparison with results from internal energy and specific heat,
// as measured with Glauber dynamics, shows that Gould & Tobochnik are correct.
/*---------------------------*/
void create_virtual_bonds(void)
{
int i, j, k, l, r;
#if TIMING_DIAGNOSTICS
clock_t start_clock, end_clock;
start_clock = clock();
#endif
nn_i = nn_j = nn_k = nn_l = 0;
mark_all_sites_unvisited(); // S&B.C
switch ( dimensionality )
{
case 2:
// Erase v_bonds2[][]
memset(&v_bonds2[0][0],0,num_sites);
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( ! ( s_spin_value = spin_at_site(i,j) ) )
// spin_at_site() returns 0 if site not occupied.
continue;
MARK_SITE_AS_VISITED_2(i,j); // S&B.C
// Does this site have non-visited nearest neighbour spins?
if ( !bonds2[i][j] )
continue;
// because it has none
get_nn_sites(i,j); // Q&A.C
// nn_sites[][] now holds the indices
// of the nearest neighbour sites in all directions
// (including those that cross the lattice boundaries).
// (no. of directions = coord_num)
for ( r=0; r<coord_num; r++ )
{
if ( bonds2[i][j] & bit[r] ) // if bond exists in direction r
{
nn_i = nn_sites[r][0];
nn_j = nn_sites[r][1];
if ( SITE_IS_VISITED_2(nn_i,nn_j) ) // ISS_DEF.C
continue; // because this possible virtual bond
// has been considered already.
if ( s_spin_value != spin_at_occupied_site(nn_i,nn_j) )
continue; // because virtual bonds opened
// only between spins with same spin value.
if ( PRNG1() < v_bond_probability )
{
// Open the virtual bond.
v_bonds2[i][j] |= bit[r];
// Mark the virtual bond as open at the other site.
v_bonds2[nn_i][nn_j] |= bit[inverse_direction((i+j)%2,r)];
num_virtual_bonds_opened++;
}
}
}
}
}
break;
case 3:
// Erase v_bonds3[][][]
memset(&v_bonds3[0][0][0],0,num_sites);
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
if ( ! ( s_spin_value = spin_at_site(i,j,k) ) )
// spin_at_site() returns 0 if site not occupied.
continue;
MARK_SITE_AS_VISITED_3(i,j,k); // S&B.C
// Does this site have non-visited nearest neighbour spins?
if ( !bonds3[i][j][k] )
continue;
// because it has none
get_nn_sites(i,j,k); // Q&A.C
// nn_sites[][] now holds the indices
// of the nearest neighbour sites in all directions
// (including those that cross the lattice boundaries).
// (no. of directions = coord_num)
for ( r=0; r<coord_num; r++ )
{
if ( bonds3[i][j][k] & bit[r] ) // if bond exists in direction r
{
nn_i = nn_sites[r][0];
nn_j = nn_sites[r][1];
nn_k = nn_sites[r][2];
if ( SITE_IS_VISITED_3(nn_i,nn_j,nn_k) ) // ISS_DEF.C
continue; // because this possible virtual bond
// has been considered already.
if ( s_spin_value != spin_at_occupied_site(nn_i,nn_j,nn_k) )
continue; // because virtual bonds opened
// only between spins with same spin value.
if ( PRNG1() < v_bond_probability )
{
// Open the virtual bond.
v_bonds3[i][j][k] |= bit[r];
// Mark the virtual bond as open at the other site.
v_bonds3[nn_i][nn_j][nn_k]
|= bit[inverse_direction((i+j+k)%2,r)];
num_virtual_bonds_opened++;
}
}
}
}
}
}
break;
case 4:
// Erase v_bonds4[][][][]
memset(&v_bonds4[0][0][0][0],0,num_sites);
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
for ( l=0; l<size; l++ )
{
if ( ! ( s_spin_value = spin_at_site(i,j,k,l) ) )
// spin_at_site() returns 0 if site not occupied.
continue;
MARK_SITE_AS_VISITED_4(i,j,k,l); // S&B.C
// Does this site have non-visited nearest neighbour spins?
if ( !bonds4[i][j][k][l] )
continue;
// because it has none
get_nn_sites(i,j,k,l); // Q&A.C
// nn_sites[...] now holds the indices
// of the nearest neighbour sites in all directions
// (including those that cross the lattice boundaries).
// (no. of directions = coord_num)
for ( r=0; r<coord_num; r++ )
{
if ( bonds4[i][j][k][l] & bit[r] ) // if bond exists in direction r
{
nn_i = nn_sites[r][0];
nn_j = nn_sites[r][1];
nn_k = nn_sites[r][2];
nn_l = nn_sites[r][3];
if ( SITE_IS_VISITED_4(nn_i,nn_j,nn_k,nn_l) )
continue; // because this possible virtual bond
// has been considered already.
if ( s_spin_value != spin_at_occupied_site(nn_i,nn_j,nn_k,nn_l) )
continue; // because virtual bonds opened
// only between spins with same spin value.
if ( PRNG1() < v_bond_probability )
{
// Open the virtual bond.
v_bonds4[i][j][k][l] |= bit[r];
// Mark the virtual bond as open at the other site.
v_bonds4[nn_i][nn_j][nn_k][nn_l]
|= bit[inverse_direction((i+j+k+l)%2,r)];
num_virtual_bonds_opened++;
}
}
}
}
}
}
}
break;
}
#if TIMING_DIAGNOSTICS
end_clock = clock();
create_v_cluster_clock += end_clock - start_clock;
#endif
}
/*-----------------------------------------*/
static void flip_swendsen_wang_clusters(void)
{
int i, j, k, l;
#if TIMING_DIAGNOSTICS
clock_t start_clock, end_clock;
start_clock = clock();
#endif
mark_all_sites_unvisited();
clear_stack();
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( !SITE_IS_OCCUPIED_2(i,j) )
continue;
if ( SITE_IS_VISITED_2(i,j) )
continue;
// About to trace out a new cluster.
s_spin_value = spin_at_occupied_site(i,j);
// Pick a random spin value.
if ( model_is_ising )
s_new_spin_value = ( PRNG1() < 0.5 ? UP_SPIN : DOWN_SPIN );
else // model is q-state Potts
s_new_spin_value = (int)(PRNG1()*q) + 1;
s_i = i;
s_j = j;
ct_size_of_cluster = 0;
flip_swendsen_wang_cluster();
// Need to call this function even if spin value unchanged
// so as to mark sites in this cluster as visited.
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
num_spins_visited_this_sweep += ct_size_of_cluster;
num_clusters_traced_this_sweep++;
}
}
break;
case 3:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
if ( !SITE_IS_OCCUPIED_3(i,j,k) )
continue;
if ( SITE_IS_VISITED_3(i,j,k) )
continue;
// About to trace out a new cluster.
s_spin_value = spin_at_occupied_site(i,j,k);
// Pick a random spin value.
if ( model_is_ising )
s_new_spin_value = ( PRNG1() < 0.5 ? UP_SPIN : DOWN_SPIN );
else // model is q-state Potts
s_new_spin_value = (int)(PRNG1()*q) + 1;
s_i = i;
s_j = j;
s_k = k;
ct_size_of_cluster = 0;
flip_swendsen_wang_cluster();
// Need to call this function even if spin value unchanged
// so as to mark sites in this cluster as visited.
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
num_spins_visited_this_sweep += ct_size_of_cluster;
num_clusters_traced_this_sweep++;
}
}
}
break;
case 4:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
for ( l=0; l<size; l++ )
{
if ( !SITE_IS_OCCUPIED_4(i,j,k,l) )
continue;
if ( SITE_IS_VISITED_4(i,j,k,l) )
continue;
// About to trace out a new cluster.
s_spin_value = spin_at_occupied_site(i,j,k,l);
// Pick a random spin value.
if ( model_is_ising )
s_new_spin_value = ( PRNG1() < 0.5 ? UP_SPIN : DOWN_SPIN );
else // model is q-state Potts
s_new_spin_value = (int)(PRNG1()*q) + 1;
s_i = i;
s_j = j;
s_k = k;
s_l = l;
ct_size_of_cluster = 0;
flip_swendsen_wang_cluster();
// Need to call this function even if spin value unchanged
// so as to mark sites in this cluster as visited.
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
num_spins_visited_this_sweep += ct_size_of_cluster;
num_clusters_traced_this_sweep++;
}
}
}
}
break;
}
if ( num_spins_visited_this_sweep != (unsigned int)num_spins )
{
printf("\nError in SW-WANG.C: "
"num_spins_visited_this_sweep != num_spins\n");
exit(1);
}
#if TIMING_DIAGNOSTICS
end_clock = clock();
flip_sw_clusters_clock += end_clock - start_clock;
#endif
}
/*
Flip a Swendsen-Wang cluster (call with 1st spin at i,j):
Flip spin i,j (if new value different)
Loop
While there is a direction r such that there is a nn spin
in direction r which is part of this virtual cluster
and which is invisited:
Push i,j on stack
Put i,j = location of spin in direction r
Flip spin i,j
End While
If stack is empty
finish
End If
Pop location from stack to i,j
End Loop
Flip spin:
If new spin value different, flip the spin
Mark spin as visited
Increase cluster count
*/
/*----------------------------------------*/
static void flip_swendsen_wang_cluster(void)
{
#if TIMING_DIAGNOSTICS
clock_t start_clock, end_clock;
start_clock = clock();
#endif
flip_swendsen_wang_spin();
switch ( dimensionality )
{
case 2:
loop
{
while ( another_swendsen_wang_spin() )
{
// direction of possible new spin is s_r
// and location of possible new spin is (nn_i,nn_j).
if ( push_onto_stack(s_i,s_j) )
{
s_i = nn_i;
s_j = nn_j;
flip_swendsen_wang_spin();
}
}
if ( stack_empty() )
{
#if TIMING_DIAGNOSTICS
end_clock = clock();
flip_sw_cluster_clock += end_clock - start_clock;
#endif
return;
}
pop_from_stack(&s_i,&s_j);
}
break;
case 3:
loop
{
while ( another_swendsen_wang_spin() )
{
// direction of possible new spin is s_r
// and location of possible new spin is (nn_i,nn_j,nn_k).
if ( push_onto_stack(s_i,s_j,s_k) )
{
s_i = nn_i;
s_j = nn_j;
s_k = nn_k;
flip_swendsen_wang_spin();
}
}
if ( stack_empty() )
{
#if TIMING_DIAGNOSTICS
end_clock = clock();
flip_sw_cluster_clock += end_clock - start_clock;
#endif
return;
}
pop_from_stack(&s_i,&s_j,&s_k);
}
break;
case 4:
loop
{
while ( another_swendsen_wang_spin() )
{
// direction of possible new spin is s_r
// and location of possible new spin is (nn_i,nn_j,nn_k,nn_l).
if ( push_onto_stack(s_i,s_j,s_k,s_l) )
{
s_i = nn_i;
s_j = nn_j;
s_k = nn_k;
s_l = nn_l;
flip_swendsen_wang_spin();
}
}
if ( stack_empty() )
{
#if TIMING_DIAGNOSTICS
end_clock = clock();
flip_sw_cluster_clock += end_clock - start_clock;
#endif
return;
}
pop_from_stack(&s_i,&s_j,&s_k,&s_l);
}
break;
}
#if TIMING_DIAGNOSTICS
end_clock = clock();
flip_sw_cluster_clock += end_clock - start_clock;
#endif
}
/*-------------------------------------*/
static void flip_swendsen_wang_spin(void)
{
#if TIMING_DIAGNOSTICS
clock_t start_clock, end_clock;
start_clock = clock();
#endif
ct_size_of_cluster++;
switch ( dimensionality )
{
case 2:
if ( s_spin_value != s_new_spin_value )
{
// Flip the spin.
if ( model_is_ising )
do_ising_spin_flip(s_spin_value,s_i,s_j);
else // model is q-state Potts
do_q_state_potts_spin_flip(s_spin_value,s_new_spin_value,s_i,s_j);
}
MARK_SITE_AS_VISITED_2(s_i,s_j);
break;
case 3:
if ( s_spin_value != s_new_spin_value )
{
// Flip the spin.
if ( model_is_ising )
do_ising_spin_flip(s_spin_value,s_i,s_j,s_k);
else // model is q-state Potts
do_q_state_potts_spin_flip(s_spin_value,s_new_spin_value,s_i,s_j,s_k);
}
MARK_SITE_AS_VISITED_3(s_i,s_j,s_k);
break;
case 4:
if ( s_spin_value != s_new_spin_value )
{
// Flip the spin.
if ( model_is_ising )
do_ising_spin_flip(s_spin_value,s_i,s_j,s_k,s_l);
else // model is q-state Potts
do_q_state_potts_spin_flip(s_spin_value,s_new_spin_value,s_i,s_j,s_k,s_l);
}
MARK_SITE_AS_VISITED_4(s_i,s_j,s_k,s_l);
break;
}
#if TIMING_DIAGNOSTICS
end_clock = clock();
flip_sw_spin_clock += end_clock - start_clock;
#endif
}
// This returns true if there is a direction r such that
// there is a nn spin in direction r which is part of this virtual cluster
// and which is unvisited.
// In this case the direction of the possible new spin on return is in s_r
// and the location of the possible new spin is (nn_i,nn_j,...).
/*---------------------------------------*/
static int another_swendsen_wang_spin(void)
{
#if TIMING_DIAGNOSTICS
clock_t start_clock, end_clock;
start_clock = clock();
#endif
switch ( dimensionality )
{
case 2:
if ( v_bonds2[s_i][s_j] )
{
// There is at least one virtual bond at this site.
get_nn_sites(s_i,s_j);
for ( s_r=0; s_r<coord_num; s_r++ )
{
if ( v_bonds2[s_i][s_j] & bit[s_r] )
// if a virtual bond exists in direction s_r
{
nn_i = nn_sites[s_r][0];
nn_j = nn_sites[s_r][1];
if ( !SITE_IS_VISITED_2(nn_i,nn_j) )
return ( true );
}
}
}
break;
case 3:
if ( v_bonds3[s_i][s_j][s_k] )
{
// There is at least one virtual bond at this site.
get_nn_sites(s_i,s_j,s_k);
for ( s_r=0; s_r<coord_num; s_r++ )
{
if ( v_bonds3[s_i][s_j][s_k] & bit[s_r] )
// if a virtual bond exists in direction s_r
{
nn_i = nn_sites[s_r][0];
nn_j = nn_sites[s_r][1];
nn_k = nn_sites[s_r][2];
if ( !SITE_IS_VISITED_3(nn_i,nn_j,nn_k) )
return ( true );
}
}
}
break;
case 4:
if ( v_bonds4[s_i][s_j][s_k][s_l] )
{
// There is at least one virtual bond at this site.
get_nn_sites(s_i,s_j,s_k,s_l);
for ( s_r=0; s_r<coord_num; s_r++ )
{
if ( v_bonds4[s_i][s_j][s_k][s_l] & bit[s_r] )
// if a virtual bond exists in direction s_r
{
nn_i = nn_sites[s_r][0];
nn_j = nn_sites[s_r][1];
nn_k = nn_sites[s_r][2];
nn_l = nn_sites[s_r][3];
if ( !SITE_IS_VISITED_4(nn_i,nn_j,nn_k,nn_l) )
return ( true );
}
}
}
break;
}
#if TIMING_DIAGNOSTICS
end_clock = clock();
another_sw_spin_clock += end_clock - start_clock;
#endif
return ( false );
}