// SPINS.C
// Functions concerned with spins
// © 2021 Peter J. Meyer
#include "iss.h"
// Given that a site is occupied, what is its spin?
// (Use this if a site is known to be occupied
// and the spin value is required.)
//
// For the Ising model, bit 0 of sites[...] has
// 0 for UP_SPIN and 1 for DOWN_SPIN.
// For the q-state Potts model, the five lowest-order bits
// hold the spin value minus one.
// q-state Potts model spin values are in the range 1 through q.
//
// This is very similar to spin_at_site().
/*----------------------------------------*/
int spin_at_occupied_site(int i, int j, ...)
{
// int k, l;
unsigned char site;
if ( i < 0 ) // Occurs only when i == NON_EXISTENT.
return ( false );
switch ( dimensionality )
{
case 2:
site = sites2[i][j];
break;
case 3:
// k = *(((int *)(&j))+1);
// site = sites3[i][j][k];
site = sites3[i][j][*(((int *)(&j))+1)];
break;
case 4:
// k = *(((int *)(&j))+1);
// l = *(((int *)(&j))+2);
// site = sites4[i][j][k][l];
site = sites4[i][j][*(((int *)(&j))+1)][*(((int *)(&j))+2)];
break;
}
if ( model_is_ising )
return ( site & spin_value_mask ? DOWN_SPIN : UP_SPIN );
else // model_is_q_state_potts
return ( ( site & spin_value_mask ) + 1 );
}
// This function returns 0 if a site is not occupied
// otherwise it returns the spin value (Ising or q-state Potts).
// It should be used when both occupancy and spin value
// information is needed.
/*-------------------------------*/
int spin_at_site(int i, int j, ...)
{
// int k, l;
unsigned char site;
if ( i < 0 ) // Occurs only when i == NON_EXISTENT.
return ( 0 );
switch ( dimensionality )
{
case 2:
site = sites2[i][j];
break;
case 3:
// k = *(((int *)(&j))+1);
// site = sites3[i][j][k];
site = sites3[i][j][*(((int *)(&j))+1)];
break;
case 4:
// k = *(((int *)(&j))+1);
// l = *(((int *)(&j))+2);
// site = sites4[i][j][k][l];
site = sites4[i][j][*(((int *)(&j))+1)][*(((int *)(&j))+2)];
break;
}
if ( ! ( site & OCCUPIED ) )
return ( 0 );
// Site is occupied.
if ( model_is_ising )
return ( site & spin_value_mask ? DOWN_SPIN : UP_SPIN );
else // model_is_q_state_potts
return ( ( site & spin_value_mask ) + 1 );
}
// What is the sum of the nn spins of a spin?
//
// This information is obtained from the bond information for the site.
// Note that an occupied nn site is not a nn spin if no open bond exists.
//
// This is used only with the Ising model.
/*------------------------------*/
int nn_spin_sum(int i, int j, ...)
{
int k, l, r, spin_sum=0;
switch ( dimensionality )
{
case 2:
get_nn_sites(i,j);
for ( r=0; r<coord_num; r++ )
{
if ( bonds2[i][j] & bit[r] ) // if bond exists in direction r
spin_sum += spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1]);
}
break;
case 3:
k = *(((int *)(&j))+1);
get_nn_sites(i,j,k);
for ( r=0; r<coord_num; r++ )
{
if ( bonds3[i][j][k] & bit[r] ) // if bond exists in direction r
spin_sum += spin_at_occupied_site(nn_sites[r][0],
nn_sites[r][1],nn_sites[r][2]);
}
break;
case 4:
k = *(((int *)(&j))+1);
l = *(((int *)(&j))+2);
get_nn_sites(i,j,k,l);
for ( r=0; r<coord_num; r++ )
{
if ( bonds4[i][j][k][l] & bit[r] ) // if bond exists in direction r
spin_sum += spin_at_occupied_site(nn_sites[r][0],
nn_sites[r][1],nn_sites[r][2],nn_sites[r][3]);
}
break;
}
return ( spin_sum );
}
// What is the sum of the kronecker deltas
// of a spin with its nearest neighbour spins?
//
// This information is obtained from the bond information for the site.
// Note that an occupied nn site is not a nn spin if no open bond exists.
//
// This is used only with the q-state Potts model.
/*------------------------------------*/
int sum_kronecker_deltas(int spin_value,
int i, int j, ...)
{
int k, l, r, sum=0;
switch ( dimensionality )
{
case 2:
get_nn_sites(i,j);
for ( r=0; r<coord_num; r++ )
{
if ( bonds2[i][j] & bit[r] ) // if bond exists in direction r
sum += ( spin_value ==
spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1]) );
}
break;
case 3:
k = *(((int *)(&j))+1);
get_nn_sites(i,j,k);
for ( r=0; r<coord_num; r++ )
{
if ( bonds3[i][j][k] & bit[r] ) // if bond exists in direction r
sum += ( spin_value ==
spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1],nn_sites[r][2]) );
}
break;
case 4:
k = *(((int *)(&j))+1);
l = *(((int *)(&j))+2);
get_nn_sites(i,j,k,l);
for ( r=0; r<coord_num; r++ )
{
if ( bonds4[i][j][k][l] & bit[r] ) // if bond exists in direction r
sum += ( spin_value ==
spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1],
nn_sites[r][2],nn_sites[r][3]) );
}
break;
}
return ( sum );
}
// This calculates the absolute magnetization from scratch.
// distinguished_spin_value is normally 1,
// but when calculating the second moment of the magnetization
// in the q-state Potts model we call this function
// with distinguished_spin_value = 1, ..., q.
/*---------------------------------------------------*/
int get_abs_magnetization(int distinguished_spin_value)
{
int i, j, k, l; // array indices
int abs_mag=0;
int spin_value;
if ( model_is_ising && distinguished_spin_value != UP_SPIN )
{
printf("\nModel is Ising, but get_abs_magnetization() called"
"\nwith distinguished_spin_value != UP_SPIN.\n");
exit(1);
}
else if ( model_is_q_state_potts &&
( distinguished_spin_value < 1 || distinguished_spin_value > q ) )
{
printf("\nModel is q-state Potts and get_abs_magnetization()"
"\ncalled with invalid value for distinguished_spin_value.\n");
exit(1);
}
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( spin_value = spin_at_site(i,j) )
// spin_value is 0 if site is not occupied.
abs_mag += ( spin_value == distinguished_spin_value ? 1 : -1 );
}
}
break;
case 3:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
for ( k=0; k<size; k++ )
{
if ( spin_value = spin_at_site(i,j,k) )
abs_mag += ( spin_value == distinguished_spin_value ? 1 : -1 );
}
}
}
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 ( spin_value = spin_at_site(i,j,k,l) )
abs_mag +=
( spin_value == distinguished_spin_value ? 1 : -1 );
}
}
}
}
}
return ( abs_mag );
}
// abs_internal_energy is the sum of the products of the spins
// in each pair of nearest neighbour spins.
// This function presupposes that nn bonds have been set up.
/*-----------------------------*/
int get_abs_internal_energy(void)
{
int i,j;
// int k,l;
int spin_value,abs_int_energy=0;
if ( model_is_q_state_potts )
{
printf("\nMeasurement of internal energy for Potts model"
"\nis not currently supported. See SPINS.C and SPINFLIP.C.\n");
exit(1);
}
switch ( dimensionality )
{
case 2:
for ( i=0; i<size; i++ )
{
for ( j=0; j<size; j++ )
{
if ( spin_value = spin_at_site(i,j) )
// This returns 0 if the site is unoccupied.
{
if ( model_is_ising )
abs_int_energy -= spin_value*nn_spin_sum(i,j);
else // model is q-state Potts
// TENTATIVE!!
abs_int_energy -= spin_value*sum_kronecker_deltas(spin_value,i,j);
}
}
}
break;
case 3:
ADD_CODE // 3d
break;
case 4:
ADD_CODE // 4d
break;
}
return ( abs_int_energy/2 );
// Divide by 2 since each energy between pairs of spins is counted twice.
}
// For measurement of autocorrelation
// we make a copy of the spin assignment.
// initial_sites[..] also used to restore spin assignment
// when num_repetitions > 0
/*---------------------------*/
void copy_spin_assignment(void)
{
switch ( dimensionality )
{
case 2:
memcpy(&initial_sites2[0][0], &sites2[0][0], num_sites);
break;
case 3:
memcpy(&initial_sites3[0][0][0], &sites3[0][0][0], num_sites);
break;
case 4:
memcpy(&initial_sites4[0][0][0][0], &sites4[0][0][0][0], num_sites);
break;
}
}
/*------------------------------*/
void restore_spin_assignment(void)
{
switch ( dimensionality )
{
case 2:
memcpy(&sites2[0][0], &initial_sites2[0][0], num_sites);
break;
case 3:
memcpy(&sites3[0][0][0], &initial_sites3[0][0][0], num_sites);
break;
case 4:
memcpy(&sites4[0][0][0][0], &initial_sites4[0][0][0][0], num_sites);
break;
}
}
// This is the same as spin_at_occupied_site() except that
// initial_sites[...] is used rather than sites[...].
/*------------------------------------------------*/
int initial_spin_at_occupied_site(int i, int j, ...)
{
int k, l;
unsigned char site;
switch ( dimensionality )
{
case 2:
site = initial_sites2[i][j];
break;
case 3:
k = *(((int *)(&j))+1);
site = initial_sites3[i][j][k];
break;
case 4:
k = *(((int *)(&j))+1);
l = *(((int *)(&j))+2);
site = initial_sites4[i][j][k][l];
break;
}
if ( model_is_ising )
return ( site & spin_value_mask ? DOWN_SPIN : UP_SPIN );
else // model_is_q_state_potts
return ( ( site & spin_value_mask ) + 1 );
}