// BONDS.C
// Functions for setting up bonds among occupied sites
// © 2021 Peter J. Meyer
#include "iss.h"
// nn_i, nn_j, nn_k, nn_l are global ints.
// This function depends on the lattice type.
/*----------------*/
int open_bonds(void)
{
int i, j, k, l, r, nn_r;
int num_opened = 0, target_num_opened;
unsigned char bonds, nn_bonds;
num_pairs_nn = 0;
if ( bond_diluted ) // bond_diluted set in SETPARAM.C.
mark_all_sites_unvisited(); // below
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( site_diluted )
if ( !SITE_IS_OCCUPIED_2(i,j) )
continue;
get_nn_sites(i,j);
// nn_sites[][] now holds the indices
// of the nearest neighbour sites in all directions
// (no. of directions = coord_num).
// Now check which nn sites are occupied
// and set bonds2[i][j] appropriately.
bonds = 0;
// Must use following "for" so that the right-shift is done 8 times.
for ( r=0; r<MAX_COORD_NUM; r++ )
{
bonds >>= 1;
if ( r >= coord_num )
continue;
if ( ( nn_i = nn_sites[r][0] ) == NON_EXISTENT )
continue;
// only with free boundary conditions
nn_j = nn_sites[r][1];
// Check if nn site is occupied.
if ( !SITE_IS_OCCUPIED_2(nn_i,nn_j) )
continue;
// nn site exists and is occupied.
num_pairs_nn++;
if ( !bond_diluted )
{
bonds |= OCCUPIED;
num_opened++; // Done twice for each opened bond.
}
else // bond diluted
{
// Check if bonds at this nn site have already been opened
// (or at least, have been considered for opening).
if ( SITE_IS_VISITED_2(nn_i,nn_j) )
{
// Yes, so check if bond has already been opened.
nn_bonds = bonds2[nn_i][nn_j];
// Get direction from the nn site back to this one.
nn_r = inverse_direction((i+j)%2,r); // DIR_TABLE.C
if ( nn_bonds & bit[nn_r] )
{
// Bond has been opened, so mark it here.
bonds |= OCCUPIED; // set bit 7
// Don't do num_opened++
}
}
else // Site not yet visited.
// Decide whether to open this bond.
{
if ( PRNG1() < bond_concentration )
{
bonds |= OCCUPIED;
num_opened++;
}
}
}
}
bonds2[i][j] = bonds;
if ( bond_diluted )
MARK_SITE_AS_VISITED_2(i,j);
}
}
num_pairs_nn /= 2; // since each pair has been counted twice.
if ( !bond_diluted )
num_opened /= 2; // since each bond has been counted twice.
// If necessary, adjust bonds to get exact concentration specified.
target_num_opened = (int)(num_pairs_nn*bond_concentration);
while ( num_opened > target_num_opened )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
if ( bonds2[i][j] )
{
// This spin has an open bond.
for ( r=0; r<coord_num; r++ )
{
if ( bonds2[i][j] & bit[r] )
{
// Bond is open, so close it.
bonds2[i][j] &= unbit[r];
// Mark the bond closed at the nn site.
get_nn_site_in_dir(r,i,j);
bonds2[nn_site[0]][nn_site[1]] &= unbit[inverse_direction((i+j)%2,r)];
num_opened--;
break;
}
}
}
}
while ( num_opened < target_num_opened )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
if ( SITE_IS_OCCUPIED_2(i,j) )
{
// Look for a closed bond.
for ( r=0; r<coord_num; r++ )
{
if ( ! ( bonds2[i][j] & bit[r] ) )
{
// Bond is closed, so open it (assuming there is
// an occupied nn site in direction r).
if ( get_nn_site_in_dir(r,i,j) )
{
// Nearest neighbour site exists.
if ( SITE_IS_OCCUPIED_2(nn_site[0],nn_site[1]) )
{
bonds2[i][j] |= bit[r];
// Mark the bond open at the nn site.
bonds2[nn_site[0]][nn_site[1]]
|= bit[inverse_direction((i+j)%2,r)];
num_opened++;
break;
}
}
}
}
}
}
break;
case 3:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
if ( site_diluted )
if ( !SITE_IS_OCCUPIED_3(i,j,k) )
continue;
get_nn_sites(i,j,k);
// See comments for 2d case.
bonds = 0;
for ( r=0; r<MAX_COORD_NUM; r++ )
{
bonds >>= 1;
if ( r >= coord_num )
continue;
if ( ( nn_i = nn_sites[r][0] ) == NON_EXISTENT )
continue;
nn_j = nn_sites[r][1];
nn_k = nn_sites[r][2];
if ( !SITE_IS_OCCUPIED_3(nn_i,nn_j,nn_k) )
continue;
num_pairs_nn++;
if ( !bond_diluted )
{
bonds |= OCCUPIED;
num_opened++;
}
else
{
if ( SITE_IS_VISITED_3(nn_i,nn_j,nn_k) )
{
nn_bonds = bonds3[nn_i][nn_j][nn_k];
nn_r = inverse_direction((i+j+k)%2,r);
if ( nn_bonds & bit[nn_r] )
bonds |= OCCUPIED;
}
else
{
if ( PRNG1() < bond_concentration )
{
bonds |= OCCUPIED;
num_opened++;
}
}
}
}
bonds3[i][j][k] = bonds;
if ( bond_diluted )
MARK_SITE_AS_VISITED_3(i,j,k);
}
}
}
num_pairs_nn /= 2; // since each pair has been counted twice.
if ( !bond_diluted )
num_opened /= 2; // since each bond has been counted twice.
// If necessary, adjust bonds to get exact concentration specified.
target_num_opened = (int)(num_pairs_nn*bond_concentration);
while ( num_opened > target_num_opened )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
if ( bonds3[i][j][k] )
{
// This spin has an open bond.
for ( r=0; r<coord_num; r++ )
{
if ( bonds3[i][j][k] & bit[r] )
{
// Bond is open, so close it.
bonds3[i][j][k] &= unbit[r];
// Mark the bond closed at the nn site.
get_nn_site_in_dir(r,i,j,k);
bonds3[nn_site[0]][nn_site[1]][nn_site[2]]
&= unbit[inverse_direction((i+j+k)%2,r)];
num_opened--;
break;
}
}
}
}
while ( num_opened < target_num_opened )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
if ( SITE_IS_OCCUPIED_3(i,j,k) )
{
// Look for a closed bond.
for ( r=0; r<coord_num; r++ )
{
if ( ! ( bonds3[i][j][k] & bit[r] ) )
{
// Bond is closed, so open it (assuming there is
// an occupied nn site in direction r).
if ( get_nn_site_in_dir(r,i,j,k) )
{
// Nearest neighbour site exists.
if ( SITE_IS_OCCUPIED_3(nn_site[0],nn_site[1],nn_site[2]) )
{
bonds3[i][j][k] |= bit[r];
// Mark the bond open at the nn site.
bonds3[nn_site[0]][nn_site[1]][nn_site[2]]
|= bit[inverse_direction((i+j+k)%2,r)];
num_opened++;
break;
}
}
}
}
}
}
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_diluted )
if ( !SITE_IS_OCCUPIED_4(i,j,k,l) )
continue;
get_nn_sites(i,j,k,l);
// See comments for 2d case.
bonds = 0;
for ( r=0; r<MAX_COORD_NUM; r++ )
{
bonds >>= 1;
if ( r >= coord_num )
continue;
if ( ( nn_i = nn_sites[r][0] ) == NON_EXISTENT )
continue;
nn_j = nn_sites[r][1];
nn_k = nn_sites[r][2];
nn_l = nn_sites[r][3];
if ( !SITE_IS_OCCUPIED_4(nn_i,nn_j,nn_k,nn_l) )
continue;
num_pairs_nn++;
if ( !bond_diluted )
{
bonds |= OCCUPIED;
num_opened++;
}
else
{
if ( SITE_IS_VISITED_4(nn_i,nn_j,nn_k,nn_l) )
{
nn_bonds = bonds4[nn_i][nn_j][nn_k][nn_k];
nn_r = inverse_direction((i+j+k+l)%2,r);
if ( nn_bonds & bit[nn_r] )
bonds |= OCCUPIED;
}
else
{
if ( PRNG1() < bond_concentration )
{
bonds |= OCCUPIED;
num_opened++;
}
}
}
}
bonds4[i][j][k][l] = bonds;
if ( bond_diluted )
MARK_SITE_AS_VISITED_4(i,j,k,l);
}
}
}
}
num_pairs_nn /= 2; // since each pair has been counted twice.
if ( !bond_diluted )
num_opened /= 2; // since each bond has been counted twice.
// If necessary, adjust bonds to get exact concentration specified.
target_num_opened = (int)(num_pairs_nn*bond_concentration);
while ( num_opened > target_num_opened )
{
i = (int)(size*PRNG1());
j = (int)(size*PRNG1());
k = (int)(size*PRNG1());
l = (int)(size*PRNG1());
if ( bonds4[i][j][k][l] )
{
// This spin has an open bond.
for ( r=0; r<coord_num; r++ )
{
if ( bonds4[i][j][k][l] & bit[r] )
{
// Bond is open, so close it.
bonds4[i][j][k][l] &= unbit[r];
// Mark the bond closed at the nn site.
get_nn_site_in_dir(r,i,j,k,l);
bonds4[nn_site[0]][nn_site[1]][nn_site[2]][nn_site[3]]
&= unbit[inverse_direction((i+j+k+l)%2,r)];
num_opened--;
break;
}
}
}
}
while ( num_opened < target_num_opened )
{
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) )
{
// Look for a closed bond.
for ( r=0; r<coord_num; r++ )
{
if ( ! ( bonds4[i][j][k][l] & bit[r] ) )
{
// Bond is closed, so open it (assuming there is
// an occupied nn site in direction r).
if ( get_nn_site_in_dir(r,i,j,k,l) )
{
// Nearest neighbour site exists.
if ( SITE_IS_OCCUPIED_4(nn_site[0],nn_site[1],nn_site[2],nn_site[3]) )
{
bonds4[i][j][k][l] |= bit[r];
// Mark the bond open at the nn site.
bonds4[nn_site[0]][nn_site[1]][nn_site[2]][nn_site[3]]
|= bit[inverse_direction((i+j+k+l)%2,r)];
num_opened++;
break;
}
}
}
}
}
}
break;
}
if ( bond_diluted )
mark_all_sites_unvisited();
return ( num_opened );
}