// S&B.C
// Sites and bonds.
// Functions for initialization of the spin model
// © 2021 Peter J. Meyer
#include "iss.h"
unsigned int total_num_opened_bonds = 0;
void initialize_lattice_and_mark_boundary_sites(void);
static int de_occupy_sites(void);
// On entry to this function the sites array has been initialized
// showing all sites occupied by up-spins.
// This function (a) marks the boundary sites of the lattice,
// (b) optionally precomputes nn sites,
// (b) in a site-diluted lattice it marks some sites as unoccupied,
// (c) it opens bonds between sites (randomly if the lattice is bond-diluted).
// Upon return the sites array shows the occupied sites
// and the bonds array shows the open bonds.
// This assumes that the PRNG has been initialized.
/*----------------------------*/
void set_up_configuration(void)
{
int i;
double site_pc, bond_pc;
initialize_lattice_and_mark_boundary_sites();
num_spins = ( site_diluted ? de_occupy_sites() : num_sites );
num_opened_bonds = open_bonds(); // BONDS.C
total_num_opened_bonds += num_opened_bonds;
if ( goal_is_equilibration )
{
site_pc = ((double)num_spins*100)/num_sites;
if ( site_diluted )
{
printf("\n%s of %s (%d",ultoa_commas(num_spins,temp),ultoa_commas(num_sites,temp1),size);
for ( i=2; i<=dimensionality; i++ )
printf("x%d",size);
printf(") lattice sites occupied (%.2f%%).",site_pc);
}
bond_pc = ((double)num_opened_bonds*100)/num_pairs_nn;
if ( bond_diluted )
{
printf("\n%s bonds opened. This is %.2f%% of the %s pairs of nn spins.",
ultoa_commas(num_opened_bonds,temp),bond_pc,ultoa_commas(num_pairs_nn,temp2));
#if false
sprintf(temp,"%.3f",(double)num_opened_bonds/num_spins);
remove_trailing_zeros(temp);
printf("\nBonds per spin = %s",temp);
#endif
}
printf("\n");
if ( multiple_temperatures_specified )
{
sprintf(temp,"%.6f",temperature);
remove_trailing_zeros(temp);
printf("Temperature = %s ",temp);
// printf(" (%d/%d)",temperature_num+1,num_temperatures);
}
}
// At this point all spins are up-spins.
// For equilibration the assignment of initial spins
// is done by assign_initial_spins() in INITSPIN.C.
}
// On entry all sites are marked as occupied by up-spins.
// This function marks some sites as unoccupied.
// It is independent of the lattice type.
/*----------------------------*/
static int de_occupy_sites(void)
{
int i, j, k, l; // array indices
int num_occupied = num_sites;
int target_num_occupied = (int)(num_sites*site_concentration);
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( PRNG1() >= site_concentration )
{
sites2[i][j] &= UNOCCUPIED;
num_occupied--;
}
}
}
// If necessary, adjust site occupancy to get exact concentration specified.
while ( num_occupied > target_num_occupied )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
if ( SITE_IS_OCCUPIED_2(i,j) )
{
// Site is occupied, so de-occupy it.
sites2[i][j] &= UNOCCUPIED;
num_occupied--;
}
}
while ( num_occupied < target_num_occupied )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
if ( !SITE_IS_OCCUPIED_2(i,j) )
{
// Site is not occupied, so occupy it.
sites2[i][j] |= OCCUPIED;
num_occupied++;
}
}
break;
case 3:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
if ( PRNG1() >= site_concentration )
{
sites3[i][j][k] &= UNOCCUPIED;
num_occupied--;
}
}
}
}
while ( num_occupied > target_num_occupied )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
if ( SITE_IS_OCCUPIED_3(i,j,k) )
{
sites3[i][j][k] &= UNOCCUPIED;
num_occupied--;
}
}
while ( num_occupied < target_num_occupied )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
if ( !SITE_IS_OCCUPIED_3(i,j,k) )
{
sites3[i][j][k] |= OCCUPIED;
num_occupied++;
}
}
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 ( PRNG1() >= site_concentration )
{
sites4[i][j][k][l] &= UNOCCUPIED;
num_occupied--;
}
}
}
}
}
while ( num_occupied > target_num_occupied )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
l = (int)(size*PRNG1());
if ( SITE_IS_OCCUPIED_4(i,j,k,l) )
{
sites4[i][j][k][l] &= UNOCCUPIED;
num_occupied--;
}
}
while ( num_occupied < target_num_occupied )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
l = (int)(size*PRNG1());
if ( !SITE_IS_OCCUPIED_4(i,j,k,l) )
{
sites4[i][j][k][l] |= OCCUPIED;
num_occupied++;
}
}
break;
}
return ( num_occupied );
}
// open_bonds() is in BONDS.C
/*--------------------------------------------------*/
void initialize_lattice_and_mark_boundary_sites(void)
{
int i, j, k;
unsigned char site_val = ( OCCUPIED_BY_UP_SPIN | BOUNDARY ); // Bit pattern 11000000
switch ( dimensionality )
{
case 2:
// First show all sites as occupied by up spins
// with all lattice bonds closed.
memset(&sites2[0][0],OCCUPIED_BY_UP_SPIN,num_sites);
memset(&bonds2[0][0],0,num_sites);
// Now mark the boundary sites.
for ( i=0; i<size; i++ )
{
sites2[0][i] = sites2[size-1][i] =
sites2[i][0] = sites2[i][size-1] = site_val;
}
break;
case 3:
memset(&sites3[0][0][0],OCCUPIED_BY_UP_SPIN,num_sites);
memset(&bonds3[0][0][0],0,num_sites);
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
sites3[0][i][j] = sites3[size-1][i][j] =
sites3[i][0][j] = sites3[i][size-1][j] =
sites3[i][j][0] = sites3[i][j][size-1] = site_val;
}
}
break;
case 4:
memset(&sites4[0][0][0][0],OCCUPIED_BY_UP_SPIN,num_sites);
memset(&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++ )
{
sites4[0][i][j][k] = sites4[size-1][i][j][k] =
sites4[i][0][j][k] = sites4[i][size-1][j][k] =
sites4[i][j][0][k] = sites4[i][j][size-1][k] =
sites4[i][j][k][0] = sites4[i][j][k][size-1] = site_val;
}
}
}
}
}
/*-------------------------------*/
void mark_all_sites_unvisited(void)
{
int i, j, k, l;
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
for ( j=0; j<size; j++ )
sites2[i][j] &= UNVISITED;
break;
case 3:
for ( i=0; i<size; i++ )
for ( j=0; j<size; j++ )
for ( k=0; k<size; k++ )
sites3[i][j][k] &= UNVISITED;
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++ )
sites4[i][j][k][l] &= UNVISITED;
break;
}
}
/*----------------------------------------------*/
void mark_all_sites_unvisited_and_clear_path(void)
{
int i, j, k, l;
unsigned char site_val = ( UNVISITED & NOT_IN_SPANNING_CLUSTER ); // Bit pattern 11001111
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
for ( j=0; j<size; j++ )
sites2[i][j] &= site_val;
break;
case 3:
for ( i=0; i<size; i++ )
for ( j=0; j<size; j++ )
for ( k=0; k<size; k++ )
sites3[i][j][k] &= site_val;
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++ )
sites4[i][j][k][l] &= site_val;
break;
}
}
#if false
// No longer needed.
// This is used in the Wolff algorithm.
/*-----------------*/
void copy_sites(void)
{
switch ( dimensionality )
{
case 2:
memcpy(&v_sites2[0][0], &sites2[0][0], num_sites);
break;
case 3:
memcpy(&v_sites3[0][0][0], &sites3[0][0][0], num_sites);
break;
case 4:
memcpy(&v_sites4[0][0][0][0], &sites4[0][0][0][0], num_sites);
break;
}
}
#endif